Collective Radiative Interactions in the Discrete Truncated Wigner Approximation

Interfaces of light and matter serve as a platform for exciting many-body physics and photonic quantum technologies. Due to the recent experimental realization of atomic arrays at sub-wavelength spacings, collective interaction effects such as superradiance have regained substantial interest. Their analytical and numerical treatment is however quite challenging. Here we develop a semiclassical approach to this problem that allows to describe the coherent and dissipative many-body dynamics of interacting spins while taking into account lowest-order quantum fluctuations. For this purpose we extend the discrete truncated Wigner approximation, originally developed for unitarily coupled spins, to include collective, dissipative spin processes by means of truncated correspondence rules. This maps the dynamics of the atomic ensemble onto a set of semiclassical, numerically inexpensive stochastic differential equations. We benchmark our method with exact results for the case of Dicke decay, which shows excellent agreement. For small arrays we compare to exact simulations and a second order cumulant expansion, again showing good agreement at early times and at moderate to strong driving. We conclude by studying the radiative properties of a spatially extended three-dimensional, coherently driven gas and compare the coherence of the emitted light to experimental results.


Introduction
The accurate description of non-equilibrium dynamics of interacting quantum spin systems is one of the major challenges of many-body theory.At the same time it is of central importance in many areas of physics.A prime example is the collective interaction of two-level atoms with the quantized electromagnetic field, which after integrating out the radiation field can be mapped onto collectively coupled spins with long-range interactions and dissipation.Collective light-matter interactions have been a central problem in quantum optics starting from the early work of Dicke [1].Dicke showed that an ensemble of closely spaced two-level quantum emitters can display intriguing collective effects in the emission of light such as suband superradiance [2,3] observed in a number of experiments [4][5][6].This collective coupling between light and atoms has recently regained substantial interest as it is at the heart of many photonic quantum technologies [7].Collective light-atom couplings are for example the basis of ensemble-based quantum memories for photons [8][9][10], quantum repeaters [11], and many concepts for realizing strongly interacting photons [12][13][14][15].Here the interplay of the nonlinear atomic response and quantum entanglement results in rich coherent many-body dynamics.
A comprehensive theoretical treatment of the collective interaction of light with quantum emitters is however only simple if the spatial extension of the emitters can be neglected as in the case of the Dicke model or in cavity QED.Spatially extended systems can only be described by solving the master equation, e.g. by Monte-Carlo Wave Function (MCWF) simulations [16], if the number of excitations is small or for small ensemble sizes.Numerical techniques based on matrix product states [17], which have proven to be extremely powerful for one-dimensional systems with short-range couplings are usually not appropriate in higher spatial dimensions and for long-range couplings.Likewise a classical treatment of collective phenomena in terms of Maxwell-Bloch equations does not capture the buildup of quantum correlations between the atoms.While some universal features of superradiance can be predicted for spatially extended systems without involved numerics [18,19], there is for example no simple access to the timing and intensity of superradiant bursts.Expanding on the classical mean-field description in terms of Maxwell-Bloch equations, cumulant expansion techniques have been employed to account for correlations [20][21][22], but generally require higher order expansions for accurate predictions.Cumulant expansions are furthermore often ill-controlled and can suffer from intrinsic instabilities.Moreover their numerical complexity grows as a power law of the expansion order n, i.e. scales as N n , where N is the number of spins, making them computationally expensive.The same holds for non-equilibrium Greens function approaches such as the one employed in [23].
In the present paper we propose an alternative, semiclassical approach that allows to describe the coherent and dissipative many-body dynamics of interacting spins, taking into ac-count lowest-order quantum fluctuations.Our approach is inspired by the success of the discrete truncated Wigner approximation (DTWA) for the treatment of unitarily interacting spin systems [24], which has recently been extended to include single-particle dissipation [25][26][27].Within the truncated Wigner approximation the dissipative many-body dynamics of spins is mapped to a generalized diffusion problem of the Wigner quasi-probability distribution in phase space.The exact relation between the dynamics of the many-body density matrix in Hilbert space and the Wigner function in phase space is given by correspondence rules, which lead to higher-order partial differential equations for the Wigner function.These are in general intractable without further approximations.A very successful approximation applicable to unitarily coupled spins is the DTWA, which can be extended to include single particle decay and dephasing [25].The approach of Ref. [25] is however not useful for collective dissipative processes such as superradiance.We here pursue a different route.We propose approximate correspondence rules which lead to Fokker-Planck type equations of motion for the Wigner quasi distribution equivalent to a set of coupled stochastic differential equations (SDEs) for the spin orientations.Since the number of these equations scales linearly in the number of spins, the solution is numerically inexpensive and allows investigating system sizes much larger than in other semiclassical approaches.In the truncated Wigner approximation quantum fluctuations are taken into account to lowest order by nondeterministic initial conditions and by collectively coupling the spins to white noise processes, which generate (weak) entanglement between the spins.
Our paper is organized as follows: In Sect. 2 we give a compact summary of the continuous Wigner phase space representation of an ensemble of two-level systems (spins).We introduce a truncated Wigner Approximation for spin ensembles with collective couplings in Sect.3. In particular we propose and motivate approximate correspondence rules and discuss general conditions for their validity.The main application of our methods are collective lightmatter couplings in free space, which we will introduce in Sect. 4. In Sect. 5 we benchmark our method for the Dicke decay, i.e. the collective emission of light from a tightly localized ensemble of two-level atoms, for which the full master equation can be solved exactly.We find excellent agreement and give a physical interpretation of the emerging collective response within the semiclassical approximation.We then study collective light-matter phenomena in spatially extended systems in Sect.6, where the full dynamics can no longer be described exactly.We consider superradiance from regular arrays of atoms and an elongated cloud of coherently driven atoms.Finally Sect.7 summarizes the results and gives an outlook to future work.

