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

We describe an efficient numerical method for simulating the dynamics and steady states of collective spin systems in the presence of dephasing and decay. The method is based on the Schwinger boson representation of spin operators and uses an extension of the truncated Wigner approximation to map the exact open system dynamics onto stochastic differential equations for the corresponding phase space distribution. This approach is most effective in the limit of very large spin quantum numbers, where exact numerical simulations and other approximation methods are no longer applicable. We benchmark this numerical technique for known superradiant decay and spin-squeezing processes and illustrate its application for the simulation of non-equilibrium phase transitions in dissipative spin lattice models.


Introduction
Large ensembles of two-level systems that can be approximately modeled as a single collective spin are of interest in many areas of physics. In quantum optics, collective light-matter interaction effects can be understood from the analysis of the Dicke model [1], which describes the coupling of many two-level atoms to a common photonic mode. In the field of ultracold atoms, the evolution of Bose-Einstein condensates in double-well potentials can be mapped onto the motion of one large collective spin [2]. In nuclear and solid-state physics the Lipkin-Meshkov-Glick (LMG) model [3] is frequently used to investigate ferromagnetic phase transitions in systems with all-to-all interactions. In many situations one is interested in the dynamics or the steady states of those systems in the presence of dephasing and decay. For example, for magnetic field sensing, spin squeezing [4,5] and related metrological applications with large ensembles of atoms [6,7] the achievable sensitivities are primarily limited by such decoherence processes. But the coupling of large ensembles of two-level systems to a common environment can also lead to new physical phenomena, such as phase-locked condensates in equilibrium [8], or superradiant [1,9] and super-correlated [10] decay. From a theoretical and computational perspective, the primary interest in collective spin models arises from their permutational symmetry. This symmetry effectively reduces the full dimension of the Hilbert space of N TLS two-level systems, d = 2 N TLS , to the dimension d S = (2S + 1) of a spin S = N TLS /2 system. At the same time, the system can still exhibit interesting many-body effects and sharp phase transitions in the 'thermodynamic limit' S 1. For this reason, dissipative versions of the Dicke [11][12][13], LMG [14,15] and related collective spin models [16][17][18] play an important role in the analysis of non-equilibrium phase transitions in open quantum systems, since an exact numerical integration of the full master equation is still possible for moderately large ensembles. However, brute force numerical simulations are no longer feasible for atom numbers encountered in many of the actual experiments and, in general, for systems involving multiple collective spins. Such scenarios appear naturally in the presence of inhomogeneous couplings or frequencies [19,20] or in extended lattice systems [19,21], which are most relevant for analyzing critical phenomena. For closed systems, the coherent dynamics of spins with S 1 is typically well-described by meanfield theory, i.e., by evaluating the dynamics of the average spin vector S only. But this approach ignores important correlations, such as spin squeezing effects, and will in general provide a poor description of the actual state in the presence of dissipation. Here quantum fluctuations associated with incoherent processes can drive the system into highly mixed states [16,17,21], where the fluctuations of the spin components become comparable to their mean values. Thus, in order to accurately model such a behavior, it requires approximate numerical techniques, which take the effect of fluctuations into account, while still being able to simulate the dynamics of large collective spins efficiently.
In this paper we describe a broadly applicable stochastic method for simulating the dissipative dynamics of systems involving either a single or many coupled collective spins. The method relies on the Schwinger boson representation of spin systems and uses an extension of the truncated Wigner approximation (TWA) [22] to map the dynamics of those bosons onto a set of stochastic differential equations in phase space. For weakly interacting bosonic systems such phase space methods based on the TWA are well established and can be used, for example, to simulate dissipative bosonic lattice systems and non-equilibrium condensation phenomena [23][24][25][26][27]. The extension of these methods to spin systems via the Schwinger representation has previously been applied for simulating the coherent dynamics of lattices of spin-1/2 systems [28,29] and that of collective spins [22]. In the latter case also alternative methods based on the discrete TWA (DTWA) [30] are very efficient. In Ref. [31], Olsen et al. showed that a stochastic sampling of the positive-P representation of the Schwinger bosons can also be used to model the collective decay of an atomic ensemble. However, the practical applications of this approach are very limited since the stochastic trajectories derived from the positive-P distribution tend to diverge after rather short times [28,29,31,32] and, to our knowledge, this method has not been developed further. Here we show how this problem can be overcome for systems with S 1 by working with the Wigner function, but performing an additional positive diffusion approximation (PDA). As a result of this approximation, the stochastic equations in phase space are well-behaved for arbitrary times, which allows us to evaluate also the long-time dynamics and the steady states of dissipative spin systems that have been inaccessible so far. Using this truncated Wigner method for open quantum spins (TWOQS), we obtain an approximately linear scaling with the number of collective spins included, in any dimension and for arbitrary interaction patterns. In a recent work [21], we have already applied this method to identify novel PT -symmetry breaking transitions in the steady state of a one dimensional spin lattice with gain and loss. Here we provide a more detailed derivation of this simulation technique and discuss and benchmark its performance in terms of several explicit examples.
The structure of this paper is as follows: In Sec. 2 we first present a general outline of the method and explain how the original master equation can be mapped, under certain approximations, onto a set of stochastic differential equations. Then, in Sec. 3 we illustrate and benchmark this procedure by studying several model systems. We also go on and show how this technique can be applied to the simulation of non-equilibrium phase transitions in dissipative spin lattice models. Finally, in Sec. 4 we present our conclusions.

