Phase-Space Methods for Simulating the Dissipative Many-Body Dynamics of Collective Spin Systems

First of all, we would like to thank the reviewer for spending the time to carefully read our manuscript for providing us with constructive feedback. We provide detailed replies to the individual reports below:

Report 2:
Referee: “1. The survey on the existing literature for phase space methods for driven-dissipative systems is basically non-existent. Works like [Phys. Rev. B 72, 125335; Phys. Rev. X 5, 041028; Phys. Rev. A 97, 013853] should be referenced somewhere.”
Reply: The application of phase space methods for simulating driven-dissipative systems is a well-established technique, which is discussed in many textbooks and has been applied in numerous papers. The purpose of this work is not to describe this method per se, but to show how it can be extended to dissipative spin systems in a practical manner, which was not known before. However, we agree with the referee that a few examples for phase space methods for bosonic many-body systems could be of interest to the reader. We have added the suggested references and few more in the introduction.

Referee: “2. The notion of the conserved quantitites of Eq. (1) is not entirely clear. For S = 1/2, the conservation of S^2 appears to be trivial, but that does not help toward reducing the complexity of the problem.”
Reply: Let us emphasize that in this work we are interested in models where each lattice site is represented by a collective spin with S>>1. In practice, each collective spin S_i might be formed by an ensemble of many two-level systems, such that the full Hilbertspace scales exponentially in S. Therefore, it is important that in Eq. (1) there appear only collective spin operators such that the Hilbertspace of each site can be restricted to (2S+1) states. For S=1/2 there is no gain, but this is also not of interest for our method.
For a lattice system the total Hilbertspace dimension still scales exponentially in the number of lattice sites N and that is why approximate simulation methods are required.

Referee: “3. There is a typo in the 3rd line of p.8, where the Q distribution should refer to k = -1.
Reply: We thank the referee for pointing out this typo, which we have corrected.

Referee: “4. It is not entirely clear what the authors mean by the P distribution being incapable of representing squeezed states. As the P distribution is an integral transform of the density matrix, the P distribution of squeezed states clearly exists, even if only in a distributional sense. It may be possible that the authors cannot represent such distributions within their numerical work, but then they should say so.”
Reply: To be more precise, we have rewritten this statement in Sec. 2.6 in the following way: “It is well-known that squeezed states, which appear commonly in interacting spin systems, cannot be represented by a positive and non-singular P-distribution and thus cannot be simulated via the stochastic equations in Eq. (13) when k=1.” A similar clarification has been added in Sec. 3.4.

Referee: “5. The claim of the method being able to accurately predict the steady state of collective spin models is not justified. Most importantly, the benchmarking against the exact solution seems to fail at the phase transition, i.e., where the dynamics is most interesting. In particularly, the non-analytic behavior of the order parameter does not seem to get captured correctly. Clearly, this point deserves more attention and discussion.”
Reply: In Figure 4 we show that our method can reproduce the steady state of a collective spin model very accurately over the whole parameter range. For \Omega<\Gamma the results are essentially exact. At and above the critical point the predictions are qualitatively correct but there are some visible quantitative differences. Here one should keep in mind that at the critical point and throughout the mixed phase, the Liouvillian gap scales as ~1/S, meaning that the model considered in this plot is particularly challenging and such deviations are not found in non-critical models or parameter settings.
In the revised version of the manuscript we have added an inset in Fig 4(a) which shows that our method improves for larger S, in which case also the non-analytic behavior at the transition point is much more pronounced. Also, for the lattice system in Sec. 3.5 we clearly see the non-analytic phase transition, diverging correlation lengths etc., which exactly match the Holstein-Primakoff predictions even very close to the phase transition point. In this sense we believe that our claim of “accurately predicting the steady state” is justified.

Referee: “6. It is not clear what the authors what to say when referring to mean-field theory and how their work constitute a beyond mean-field result. This is particularly striking when comparing with Ref. [33], because the mean-field theory employed there can certainly be applied to highly mixed states (whether it is accurate is another question).”
Reply: We do not claim that mean-field theory cannot be applied to the models considered in this work, but it is obvious that a mean-field model, which only accounts for the average values of the spin components, cannot accurately represent states that are dominated by fluctuations. Our method accounts for (quantum) fluctuations and therefore, by definition, goes beyond mean-field approaches. This allows us to predict, for example, spin squeezing effects or the spin variances in the highly mixed phase in Fig. 4. As discussed in more detail in [21], the mean field theory in Ref. [33] (Ref. [38] in the revised version) is also not able to predict the properties of mixed phases in spin lattice models, even when S>>1.
To avoid any confusion, we have clarified the meaning of “beyond mean-field” in the introduction.

Referee: ”7. Contrary to what the authors say, there are works on the dissipative Heisenberg model in 2D using tensor network methods [Nature Commun. 8, 1291; arXiv:2012.03095].”
Reply: We do not claim that there are no works on 2D dissipative Heisenberg models, but that tensor network methods scale unfavorably in 2D, in particular for higher spin quantum numbers. Since the first submission of our work, two relevant papers on this topic appeared on the arXiv [arXiv:2012.03095, arXiv:2012.12233], which we have cited in the revised version together with the previous result published in Nature Commun. 8, 1291.

Report 1:
Referee: “Below eq. (37), the condition for the appearance of the polarized state and the mixed state is the same. I assume in the first equation it should be > instead of < (i.e., large dissipation)."
Reply: We thank the referee for pointing out this typo, which we have corrected.