Wigner Representation for Spins
An approach widely used in quantum optics to describe the dynamics of interacting, drivendissipative many-body systems beyond the mean-field level is the truncated Wigner approximation (TWA) [28][29][30][31][32].It describes interactions on a mean-field level but allows taking both thermal and leading-order quantum fluctuations into account by averaging over nondeterministic initial conditions and by coupling to stochastic noise sources.In the following we will give a compact summary of the Wigner representation of an ensemble of two-level systems or spins, but refer to Refs.[33,34] for a more general introduction to phase-space representations in quantum mechanics.We formulate the TWA by studying the correspondence rules [35], which translate the action of an operator in Hilbert space to a differential operator in phase space, and show that they have a simple asymptotic limit for collective processes.

Wigner representation of two-level systems
The connection between Hilbert space and Wigner phase space, spanned by some c-number variables Ω is given by the representation of an operator Ô in terms of a complex function W Ô(Ω), called the Weyl symbol Here ∆(Ω) is the so-called phase point operator or Wigner kernel.Inversely the Weyl symbol can be expressed explicitly in terms of the operator by Of particular interest is the Weyl symbol W ρ (Ω) of the density operator ρ, which is called the Wigner function or Wigner (quasi-probability) distribution.The latter notion is due to the fact that W ρ (Ω) ∈ and but can be negative.
Originally formulated for continuous degrees of freedom the concept of phase space representations can be extended to systems with finite-dimensional Hilbert spaces [36] such as spin- 1  2 systems.There is however some freedom in choosing the phase point operators.A specific discrete representation has been introduced by Wooters in [36], which is the foundation of the discrete truncated Wigner approximation [24].Here we adopt however a different, continuous representation of spin-1/2 states ρ through rotations of the particular discrete phase point operator ∆0 = 1 2 ( ˆ 2 + 3 σz ) which was shown in [25] to be more appropriate to describe dissipative spin systems.Here U(θ , φ, ψ) = e −i σz φ/2 e −i σ y θ /2 e −i σz ψ/2 , are the SU(2) rotation operators with Euler angles (θ , φ, ψ), which gives Note that in [25] we have used a slightly different definition of the Wigner kernel that is obtained by letting θ → π − θ and φ → −φ.Based on the same work, we have shown that, using a gauge freedom, the most relevant states | ↓〉 and | ↑〉 can be represented by simple, positive Wigner functions.Namely which are straightforwardly verified by substituting them into Eq.(1).The above Wigner function can be sampled by a fixed θ = θ ψ and drawing 0 ≤ φ < 2π from a uniformly random distribution.We note furthermore that the vector s (θ , φ) appearing in Eq. ( 5) is just the Weyl symbol of the Pauli spin-matrices Similarly, the Weyl symbols of the creation-and annihilation operators σ± = ( σx ± i σ y )/2 are given by The kernel in Eq. ( 5) can easily be extended to a system of N spin-1/2 systems via s → s n and Ω → Ω = {Ω n }, with j = 1, 2, . . ., N labelling the spins.

Time evolution of the Wigner function
Our goal is to find an approximate solution of the master equation of many-body spin systems where the many-body Hamiltonian Ĥ and the Lindblad operators Lµ , describing Markovian dissipative processes, are some functions of the spin-1/2 operators σµ j .Generically Ĥ and/or the Lµ describe interactions between spins which are higher dimensional, i.e. have couplings that cannot be reduced to a one-dimensional topology.The latter excludes in general efficient descriptions in terms of matrix product states [17].
To develop an approximate, semiclassical approach we need to translate the master equation of the density operator ρ into an equation of motion for the Wigner function W ρ (Ω).As the terms on the right hand side of Eq. ( 9) can be decomposed into products of spin operators and the density operator, this requires expressing the Weyl symbol of a composition of operators as emerging on the r.h.s. of ( 9), e.g.Ĥ ρ → W Ĥ ρ in terms of the individual symbols W Ĥ and W ρ .In phase space the Weyl symbol of a product does not correspond to a simple multiplication W ÂB ̸ = W Â • W B of the scalar functions.Instead, the composition is given by the Moyal product or star product which also has a differential form [33,37].The so called correspondence rules allow us to express the star product of Weyl symbols involving a spin operator, such as W σµ j ⋆ W Â, as differential operators acting on W Â. With these rules we can iteratively translate compositions of operators as they appear in the master equation of ρ into a partial differential equation for the Wigner function W ρ .
A more direct approach for deriving phase space equations that we have recently considered [25] is based on a simple observation: For a continuous phase space representation of a single spin, generated by the kernel given in Eq. ( 5), the matrices ∆, ∂ θ ∆, ∂ φ ∆ and ∂ 2 φ ∆ span the Hilbert space.Hence any product of operators Ô ∆ or ∆ Ô can be expressed as a differential operator acting on ∆.Therefore, as can be seen from Eq. ( 1), the same operator acting on ρ can be converted into a differential operator acting on the Wigner function.However, the infinitesimal volume elements dΩ n are not constant due to the curved phase space.It is therefore instructive to express the correspondence rules in terms of the contravariant coordinates (x 1 , x 2 ) = (θ , φ) of the phase space of a single spin and the metric tensor g µν which is given by The derivatives acting on the Wigner function are then given by covariant derivatives with x n = θ n , φ n .This yields correspondence rules such as [25] σz ρ ↔ A full list of these rules, but in the aforementioned different angle convention, is given in Ref. [25].
For general spin-j systems, exact expressions for the correspondence rules are known [35], but are complicated.They do have a simple semiclassical form in the limit j → ∞, but for j = 1/2, which is by far the most commonly considered case in many branches of physics, this semiclassical limit is not directly applicable.