Outline of the Method
We are interested in the open system dynamics of i = 1, . . . , N coupled spin-S systems, which can be modeled by a master equation for the system density operator ρ, Here H is the many-body Hamiltonian describing the coherent evolution and the Lindblad superoperators, where D[c]ρ = 2cρc † − c † cρ − ρc † c, account for incoherent processes with jump operators c n and rates Γ n . In the following we assume that H and all c n can be written in terms of products of the collective spin operators S z i and S ± i = (S x i ± iS y i ), which obey the usual spin commutation relations, [S z i , S + j ] = δ ij S + i and [S + i , S − j ] = 2δ ij S z i . Equation (1) conserves the length of each individual spin, ∂ t S 2 i = 0, and therefore the dynamics of each subsystem can be restricted to a d S = (2S + 1) dimensional subspace. However, the dimension of the full density operator, d ρ = (d S ) 2N , still scales exponentially with the number of subsystems or lattices sites N . This scaling makes an exact numerical integration of Eq. (1) impossible when S or N are large. Here we introduce an approximate method, the TWOQS, to simulate such systems in the limit S 1, which only scales linearly with the system size N . The derivation of this method consists of four main steps: 1. The N spins are mapped onto a set of 2N bosonic modes using the Schwinger boson representation.
2. The master equation for the bosons is mapped onto an equivalent partial differential equation for the Wigner phase-space distribution.
3. We use the TWA and the PDA to obtain a Fokker-Planck equation (FPE) for the Wigner function with an (almost) positive diffusion matrix.
4. This FPE is mapped onto an equivalent set of stochastic Ito equations, which can be efficiently simulated numerically.
In the following, we first give a brief general outline of the individual steps in this derivation, while the application of this method for concrete examples is discussed in more detail in Sec. 3.

Bosonization
In a first step, we use the Schwinger boson representation to map each of the spins, S i , onto two independent bosonic modes, a i and b i , by identifying One can easily show that this transformation preserves all the spin commutation relations given above. For all models constructed from collective spin operators only, the total number of excitations at each site, The initial condition can then be chosen such that to simulate spins of different lengths. This is more useful than mapping each site to a single Holstein-Primakoff boson [33], since the transformation above does not involve any operator square roots, which can be numerically difficult to work with.