An example for an exact FPE: spontaneous emission of a single two-level atom
Let us start by applying the correspondence rules such as Eq. ( 13) to the important simple example of spontaneous decay, where an exact FPE can be derived.Two-level atoms in free space can undergo spontaneous relaxation to the energetically lower state by emission of light quanta.This is due to the fundamental coupling of the atoms to the quantized electromagnetic field.The description of this phenomenon can be drastically simplified by assuming that the field is in equilibrium (which is the vacuum at optical frequencies) and by subsequently integrating out the field's degrees of freedom.A Born-Markov approximation then yields the effective Lindblad master equation [38] for the density operator ρ(t) of an individual atom.The rate Γ 0 is the Einstein A coefficient.Following the arguments of [25] we translate the master equation for ρ into an Fokker-Planck equation for the Wigner function W ρ (Ω, t): It has an equivalent set of Itô SDEs [39,40] dθ where dW φ is the differential of a Wiener process.These equations are exact and solving them is numerically inexpensive without further approximation.

The general case: Truncated Wigner Approximations (TWA) as diffusion approximations
Knowing the exact phase space formulation of the master equation shifts the quantum manybody problem of solving a large matrix differential equation of the density operator ρ(t), Eq. ( 9), to solving a high-dimensional partial differential equation with possibly infinitely many orders of derivatives for the c-number quasi-distribution W ρ (t).Except for special cases, such as the one discussed in Sec.2.3, both formulations are useless without the introduction of further approximations.
From the perspective of complexity, the core idea behind different variants of the TWA consists of neglecting higher order terms of the equation of motion of W ρ (Ω, t) such that the remaining expression is a covariant Fokker-Planck equation (FPE) in terms of suitable phase space variables Ω where ×2N is a positive semidefinite matrix.It then can be equivalently expressed by the set of SDEs [41] where d W ∈ 2N is a multivariate differential Wiener process.This and all further SDEs will implicitly be stated in the Itô calculus.
In a numerical implementation, we can efficiently compute N traj independent solutions of the SDEs [42], which we call trajectories.All relevant expectation values can then be directly calculated in the Wigner phase by using the relation The bars indicate the stochastic average where Ω (n) refers to the phase space coordinate of the n'th trajectory and where the approximation due to a stochastic error vanishes as we let N traj → ∞.In summary, this means that the time evolution of the Wigner function in TWA is governed by a diffusion process in the spherical Wigner phase space.From a physical standpoint, this truncation should be formulated in a systematic fashion which elucidates its validity in terms of a small parameter.
When adding spin-spin interactions, such as Ising-type couplings, the resulting equation for W ρ (Ω, t) is no longer of Fokker-Planck type and approximations are needed.A commonly used approximation is the discrete truncated Wigner approximation (DTWA) [24], which essentially amounts to a mean-field factorization of the Wigner function.This approach always produces deterministic equations and cannot account for the noise expected in dissipative systems.
As shown in [26,43] independent dephasing of spins can be incorporated in the DTWA, but the description of decay requires some ad-hoc modelling [44], which is not justified in general.Therefore it is not surprising that the standard DTWA cannot be applied to collective decay processes such as superradiance.We recently developed an alternative approach, termed hybrid continuous-discrete truncated Wigner approximation (CDTWA), which describes the time SciPost Physics Submission evolution of the many-body density operator by a continuous representation of the many-body spin Wigner function but samples the initial distribution from a discrete representation [25].The CDTWA incorporates (uncorrelated) decay and dephasing of the spins in a consistent way, but in the context of collective decay does not generally reveal which correlated terms can be neglected or not (see Sec. IV.E of Ref. [25]).

Truncated Wigner Approximation for Large Spin Ensembles with Collective Couplings
Neither the standard DTWA nor the CDTWA mentioned in the previous section are suitable for describing problems of collective couplings among spins.We now present an alternative approach based on an approximate form of the correspondence rules for collective spin processes and derive conditions for their validity.