Phase Space Distributions
The main advantage of switching to a representation expressed in terms of bosonic modes is that the master equation, Eq. (1), can be mapped onto an equivalent partial differential equation for a class of phase space distributions, which contain the same information as the density operator [32]. We parameterize the set of distributions by the variable k = −1, 0, 1 and define where v = (a 1 , b 1 , a 2 , b 2 , . . . , a N , b N ) is a vector of all 2N bosonic annihilation operators and α and λ are vectors containing the same amount of complex numbers. When k = 0 this phase space distribution corresponds to the Wigner function, for k = 1 it is the Glauber-Sudarshan P -representation and when k = −1 we obtain the Husimi Q-function. We can use this definition to calculate what form each term in the master equation takes in the equation for F k [32]. For example, for a single mode one finds the mapping This translation lets us recast the master equation for ρ in the form of a partial differential equation for the phase space distribution, with L some linear differential operator that depends on the specific problem under consideration.

Truncated Wigner Approximation
The result in Eq. (9) is still exact and therefore in general not very useful. In particular, the differential operator L may contain third-or higher-order derivatives, which prevent an efficient stochastic sampling of F k . For example, for the coherent dynamics generated by the master equationρ = −iΩ[S 2 x , ρ], the corresponding partial differential equation To proceed we neglect all third-and higher-order derivatives, which in this example corresponds to omitting all terms in the second line of Eq. (10). This approximation is just the usual TWA [22] applied to arbitrary distribution functions. For spin systems we expect this approximation to become accurate in the limit of large S, since terms proportional to αF k or βF k scale as ∼ √ S compared to derivatives such as ∂F k /∂α ∼ O(1). After performing the TWA we obtain a FPE of the form with a drift matrix A and a diffusion matrix D. Here we have assumed Einstein's sum convention, where the indices i and j run over the 4N components of the vector

Positive Diffusion Approximation
For stochastic simulations, performing the TWA is not enough since in general the diffusion matrix D obtained in this way is not positive semi-definite. This can already be seen from the underlined terms in Eq. (10). Similarly, we find that an incoherent decay process,ρ = ΓD[S − ]ρ, is mapped under the TWA onto the FPE and there are again second-order derivatives that can lead to negative diffusion rates. Thus, in a second step we perform a PDA by neglecting some of these diffusion terms. In the two examples above this approximation amounts to omitting all the underlined terms in Eq. (10) and Eq. (12), while keeping the diffusion terms in the second line of Eq. (12). This choice cannot be justified by simple scaling arguments and in Sec. 3 we discuss and verify the applicability of this approximation in terms of several explicit examples. In general, the PDA can be motivated by the fact that it eliminates the dominating negative contributions to D, while conserving the total spin S and leaving the equations of motion for the mean values S k unaffected. The price we pay for this last requirement is that for k = 0 the resulting diffusion matrix can become negative for certain values of α. However, the corrections scale as ∼ 1/S compared to other terms and for S 1 the residual negative contributions do not affect considerably the stochastic sampling of trajectories in actual simulations.
Before we proceed let us remark that the problem of non-positivity can also be overcome by working with a positive-P representation, where α i and α * i are replaced by a pair of independent complex variables [31,32]. In this case, a positive semi-definite diffusion matrix can be obtained for this larger set of variables without neglecting any terms. However, it is known that the resulting stochastic equations are often not well-behaved [32]. In particular, the appearance of "spikes", where individual trajectories diverge at a finite time [28,29,32], often prevents the simulation of the long-time behavior of a system or its steady state.

Stochastic Simulations
After applying the TWA and the PDA, we end up with a FPE with an (almost) positive semidefinite diffusion matrix D. This FPE can be mapped onto an equivalent set of stochastic (Ito) differential equations [34], where dW i are real-valued independent Wiener processes with dW i dW j = δ ij dt and B( x) is the factorized diffusion matrix with B( x)B( x) † = D( x). This set of equations can be efficiently simulated with the Euler-Maruyama method [34]. This means that we do not calculate the full probability distribution, but instead obtain the required expectation values by averaging over n traj trajectories of these stochastic equations. Note that for a closed system, where all second-and higher-order derivatives have been neglected, the amplitudes α i evolve according to the mean-field equations of motion,ẋ consistent with the usual applications of the TWA [22]. In the presence of dephasing or decay our approach accounts for the corresponding damping terms and the associated amount of quantum fluctuations in a consistent manner. For sufficiently many trajectories and with initial values sampled according to the distribution F k ( α, t = 0), these stochastic averages provide accurate approximations of the corresponding quantum mechanical expectation values where, depending on the chosen distribution function, . . . P |Q|W denotes the normallyordered, anti-normally-ordered or symmetrically-ordered expectation value. All expectation values of the original spin system can then be obtained using the relations in Eq. (2).