Semiclassical limit of the correspondence rules for collective operators
If an ensemble of two-level atoms is confined to a small volume comparable in size with the wavelength of the dipole transition between the two states, the coupling to the quantized electromagnetic field leads to a correlated emission of photons known as sub-and superradiance.
Collective processes in an ensemble of N spins can be described in terms of collective operators where the "degree of cooperativity" is encoded in the weights {c n } = (c 1 , c 2 , . . . ) ∈ N .For describes the maximally cooperative case of an all-toall coupling, relevant e.g. for modelling Dicke superradiance, see Sect. 4, while a distribution of the c n 's peaked for some index n = j corresponds to the low-cooperativity case of short-range interaction.The action of Ŝ({c n }) on the state ρ can be exactly expressed as a differential operator acting on the Wigner function W ρ (Ω) in the phase space, however this differential operator does not have a simple form [34,35].For the resulting equation of motion for the Wigner function to be of practical use, we propose instead truncated correspondence rules with where s n (Ω) is given by Eq. ( 7) and L n is the angular momentum differential operator expressed in terms of covariant derivatives.Similarly we find the action for operators acting from the right-hand side.Note that the above contributions are the first and third term on the right-hand side of Eq. ( 13), whereas the remaining terms denote "quantum corrections" to this lowest order contribution.The intuition behind this truncation is simple: For a single spin-j, the same semiclassical limit can be obtained by letting j → ∞.This reveals that classical and quantum contributions separate in the Wigner phase space.
We note that this approximation leads to a Fokker-Planck equation for W ρ without higherorder derivatives if the master equation is at most bilinear in the collective operators.This allows for an efficient simulation in terms of SDEs.Eqs.(22) are the central element of our approach and form the basis of the simulations of collective decay phenomena discussed in Sects.4 and 6.

Validity of the approximate correspondence rules
We now discuss the range of validity of the truncated correspondence rules, Eqs.(22).To this end we first note that the density operator ρ of a system of N spins has the general form and µ n = (0, x, y, z), with σ0 = ˆ and s 0 = 1, from which we can immediately deduce Note that this expression, while being exact, is only of formal use as the sum contains an exponentially large number of terms.It does allow us, however, to explicitly calculate the exact Weyl symbol of operators such as Ŝz ({c n }) ρ through direct evaluation of Eq. ( 2) via Eq.( 25): Applying the spin algebra of the Pauli matrices and evaluating the individual Weyl symbols yields where ϵ i jk is the Levi-Civita symbol.The truncation approximation in Eqs.(22) of the same Weyl symbol is obtained, on the other hand, by applying the z-component of Eq. (22a) to the Wigner function in Eq. ( 25), which yields: To determine the error of Eq. ( 27) made by the truncated correspondence rule we define its difference to Eq. ( 26) In a similar way we can proceed with the xand y-components Ŝ x ({c n }) ρ and Ŝ y ({c n }) ρ.This gives the full difference vector where the superscript µ n indicates the µ n 'th row of the given matrix.Finally, the error can be quantified by the norm of the vector δ(Ω) We now evaluate the diagonal and non-diagonal parts separately.We find for the diagonal contribution: which follows from dΩ s as each contains factors we see that When ρ is the completely mixed state, we have µ n = 0 for every n and therefore the truncated correspondence rules are exact.For general states we can infer the error to scale as One recognizes that if the coefficients c n all have comparable magnitudes we have A necessary condition for the asymptotic correspondence rules to be valid is that the relative deviation to the exact Weyl symbol is small, i.e.

||δ(Ω)|| ||W Ŝ({c
Note that This expression can maximally scale as O(N ), in which case the truncation approximation Eq. ( 36) is satisfied for large ensemble sizes N .We now argue that this is the case if the dynamics of the system takes place in the subspace of states with large cooperativity.To this end consider the totally symmetric operators with c n ≡ 1.The total angular momentum operator Ŝ2 = Ŝ † Ŝ has eigenstates | j, m, α〉 with so-called cooperativity 0 ≤ j ≤ N 2 , projection |m| ≤ j on the z-axis and the parameter α distinguishing degenerate states.If ρ = | j, m, α〉〈 j, m, α| and j ̸ = 0, then Eq. ( 36) yields i.e. the cooperativity of the spin ensemble determines the validity of the asymptotic form of the correspondence rules.

Two-body interactions and collective dephasing
Before turning to specific applications of our TWA approach, let us discuss two special cases of collective spin-spin interactions and collective dissipative processes which are relevant e.g. for ensembles of two-level atoms coupled via a cavity field.
To describe the time evolution under the action of a collective interaction we can use the truncated correspondence rules of Eq. ( 22) resulting in The equivalent stochastic differential equations in θ n , φ n are in fact deterministic and quantum fluctuations enter only through the averaging over the Wigner distribution of the initial state.
A change of variables θ n , φ n → s n to Cartesian coordinates then gives equations of the type where S µ ({c n }) = m c m s µ m e µ with µ = x, y, z.This is a Larmor precession of the vectors s n about the cumulative magnetic field 2c n S µ ({c n }) and is equivalently predicted by a mean-field approximation and the standard DTWA.
In addition to unitary interactions described by a von Neumann equation, collective dissipative processes described by Linblad master equations are oftentimes of interest as well.A particularly simple case is that of collective dephasing for which we find an exact mapping to a FPE [34] γ This equation has the equivalent set of very simple SDEs It is not surprising that all angles φ n couple to the same noise dW , as their time evolution can equivalently be generated by a dynamic Hamiltonian contribution Ĥ = γ Ŝz ({c n })η(t) where η(t) is a white noise process with identical properties as dW .

TWA Description of Collective Light Emission
Let us consider N two-level atoms with arbitrary but non-overlapping positions r n coupled to the quantized electromagnetic field at distances comparable to the wavelength λ e of the two-level transition.In contrast to the case of atoms spaced at distances much larger than λ e , which allows a formal elimination of the coupling to the radiation field for each atom individually, leading to the effective Lindblad master equation ( 14), here radiative couplings between the atoms need to be taken into account.In addition we allow for a driving of the atoms by an external coherent light field which is polarized along the unit vector e c and has the wave vector k c = e n ω c /c with e n •e c = 0.The corresponding Rabi frequency for the j'th atom is Ω j = p • e c E(r j ), where p is the atomic transition dipole moment, which is assumed to be identical for all atoms.We denote the detuning between the classical field and the atoms as ∆ = ω c − ω e .Formally integrating out the electromagnetic field and using a Born-Markov approximation results in a master equation of the N atom system which reads [45] The effective Hamiltonian describes the coupling to the external coherent drive as well as the radiative coupling between the two-level atoms with rates J mn .These rates as well as the positive definite decay matrix Γ = GG T ∈ N ×N , where the choice of G is unique up to a unitary rotation, are given by the free space Green's tensor G E (r m , r n , ω e ) of the electric field Their explicit expressions are where r mn = r m −r n and e p (e r mn ) is the unit vector along the polarization p (the position r mn ) and k e = ω e /c.The diagonal elements Γ nn = Γ 0 are given by the Einstein A coefficient and we set J nn = 0, thereby absorbing it into the atomic detuning ∆.This contribution corresponds to the Lamb shift, which is however not correctly described within the dipole approximation of the atom-light coupling.In fact J nn diverges since r mn → 0 for m = n.In the following sections we will assume resonant driving of the atoms and therefore set ∆ = 0.An exact mapping of the master equation to phase space would go beyond a Fokker-Planck description, however the asymptotic correspondence rules of Eqs. ( 22) reduce it to where φ mn = φ m − φ n .Note that applying the approximate correspondence rules to terms such as Ŝ({c n }) − ρ Ŝ({c n }) − from either left-to-right or right-to-left produces a small imaginary contribution even though the master equation is real-valued.This produces a neglectable error but can be alleviated by taking the symmetric average of both variations which corresponds to just taking the real part of either one.The above FPE possesses an equivalent set of SDEs given by sin θ m J mn sin φ mn + Γ mn 2 cos φ mn d t sin θ m −J mn cos φ mn + Γ mn 2 sin φ mn d t with 2N independent Wiener increments.
Note that Eq. ( 48) does not reduce to Eq. ( 15) when taking the limit N = 1.The latter is exact due to a representation in terms of a suitable choice of operators in the single-spin Hilbert space whereas the former uses an expansion in the total angular momentum of the ensemble of spins and therefore becomes invalid at small N .
The single-particle terms can be treated exactly and yield which gives the following additional deterministic contributions to the above SDEs.
In the following sections we will investigate specific examples of atomic matter coupled to quantized light fields and demonstrate the strengths and weaknesses of the TWA by comparing its predictions of several observables to numerically exact results.An experimentally available observable is for example the total photon emission rate into all spatial directions.It is typically easier to detect the intensity of the emitted light into a solid angle with a direction defined by the unit vector e k or small areas obtained from an integration over some geometric configuration thereof.The photon emission rate along the unit vector e k is given by [45] where I 0 e k is the enveloping emission profile of a single atom.Note that this can be rewritten in terms of collective operators such that the Weyl symbol required for the calculation of the expectation value follows from applying the truncated correspondence rules from right to left.The aforementioned total radiated intensity is thus explicitly given as Unlike cumulant expansions, which only allow the investigation of expectation values up to the truncation order, the TWA can in principle be used to calculate arbitrary expectation values.
A prime example for this is the second order correlation function the detected light g (2)   e k which similarly evaluates to Just like in Eq. ( 48) the truncation produces a small imaginary component for the Weyl symbol if complex c n are considered.This is again alleviated by taking the symmetric average over the left-to-right and right-to-left application of the correspondence rules which effectively produces just the real part of the above symbol.Moreover we consider the spin squeezing parameter ξ 2 defined as where Ŝ = Ŝ x , Ŝ y , Ŝz T is the collective spin operator and ance of the operator Ŝe n = e n • Ŝ projected onto an axis that is orthogonal to the mean spin, i.e. 〈 Ŝ〉 • e n = 0.This minimal variance is not only of interest in quantum metrology, but furthermore a squeezing of ξ 2 < 1 implies entanglement [46].

Dicke Decay
To benchmark our method and to illustrate its strengths, we will first study the case where all atoms couple with identical rates Γ mn = Γ 0 , and where the unitary couplings are ignored J mn ≡ 0. This model was proposed by Dicke as an approximation to the radiative coupling of a free gas at very strong confinement [1].The model also typically arises in cavity-and waveguide QED.We fix the non-unique choice of G in Γ = GG T to G mn = Γ 0 δ n,1 .Substituting this into Eqs.(49) reveals that the phase space angles only couple to 2 of the possible 2N white noise processes.This is intuitive, e.g. from cavity QED, where these two degrees of freedom represent the noisy coherent amplitude dα = dW x + idW y of the bosonic cavity mode that adiabatically follows the state of the atoms.If the system is initially in the inverted state |e 1 e 2 . . .e N 〉 = | j = N /2, m = N /2〉, it can only decay along the states of maximal cooperativity j = N /2.Hence we expect the TWA be a good approximation at large N .The ensemble descends this ladder of states with initially increasing and then decreasing rates.This gives rise to the effect of superradiance, i.e. the emission of light at a rate faster than that of a single atom [1].Furthermore the restriction to just N + 1 states means that an exact and efficient numerical integration of the master equation in terms of rates is possible [3].
In Fig. 1 we compare the TWA prediction of the number of excitations and the total emission rate to exact results.The TWA results were produced using an Euler-Maruyama integration scheme [42] with a timestep ln(N )Γ 0 ∆t/N = 10 −3 and an averaging over 64 • 10 3 trajectories.They accurately reproduce the exact results.Even at small ensemble sizes of Since only the states of maximal j = N /2 couple to the vacuum state |g 1 g 2 . . .g N 〉, other initial states cannot fully emit their excitations.This gives rise to the effect of excitation trapping.For simplicity consider even N .If we assume the initial state to be the completely mixed state, given by the factorized Wigner function the steady state population can be determined by summing over the (2 j + 1)d j states in each j-ladder with degeneracy d j and with probability 2 −N each and multiplying by the population − j of the bottom state, leading to In Fig. 1 d) and e) we again see a very good agreement of the TWA with the exact results that improves as N increases.Furthermore, the Dicke decay is a prime example for revealing how superradiance emerges within a semiclassical framework.With the assumption that J mn = 0 and Γ mn = Γ 0 we can see that the SDEs of Eqs.(49) are closely related to the Kuramoto model [47] which describes harmonic oscillators with frequencies ω n and pairwise coupling rates K mn .If we compare this to the equations of the relative phases φ n of the two-level states, we can identify ω n = −∆ = 0 and K mn = 1 2 Γ 0 sin θ m cot θ n .The coupling is long-ranged and, due to the appearing θ m terms, time-dependent.Additionally, the phases are subjected to non-diagonal and non-linear noise.Nevertheless the origin of superradiant bursts in the Dicke model is related to the phase transition in the Kuramoto model from a completely incoherent state where all {φ n } are uniformly distributed to that of spontaneous synchronization.This emergence of synchronization φ m = φ n causes a dynamic shift of the changes dθ n and therefore of the total number of excitation 〈 Ŝz 〉.As a result, the photon emission rate will transition from individually radiating atoms γ(t = 0) ∼ N to collectively enhanced emission ∼ N 2 .In the Kuramoto model, synchronization is quantified by the order parameters where 0 ≤ r ≤ 1 is the coherence and ψ is the average phase.Individual trajectories of the Dicke decay in TWA, denoted by the subscript ( j), indeed share the feature of emerging transient coherence as is shown in Fig. 2.Even though the coherences r ( j) approach zero at short and long times, there is an intermediate window where they peak significantly.Around this peak, the change of the phases in time vanishes and the signal-to-noise ratio is strongly enhanced.At the same time, the slope of the number of excitations is minimal, i.e. a photon emission burst occurs.At first glance this coherent locking of phases might be surprising when compared to the rate equation of the density operator which does not show such an effect.We note however that the emerging average phase ψ ( j) of a single trajectory during the burst is uniformly distributed.By taking an additional trajectory average before computing the order parameters, the coherence vanishes at all times.