Initial conditions
In many situations of interest the initial state can be chosen as a fully polarized state with S z i = −S at each site. This corresponds to a state where one of the two Schwinger bosons is prepared in the vacuum state |0 , the other one in the Fock state |2S . For k = −1 this state is described by the Q-function which is positive everywhere and can be used as an initial probability distribution for the trajectories. For k = 1 and k = 0 the corresponding P -and Wigner distributions for Fock states are singluar or have negative values. It is thus necessary to approximate the initial state by replacing the Fock state |2S by a coherent state with the same mean amplitude. The corresponding initial conditions are then given by and respectively. This approximation introduces an uncertainty in the spin quantum number S, which, however, scales only with √ S and becomes negligible in the limit of interest, S 1. In order to initialize the system in an arbitrary spin coherent state |θ, φ on the Bloch sphere we can simply rotate this state by the angle θ around the y-axis and φ around the z-axis. This amounts to replacing α and β by the rotated amplitudes α = e iφ (cos(θ/2)α − sin(θ/2)β), i.e., W θ,φ (α, β) = W 0 (α,β).
Up to now we have kept our analysis completely general and derived all the results for arbitrary distribution functions F k ( α). This raises the question of which distribution to choose in an actual simulation? 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 given in Eq. (13) when k=1. The Q-distribution (k = −1) has the obvious advantage that it can represent spin states with a well-defined spin quantum number, i.e., there is no need to approximate the initial state. Further, as can be seen from Eq. (12), after the PDA the diffusion matrix for the Q-distribution is strictly positive semi-definite. However, it turns out that for models that include (S x ) 2 or similar interaction terms in the Hamiltonian, performing the PDA eliminates relevant contributions to the coherent dynamics. As can be seen from Eq. (10), this is not the case for the Wigner distribution (k = 0), since there are no second-order derivatives in the Hamiltonian dynamics and the PDA only affects incoherent processes. This is true for all quadratic coupling terms in the Hamiltonian ∼ S ν i S µ j (ν, µ = z, ±), which already includes the most common types of spin-spin interactions. Therefore, while below we will also discuss several basic examples where the P -or the Q-distribution yield equally accurate results, we find that for generic interacting systems it is necessary to work with the Wigner function, which reproduces most accurately the Hamiltonian part of the dynamics.

Examples and Applications
In this section we will present several explicit examples, to show how our method can be applied to simulate some of the most frequently encountered interactions and decoherence processes. To do so we will mainly focus on systems with a single collective spin, where all the results can still be compared with exact numerical results. This will allow us to test the validity of the approximations described above and make a comparison of how the different phase space representations perform under different circumstances. In Sec. 3.5 we will then extend these results and discuss the simulation of a whole chain of collective spins, for which exact methods are no longer available.