Dynamics of Spatially Extended Systems
Let us now turn to the more realistic spatially extended systems, where idealizations such as the Dicke decay are no longer sufficient.First, we consider the recently developed lightmatter interfaces based on regular arrays of atoms [48] with sub-wavelength lattice constants.In these arrays interference from the precisely positioned atomic emitters leads to pronounced collective responses despite a comparatively small number of atoms.Using atomic configura- tions inspired by these recent experimental advancements, we benchmark the performance of the TWA based on the truncated correspondence rules with numerically exact results.
The atomic ensembles are treated according to their full master equation of Eq. ( 44) including dissipation and the dipole-dipole interactions.The numerical predictions are obtained from solving the SDEs of Eqs.(49).For all examples we choose a timestep Γ 0 ∆t = 10 −3 and N Traj = 64•10 3 trajectories for the TWA simulations.We compare the semiclassical predictions to Monte Carlo wavefunction (MCWF) simulations obtained by using the QuantumOptics.jlpackage [49].These are, apart from stochastic fluctuations due to a finite number of trajectories, exact.All MCWF expectation values were calculated from N Traj = 10 3 trajectories.
Finally, we consider a dense elongated cloud of harmonically trapped atoms driven by a laser.We model the geometry and coherent drive after a recent experimental investigation [50] and study the coherence of the light emitted by the cloud.

Driven atomic arrays
Let us first consider an array of N = 16 atoms in a 4 × 4 quadratic lattice in the x-y-plane with lattice constant a = 0.8λ e .
Here and in the next section, each atom is assumed to have a dipole allowed transition from the ground state |g〉 to the excited state |e〉 with circular polarization e p = (1, i, 0, ) T / 2 along the z-direction.They are initially in the collective ground state |g 1 g 2 . . .g N 〉 and are driven by a plane wave which propagates perpendicularly to the array such that the Rabi frequencies are simply reduced to Ω n = Ω.
In Fig. 3 we see that the interplay of driving and dissipation leads to damped Rabi oscillations in the number of excitations and finally to a non-trivial steady state.At small driving Ω/Γ 0 < 1 the TWA does not reproduce the transient dynamics and steady state values obtained Figure 4: Dynamics of an initially inverted array at a spacing of a = 0.2λ e and without an external drive.a) Number of excitations as a function of time as predicted by our TWA method (blue), a second order cumulant expansion (purple, MF 2) and by an exact MCWF simulation (green).The shaded areas variances ∆S z /2.b) TWA (solid lines) and second order cumulant expansion (dashed lines) predictions of the number of excitations starting from the fully inverted state (green) and the collective ground state (blue).c) Distribution of the eigenvalues γ i of the dissipation matrix Γ and | j i | of the dipole-dipole interaction matrix J. Stacked points denote degeneracies.d) Weights of the most superradiant and subradiant eigenvectors of Γ .from MCWF simulations, which are still feasible for this small number of emitters.
As the driving increases, the match between semiclassical and exact dynamics improves.In the moderate to strong driving regime Ω/Γ 0 ≥ 1 we see a very good agreement across all observables.Most notably, the spin squeezing parameter is matched closely, suggesting an overall excellent prediction of general second moments.
The ever improving performance in the strong driving regime Ω/Γ 0 ≥ 1 can be explained by the competition of the driving and the dissipation.Here, only the cooperative, i.e. superradiant, modes significantly contribute to the dynamics.On the other hand, subradiant modes with their weak rates become insignificant for the overall dynamics and the steady state.