Spontaneous Emission
As a first example we consider the collective decay of a large ensemble of two-level systems, which can be described by the master equatioṅ Note that here we rescale the emission rate by a factor 2S in order to obtain the same time scale for the dynamics for different values of S. After performing the TWA the resulting FPE for this model is already given in Eq. (12) above. The PDA then corresponds to neglecting the underlined term in this equation, after which we can map the FPE onto the following set of stochastic Ito equations where the dW n are real-valued and independent Wiener processes with dW n dW m = δ nm dt.
In Fig. 1 we plot the outcome of a stochastic simulation of this coupled set of equations for k = 0, ±1 and for two different spin quantum numbers, S = 10 and S = 100. In these examples it is assumed that the spin is initially prepared in the maximally excited state with S z |S = S|S , which we represent by initial distributions as given in Sec. 2.5.1. For the considered values of S we can also solve the full master equation exactly and use these results to benchmark our approximate approach. We find that for about n traj = 1000 trajectories the TWOQS reproduces very accurately the superradiant decay of a large ensemble, with higher accuracy for larger values of S. For this example we find almost no visible differences between the three different distribution functions. However, a closer inspection shows that in the case of the Wigner function (k = 0), the square root in Eq. (23) can become negative for some trajectories. This becomes a crucial problem for very small values of S and restricts stimulations to short integration times, since at longer times these unphysical trajectories can dominate the dynamics. For larger spins, this error is suppressed by 1/S and becomes a negligible effect for S 100, as shown in Fig. 1. In a simulation, possible errors arising from the negative diffusion term can be easily tracked by monitoring the change of the total spin, i.e., |α| 2 + |β 2 | , over time. This example illustrates that even for the Wigner function, residual negative diffusion terms are not a practical limitation for simulating dissipative processes in collective spin systems when S is large. Instead, when using the exact positive P -representation [31], the same simulation would be limited to times of about t Γ −1 , before the appearance of spikes prevents any converging results. Note that the same conclusions also apply to master equations with a gain term, D[S + ], which can be described by simply exchanging the two bosonic modes, i.e., α ↔ β in Eqs. (22) and (23).

Dephasing
We now proceed with the derivation of the stochastic equations of motion for a collective spin which is subject to dephasing. In the absence of any other interactions, dephasing can be described by the master equationρ The bosonized form of this equation is obtained by substituting S z → (a † a − b † b)/2 and under the TWA the resulting FPE reads t). (25) Although also in this case there are second-order derivatives with negative prefactors, a straight-forward diagonalization of the diffusion matrix shows that D(α, β) is already positive semi-definite for all α and β. In this case the PDA is obsolete and we can factorize the diffusion matrix as D(α, β) = B(α, β)B(α, β) † , where Note that this factorization is not unique, but with the current choice we obtain a very simple and symmetric form for the stochastic equations, where dW is a single real-valued Wiener processes. These equations are independent of k and there is no preferred phase space distribution to simulate dephasing. In the example plotted in Fig. 2, which shows the dephasing of a spin that is initially polarized along the x direction, the stochastic averages for all distributions agree within the statistical errors with the exact dynamics, keeping in mind that for k = 1 and k = 0 the initial distributions are only approximate.

Dynamics and Steady States of Driven Spin Systems
We now consider slightly more complicated models in which there is an interplay between coherent driving and incoherent decay. The simplest model in this class is that of a collective spin driven by a transverse field of strength Ω and including a collective decay with rate Γ. The corresponding master equations readṡ with a Hamiltonian H D = ΩS x .

Transient dynamics
In Fig. 3 we show again a comparison between the TWOQS and the exact numerical simulations of this master equation for all three distribution functions and for the spin quantum numbers S = 10 and S = 100. For S = 10, we find clearly visible deviations from the exact oscillations, which can in part be traced back to the approximation we made in the initial condition (see Sec. 2.5.1). For this reason, sampling of the Q-function is most accurate in this situation. However, these deviations become negligible when we consider higher spins and already for S = 100 all distribution functions reproduce very precisely the exact spin dynamics over many oscillation periods.