Inverted atomic arrays
Now we analyze the relaxation of an initially fully inverted state |e 1 e 2 . . .e N 〉 in the absence of a classical driving field, i.e. we set Ω 0 = 0.The atoms again form a 4 × 4 quadratic array in the x-y-plane with a smaller lattice constant of a = 0.2λ e .
In Fig. 4 a) we see the comparison between the TWA prediction, second order cumulant expansion (MF 2) and a MCWF simulation.The second order cumulant expansion [20][21][22] result was produced using the QuantumCumulants.jlpackage.The dynamics can be split into two time regimes: The superradiant regime Γ 0 t ≲ 1 and the subradiant regime Γ 0 t ≳ 1.In contrast to the second order cumulant expansion, the TWA closely matches the exact results during the superradiant burst.As the system transitions into subradiance, the semiclassical prediction converges to a finite number of excitations n 〈 σ+ n σ− n 〉 ≈ N /10, i.e. it gets stuck at a subradiant plateau which is unstable in a full quantum treatment.We verify that this is a steady state according to the TWA by comparing it to another semiclassical evolution starting from the collective ground state as shown in b) where, even in the absence of any excitation processes, the system evolves out of the collective ground state.The second order cumulant expansion does not suffer from this problem, however fails to match the early time dynamics and similar unphysical creation of excitations in systems with subradiant modes, though not as pronounced as the TWA, has also been observed in other systems [20,21].
This can be explained by the distribution of the eigenvalues of Γ and J as can be seen in c).At the timescale Γ 0 t ≤ 1 the rates γ i of the cooperative superradiant eigenmodes of the matrix Γ = U • diag({γ j }) • U T dominate the dynamics, whereas the eigenvalues j i of the dipole-dipole interaction matrix J only significantly contribute at Γ 0 t ≳ 1.
In d) we show the amplitude distribution of the most superradiant (green, j = N ) and subradiant eigenvectors (blue, j = 1) of U i j .The superradiant mode is strongly cooperative as all coefficients have the same sign and are of approximately equal magnitude.In contrast, the subradiant mode has alternating signs and couples the atoms with more varying magnitudes.The validity criterion of Eq. ( 36) is satisfied when the superradiant mode acts on the collective ground state, but for the subradiant mode it is not.The presence of several modes with eigenvalues 10 −2 < γ i /Γ 0 < 10 0 shows that such low-cooperativity effects already become relevant at the timescale of the simulation and therefore the TWA fails to escape the plateau.

Superradiance from an extended atomic cloud driven by an external laser
Lastly, we consider the radiative properties of an ensemble of N = 1400 atoms in an extended, three-dimensional harmonic trap, see Fig. 5 a).An experimental realization of this configura-tion was recently investigated [50].The authors claim that the cloud of atoms behaves like an effective driven Dicke model (DDM) along the elongated x-direction, i.e. it can be described by the master equation However, the disordered atomic positions lead to a reduced cooperativity which is captured by an effectively smaller ensemble size Ñ = µN .They found that µ ≈ 0.005 yields a good match between the experimental observations and the DDM.
To describe the experiment within the DTWA, the positions of the atoms are normally distributed with standard deviation ξ = (10, 0.25, 0.25)λ e along each dimension and vanishing mean, such that the cloud is cigar-shaped and has the reported 1/e 2 -radial widths.
The atoms are assumed to have circular polarized transitions in the y-z-plane, i.e. we choose e p = (0, 1, i) T / 2. They are initially in the collective ground state |g 1 g 2 . . .g N 〉 and are are driven by a plane wave with k c = 2π λ e e z such that we obtain Rabi frequencies Ω n = Ω 0 e ik c •r n .We perform simulations at varying driving strengths Ω 0 /Γ 0 = 1, 2, . . ., 10.
In Fig. 5 b) and c) the dynamics of the cloud at Ω 0 = 10Γ 0 , i.e. in the superradiant regime, is compared to that of a single atom in free space.The incoming field drives damped Rabi oscillations which are further suppressed due to collective effects.The photon emission rate I e y (t) along the short side also shows collectively-damped oscillations, but converges to the singleparticle emission rate.In stark contrast to this, the emission into the x-direction shows a strong initial superradiant burst.The steady state photon emission rate per atom I e x (Γ 0 t = 10)/N is also greatly enhanced.The second order correlation function of the light in the y-direction immediately saturates to that of a single thermal mode, demonstrating that the atoms emit independently into this direction and that no coherent locking occurs.In contrast, photons emitted into the elongated direction show a much stronger degree of coherence.
In d) and e) we investigate the steady state radiation according to the DTWA.Here, at β = 2Ω 0 /Γ 0 Ñ ≈ 1.4 the transition from the magnetic to the superradiant regime of the DDM occurs.We compare the coherence of the emitted light to the experimental results and the corresponding DDM prediction [50] and see a much closer agreement with the DDM.The higher value of g (2) in the experimental data suggests additional dephasing effects, e.g.due to the thermal motion of the atoms, which is not included in the DDM and DTWA descriptions.