Steady states
A specific interest in the model given in Eq. (29) arises from the fact that it exhibits a nonequilibrium phase transition at a driving strength of Ω = Γ [17,35,36]. At this point the steady state of this system changes from a spin coherent state on the lower half of the Bloch sphere to a highly mixed state with S z = 0.
From the analysis of coherent bosonic or spin systems it is known that the TWA often leads to inaccurate results for long simulation times [22]. The same problem is encountered when the TWOQS is used to simulate, for example, the oscillations shown in Fig. 3 for much longer times. However, the timescale beyond which significant errors occur increases with S and for many practical applications the system reaches a steady state before problems arise. This is demonstrated in Fig. 4, where we use our stochastic approach to simulate the master equation for a driven spin with S = 1000 up to a time t = 50Γ −1 . Note that for Eq. (29) there still exists an analytic solution for the steady state [17,35,36], which allows us to compare these simulations with the exact results for the mean values and the fluctuations of the spin components.
In the polarized phase, Ω < Γ, we find that both the mean values as well as the fluctuations of all spin components agree almost perfectly with the exact results. For the considered example of S = 1000 there are still some visible differences for the predicted spin fluctuations at and above the transition point, Ω/Γ = 1. However, as shown in the inset of Fig. 4 [17,35,36], while the crosses, diamonds and circles are obtained from a stochastic sampling of the P -, the Q-and the Wigner-distribution. The inset in (a) shows the simulations of the Wigner distribution for even larger spin numbers S around the transition point Ω/Γ = 1. The steady state was obtained by time averaging after t = 40Γ −1 for another period of ∆t = 10Γ −1 and for n traj = 2500. the non-analyticity at the phase transition point becomes more pronounced and closer to the exact result by increasing the spin quantum number S. We emphasize that in the whole mixed phase, Ω/Γ ≥ 1, the Liouvillian gap of the considered model, i.e., the smallest decay rate in the problem, scales as ∼ 1/S. This means that in the mixed phase this system is particularly challenging to simulate and oscillations around the steady state can persist for very long times. Nevertheless, we see that by simply approximating the steady state at a fixed time t = 40Γ −1 by an average over a time span of ∆t = 10Γ −1 , all the essential features of the model are already rather accurately reproduced. In particular, for Ω Γ, all the fluctuations are around (S k ) 2 ∼ S 2 /3, indicating that the system is close to a fully mixed state. This and other examples show that by using the TWOQS it is possible to access the steady states of driven-dissipative collective spin models.

Spin Squeezing
Spin squeezing is an important non-classical effect in quantum metrology, which reduces the variance of one spin component below the value of S/2 obtained for N TLS independent twolevel systems. In the presence of collective decay and dephasing, the effect of spin squeezing can be described by the master equatioṅ where the Hamiltonian term ∼ S 2 x has already been discussed as an example in Sec. 2. Therefore, under the TWA and the PDA we obtain the stochastic equations, where the last two terms in each line account for the decay and dephasing processes described by Eqs. (22)- (23) and Eqs. (27)-(28), respectively. In Fig. 5, we use the approximate stochastic equations to simulate the spin squeezing parameter ξ as a function of time. For a state pointing in the z-direction this parameter is defined as [4] where (∆S φ ) 2 = (S φ ) 2 − S φ 2 and S φ = cos(φ)S x + sin(φ)S y . Note that a squeezing parameter below unity, ξ < 1, requires a finite amount of entanglement between the two-level systems [37]. Compared to all the previous examples, we now see a clear difference between the results obtained for different distributions. For k = 1 the value of the squeezing parameter is ξ ≥ 1 for all times, since squeezed states can only be represented by a non-positive P -distribution. Therefore, these results have not been included in Fig. 5. For the Q-distribution we obtain a finite amount of squeezing, but the predicted values for ξ do not match at all the exact results. This discrepancy can be traced back to the fact that in the PDA we neglect essential contributions to the coherent dynamics, which appear whenever there are spin-spin interactions. Therefore, in such cases neither the P -nor the Q-distribution give reliable predictions.
For simulations based on the Wigner function, we find very accurate results for ξ at short times, but considerable deviations from the exact behavior for longer simulations when Γ is small. This is consistent with the general observation that the TWA is not well suited to simulate coherent dynamics over longer times. However, these discrepancies are significantly reduced for larger dissipation rates and for larger spin quantum numbers. Importantly, Fig. 5 shows that already for S = 100 the dissipative evolution into an entangled quantum state with ξ 2 ≈ 0.05 − 0.5 can be accurately simulated with our method. As further demonstrated in the lower two panels of Fig. 5, this level of accuracy is sufficient to predict optimal squeezing parameters in open quantum systems, as relevant for metrological applications. Very similar conclusions can be obtained from the investigation of squeezing in the presence of dephasing, as summarized in Fig. 6. In general we find that dephasing processes are more accurately captured by our method than decay.