Conclusion
We here presented a semiclassical, numerically efficient approach to describe the many-body dynamics of spins with collective interactions and dissipation.The approach is an extension of the discrete truncated Wigner approximation [24], which approximates the time evolution of a physical state in the Wigner phase space by a diffusion-like process taking into account classical and leading-order quantum fluctuations.The equation of motion of the Wigner distribution of the many-body density matrix can be cast into a differential form by applying correspondence rules, i.e. the action of an operator on the state translated into phase space.We proposed a specific truncation of said correspondence rules by only keeping lowest-order contributions which maps the Lindblad master equation of interacting two-level systems to a Fokker-Planck equation with positive diffusion.The latter allows for a numerically inexpensive propagation in time by solving only linearly many stochastic differential equations.
We derived quantifiable conditions for the validity of the approximation in terms of an upper bound on the error that is introduced by the truncation.We showed in particular that the truncation becomes exact if the many-body dynamics is dominated by degrees of freedom with high cooperativity in a large ensemble.Thus the method is ideally suited for the analysis of collective processes such as superradiant emission of light in atomic ensembles.We benchmarked our method against exact results for the Dicke decay, which can be obtained without further approximations and found excellent agreement that improves with the number of atoms in the ensemble.In the case of small atomic arrays, we compared predictions from our semiclassical approach with exact Monte Carlo wavefunction results and showed that early superradiant timescales are well captured, however longer subradiant timescales cannot be reliably described.When the array is coherently driven with Rabi frequencies at or above the single-particle linewidth, the influence of the subradiant modes becomes negligible and the emerging dynamics is again well captured within the semiclassical approximation, while a second order cumulant expansion shows major deviations from exact results.Furthermore we study the dynamics of a driven, harmonically trapped, spatially extended ensemble of quantum emitters and calculated its population dynamics, direction resolved photon emission rates and their corresponding degree of coherence expressed in terms of g (2) correlations.The experimental detected light is more thermal than the TWA simulation and the prediction from a simpler theoretical model, which suggests the existence of dephasing mechanisms not yet present in the current theory.
Our approach paves the way for studying strongly cooperative effects in large and spatially extended ensembles of two-level systems.Specifically, recent light-matter interfaces such as trapped gases and atomic arrays and their non-linear response to incoming coherent light [48,50] can be studied with ensemble sizes much larger than N ≃ 50 as is considered in this work.Typically, not much analytical progress can be made in such systems and methods based on tensor networks do not work reliably due to the high dimensionality of these setups and intrinsic long-range interactions.
Future works will investigate whether the truncation, which so far is a diffusion approximation, can be improved by extending it to a jump-diffusion approximation by including classical Poissonian jump processes.This might allow for an extension of the validity of our theory to regimes of moderate or even low cooperativity.Motivated by the excellent agreement of spin squeezing and coherence of the emitted light, we believe that our method can be used to analyze the generation of non-classical states of realistic, experimentally realizable atomic configurations and their emitted light.

Figure 1 :
Figure 1: Dynamics of the Dicke decay for varying ensemble sizes N .For an initially inverted system the a) number of excitations and b) emission rate of are shown as predicted by the TWA (solid lines) and exact results (dashed lines).c) Absolute difference of the exact excitation number and TWA prediction.For the initially fully mixed state d) depicts population trapping of the excitation number.The dashed lines indicate the exact steady state populations.e) Time evolution of the relative deviation of the TWA population prediction to the exact steady state.

Figure 2 :
Figure 2: Time evolution of three sample trajectories according to the Dicke decay in TWA with N = 256 atoms.Depicted are a) the number of excitations, b) coherences and c) average phases of the atomic ensemble.The dashed vertical lines denote times of peak coherence.

Figure 3 :
Figure 3: Dynamics of a coherently driven 4 × 4 quadratic atomic array with lattice spacing a = 0.8λ e .a) Sketch of the array that is aligned in the x-y-plane and driven by a plane wave propagating along the z-axis.b) Total excitation number, c) total photon emission rate into free space and d) spin squeezing for a coherently driven atomic array as a function of time.The TWA predictions (solid lines) are compared to MCWF simulations (dashed lines).The different colors denote varying Rabi frequencies Ω.

Figure 5 :
Figure 5: Dynamical and steady state radiation of a driven, cigar-shaped cloud of N = 1400 atoms.a) Sketch of the geometry of the driven cloud.b) Photon emission rate into the x-and y-direction and c) corresponding coherence of the emitted light at Ω 0 = 10Γ 0 as a function of time.Steady state values of d) the photon emission rate in x-direction and e) coherences as measured (green circles), predicted by a driven Dicke model (blue solid line), both taken from [50], and the DTWA (purple diamonds) at varying Ω 0 .
r mn | 2 cos(k e r mn ) k e r mn − 1 − 3|e p • e r mn | 2 sin(k e r mn )(k e r mn ) 2 +