Spin Chains
In all the examples so far we have considered the dynamics of a single spin, where for S ≈ 100 the full master equation can still be solved exactly. This is no longer possible for systems involving N 2 collective spins, while the TWOQS scales only linearly with N . This feature becomes highly relevant, for example, for the study of non-equilibrium magnetic phases in driven-dissipative spin chains. In this context, one typically considers generic Heisenberg models of the form [38,39] where in addition each spin is subject to decay. Thus, the master equation for this system readsρ whereJ k = J k /(2S) andΓ = Γ/(2S) are the rescaled coupling strengths and the rescaled dissipation rate for general spin-S systems. For S = 1/2, Eq. (35) can still be simulated for large 1D chains using numerical techniques based on matrix-product operators [39]. However, in this case one does not observe any sharp phase transitions for finite Γ, while the reliability and applicability of related techniques for 2D systems are still under investigation [40][41][42]. Both in 1D and 2D, such tensor network methods have very unfavorable scaling for larger S. The current method allows us to address the limit S 1, where already in 1D distinct non-equilibrium phases and sharp transitions between them are expected. In a previous work [21] we have already applied this approach to study PT -symmetry breaking transitions in spin chains with both gain and loss, which can be mapped back onto a loss-only model withJ x = −J y andJ z = 0. Here we outline the implementation of this method for the general Heisenberg model in Eq. (34). Since we are dealing with in interacting spin system we must use the Wigner function, i.e., k = 0. After carrying out the general procedure described in Sec. 2 we obtain the stochastic equations and Depending on the relations between all the coupling parameters and the dissipation rate, the model in Eq. (34) exhibits many different stationary phases, which have been analyzed in Ref. [38] using mean-field theory. As a proof-of-concept demonstration of the TWOQS we consider here the case J z = 0. Then for J x J y > −Γ 2 the steady state of the system is the fully polarized state along the z-direction and we can use a Holstein-Primakoff approximation to study the fluctuations around this state, similar to the analysis in [21,38]. Beyond the transition point, e.g. for J x > 0 and J y < −Γ 2 /J x , we expect a strongly mixed phase, but in this regime mean-field theory and linearization techniques are no longer applicable. In Fig. 7 we show the results of a stochastic simulation of a spin chain with N = 100 sites and S = 5000. This simulation confirms that in the limit of large S there is a non-equilibrium phase transition between a polarized and a highly mixed phase, even in 1D. At the transition point the mean value of S z and the fluctuations of all spin components exhibit a sharp jump and spin-spin correlations along the chain diverge. In the polarized phase we can still use the Holstein-Primakoff approximation to benchmark the simulations also in this extended chain and we find almost perfect agreement. Importantly, the TWOQS also allows us to explore the non-polarized phase, where the strong fluctuations cannot be captured by a Holstein-Primakoff or mean-field approximation. While a detailed analysis of this phase is outside the scope of this work, we find many similarities with the pseudo PT -symmetric phase described in Ref. [21], where further discussions about its physical properties can be found.

Conclusion
In summary, we have introduced a new numerical method for simulating the dissipative dynamics of collective spin systems. This method works best in the limit of large spin quantum numbers, where exact simulations are no longer possible. At the same time the TWOQS goes beyond mean-field theory by taking the most relevant quantum fluctuations associated with dephasing and decay processes into account. A crucial step in the derivation of the stochastic differential equations is the PDA, which enforces the positivity of the diffusion terms.
Although seemingly a very crude approximation, it does not affect the accuracy of actual simulations for large S and allows us to access the long-time dynamics and steady states of open spins systems. This was not possible using previous approaches based on the otherwise more accurate positive P -distribution.
We have illustrated and benchmarked the application of this method for various spin models with dephasing and decay. Since the accuracy of the method improves with increasing S and only scales linearly with the number of spins, these simulations can be directly applied for atomic ensembles with sizes encountered in real experiments or be extended to simulate dissipative spin models in two or even three dimensional lattices. Finally, this technique can be readily combined with existing TWA simulations for bosonic systems and therefore be applied as well for simulating Dicke-type models, where collective spins are coupled to single or multiple bosonic modes.