Phase space theory for open quantum systems with local and collective dissipative processes

In this article we investigate driven dissipative quantum dynamics of an ensemble of two-level systems given by a Markovian master equation with collective and non-collective dissipators. Exploiting the permutation symmetry in our model, we employ a phase space approach for the solution of this equation in terms of a diagonal representation with respect to certain generalized spin coherent states. Remarkably, this allows to interpolate between mean-field theory and finite system size in a formalism independent of Hilbert-space dimension. Moreover, in certain parameter regimes, the evolution equation for the corresponding quasiprobability distribution resembles a Fokker-Planck equation, which can be efficiently solved by stochastic calculus. Then, the dynamics can be seen as classical in the sense that no entanglement between the two-level systems is generated. Our results expose, utilize and promote techniques pioneered in the context of laser theory, which we now apply to problems of current theoretical and experimental interest.


I. INTRODUCTION
Quantum optical models, such as the Dicke model, where an ensemble of identical twolevel atoms is coupled to a cavity mode, have been among the first systems where emergent collective behavior such as super-radiance has been investigated [1]. The cooperative effects can be described in terms of collective atomic operators in such a way that the ensemble of N two-level atoms is treated as a large spin of length j = N 2 . Interestingly, such a cooperative state can still possess non-trivial quantum correlations [1]. The proposal on how to implement the Dicke model in an optical cavity QED system [2] and the experimental realization using a super-fluid gas trapped inside an optical cavity [3] are important milestones, which have inspired more studies on the role of collective effects in quantum phase transitions, such as [4,5].
Purely collective dynamics cannot always be achieved due to experimental conditions, for example when the ensemble of two level systems is inhomogeneous [6] or when the individual two level systems couple to different reservoirs [7]. In these cases, the dimensionality of the relevant Hilbert space in general grows exponentially with the system size. An interesting alternative regime exists when in addition to collective processes only local processes occur which preserve permutation invariance. This holds if all non-collective processes are identical for each constituent in the ensemble. In this regime, permutation invariance can be utilized to reduce the effective dimension of the problem [8]. Formally, this can be understood as follows. In an ensemble of N two-level systems, collective processes confine the dynamics on a subspace H j ⊂ C 2N spanned by the Dicke states |j, m , where −j ≤ m ≤ j, whereas the permutation invariant local processes couple the different subspaces H j spanned by Dicke ladders labeled by differing j . The reduction in the effective dimension emerges due to high degeneracy of the Dicke states [9].
When local and collective processes coexist, a host of interesting phenomena can be studied. For example, incoherent pumping on a super-radiant ensemble of two level systems can lead to robust steady state super-radiance [10]. Competing collective phenomena, such as driving and dissipation can lead to non-equilibrium phase transitions [11,12] and recently the robustness of such transitions against local dephasing has been investigated [13].
In this article we will investigate a paradigmatic driven, dissipative open system consisting of N two level atoms which is affected by both collective and permutation invariant local dissipative processes. We describe the dynamics of such a system with the following GKSL type master equation [14,15] where collective driving is given by with B ∈ R 3 , σ i λ = σ i for i = z, ± are Pauli spin matrices σ z , σ ± = σ x ± iσ y acting on system λ, and J is the collective spin operator with components J z and J ± = J x ± iJ y . The collective dissipation is given by whereas the permutation invariant local dissipation is given by We assume that all of the local and collective rates, γ k and κ k respectively, are non-negative, so that the dynamics generated by the master equation is completely positive. Note that the non-collective term does not commute with the total angular momentum J 2 . It therefore couples different eigenspaces of J 2 in contrast to the collective terms.
Instead of using the basis of Dicke states directly, we will employ a phase space approach by representing the state of the open system in terms of generalized spin coherent states and the associated P -function [16,17]. This type of approach to collective phenomena has been already discussed in the quantum optics community during the 70s and 80s of the previous century, in the context of cooperative fluorescence [18][19][20][21]. The generalized coherent state approach that we take has been pioneered in [22][23][24] and used also more recently in [25].
Here we show that in certain parameter regions this method can be used to map the solution of master equation (1) to a Fokker-Planck equation for any system size. Then the dynamics is classical in the sense that no quantum correlations between the two-level systems build up.
In the method the system size is merely a parameter determining the strength of diffusion in phase space. Therefore, it allows to analytically interpolate between finite-size and meanfield theory.
The outline of the article is the following. First we introduce a generalization of spin coherent states in section II. In section III we derive the equation of motion for the associated P -function. The phase space approach provides immediately a consistent mean field theory, which we explore in section IV. In section V we focus on the exact semi-classical regime of our model, where the equations of motion can be efficiently solved. At last, we conclude with discussion and outlook in section VI.

II. GENERALIZED SPIN COHERENT STATES
The main technical tool of this paper is to expand the state of all identical two level systems in terms of a generalized P -representation. More precicely, consider the density operator α( r) of the simplest product state where all two level systems are identical This state is positive iff the Bloch vector r has a length smaller or equal to 1. It is pure iff | r| 2 = 1. We show later on that in this case the definition (5) reduces to the standard spin coherent states in the symmetric Dicke subspace with j = N 2 . Due to this property, we call these states the generalized coherent states. However, we stress that these are not pure states and do not have the same group theoretical interpretation as regular coherent states [26]. We will later express the state in terms of a diagonal representation using a P -function associated with these generalized coherent states. Firstly, let us discuss how certain operators act on these states. A direct calculation shows that where the derivatives are with respect to r, i.e. ∂ i = ∂/∂r i , and summation over repeated indices is implied. With these relations we can replace the action of operators acting on generalized coherent states by differential operators. One example is the total angular momentum operator J This gives us the algebra of spin coherent states. It is analogous to well known relations such as a |z = z |z in the case of bosons. With these rules we can express the action of the Hamiltonian part, as well as the collective dissipative part of the Lindbladian (1), onto a coherent state as a differential operator with respect to the coherent state label r. Due to our more general definition of spin coherent states, this is also true for the local dissipators.
Consider for instance identical local dephasing (in the σ z basis with rate γ z /2) of all two level systems Now the action of any superoperator on a two level system can always be written as a first order differential operator with respect to the Bloch vector. This is obvious because the state itself is linear in r. In the example of dephasing we find with canonical notation r 1 = x, r 2 = y, r 3 = z. Inserting this in the above expression and using the product rule one finds for the coherent state Of course, the same calculations can be done for local decay and pump channels. In summary, we can now express the action of the entire Lindblad superoperator on a generalized spin coherent state as a differential operator with respect to the state label, that is the vector r.
Let us briefly make the connection to the standard spin coherent states. These are usually defined in an eigenspace of J 2 with eigenvalue j = N 2 (symmetric Dicke states). They are rotations of the 'fully polarized' state α( e z ). We recover these states from our generalized coherent states by setting | r| = r = 1. For better illustration, let us parameterize the vector with spherical coordinates where one may identify η = cos θ. The differential representation of the dissipators in these coordinates are provided in appendix B. The main insight is that for r = 1 in the collective operators, no derivatives with respect to r appear. The radius is preserved because the dynamics does not lead out of the symmetric j = N/2 subspace. This is different in the case of local dissipation. The local dissipators are not confined to this subspace which is reflected in the occurence of derivatives with respect to r. If the dynamics remains in the Dicke subspace however, the equations of motions which we derive in the next section will automatically conserve r = 1 so that both cases are covered in the formalism. Now as an ansatz we express the state with a diagonal representation in terms of the generalized coherent states with a generalized P -function P ( r), which is a quasiprobability distribution on the unit ball in R 3 centered at the origin. As α( r) is a polynomial in r of order N , it is obvious that P is not unique. If ρ is a normalized state, the distribution is normalized as well If P ( r) is a positive function, the state is also positive by construction. However, P ( r) does not need to be positive. In fact, if P ( r) is positive for all r, then the state is separable by definition, as it is a classical mixture of separable states. An entangled state necessarily has a negative P function. Once the P function of the state is known, all observables can be computed easily For example, in the case of the angular momentum operator one finds tr α( r) J = N/2 r so that

III. EQUATION OF MOTION FOR THE P-FUNCTION
With the differential form of operators at hand, it is a straightforward task to derive an equation of motion for the P -function of a state evolving under the GKSL equation (1). One simply plugs in the expression (13) in the equation of motion Using the differential representation, we can express the action of the GKSL-generator as a second order differential operator To find the "forward in time" equation of motion for P one performs partial integration under the integral, neglecting boundary terms. Comparison of the integrands on the left and right hand sides gives the following partial differential equation if the dynamics does not generate entanglement between the two level systems. Indeed, we do find certain parameter regimes in our model where this holds true. Then the partial differential equation is well behaved and can easily be integrated numerically using stochastic differential equations. We give a detailed description of this in section V.

IV. MEAN FIELD THEORY
From equation (19) follows directly a 'mean field' approximation by neglecting the diffusion term. In general, this approximation is exact in the thermodynamical limit N → ∞ because then the diffusion vanishes exactly due to the prefactor 1 N . The remaining equation of motion is first order and can be solved by the method of characteristics The solution is given by the time evolution of the initial distribution along the mean-field trajectories P t ( r) = d 3 r 0 P 0 ( r 0 )δ( r − r(t)) , r(0) = r 0 ,˙ r(t) = a( r(t)) .
For our master equation the mean field equations of motion are explicitlẏ where we have neglected subleading terms of order 1/N consistent with neglection of the diffusion. Note that collective dephasing does not influence the mean field dynamics. Except for some fine tuned cases [12,27], mean field theory characterizes the steady state phases of the model [28]. Phase diagrams can be derived by finding the fixed points of (22)

V. EXACT SEMI-CLASSICAL REGIME
Going beyond mean-field one must include the diffusion term in the P -function evolution equation. This term is generated by collective dissipation in the master equation (1). Using   spherical coordinates, the diffusion matrix reads If the diffusion matrix is not positive semidefinite, the partial differential equation is no longer parabolic and numerical solution with, for example, finite element methods, is unstable.
There exist parameter regimes where this matrix has only positive eigenvalues. In this case the dynamics is described by a Fokker-Planck equation and can be solved with stochastic methods. An obvious example is the case of dephasing only, i.e. κ − = κ + = 0. Dephasing can be realized by a stochastic unitary evolution [29][30][31]. The corresponding stochastic Hamiltonian is non-interacting, so that the full dynamics can be mapped to a classical stochastic process. This is reflected by the positivity of the diffusion matrix. More generally, whenever the loss and the pump rates are identical κ + = κ − , the diffusion is positive.
Let us now focus on the purely collective case with no local decay processes. Then the Fokker-Planck equation preserves the radial direction r, so that we can set r = 1 and recover the canonical spin coherent state P -function. Then, the matrix (23) has two non-zero eigenvalues given by We see that Eq. (19) is a proper Fokker-Planck equation if both λ 1 , λ 2 are positive. λ 1 is always positive since −1 ≤ η ≤ 1 and λ 2 is positive for all η and φ if Satisfying this condition, the open quantum system evolution described by Eq. (1) can be mapped to a classical diffusion process on the sphere. Nevertheless the state lies in a finitedimensional Hilbert-space and contains quantum fluctuations due to the non-zero overlap of coherent states. Relation (26) can be understood in the sense that one simply has to add enough dephasing to compensate the positivity violation of the diffusion matrix due to collective losses. We point out that adding dephasing does however not influence the mean-field theory and thus the different phases of the model.
The standard way of solving the Fokker-Planck equation is to consider the set of stochastic differential equations where dξ i is a real Ito increment with properties dξ i dξ j = δ ij dt and the drift terms are given in appendix B. In the absence of the noise terms we obtain again the mean-field theory, as in the last section. The stochastic equations provide an efficient way of calculating the P -function of the system for any system size N when condition (26) is satisfied. Since in Eq. (27) N is merely a parameter determining the strength of the diffusion, the numerical effort for computing the state with this stochastic method is independent of the dimension of the underlying Hilbert-space. The method becomes more efficient than direct numerical integration of the master equation for moderate system sizes, i.e. when the mean field approximation is not yet applicable. We find that expectation values converge quickly with respect to the number of sampled trajectories, as we can see from the following example.
We consider the cooperative resonance fluorescence model [18][19][20][21] with additional collective dephasing. The model is described by the master equation (1) with κ − = B y = B z = 0, κ + = κ z and without non-collective terms γ + = γ − = γ z = 0. It features a mirror symmetry J x → −J x which is spontaneously broken when the driving exceeds the collective damping |B x | > κ + . If the damping is large, as in Fig. 2 (a), a localized steady state is reached as predicted by mean-field theory. In the symmetry broken regime the dynamics is characterized by damped oscillations towards a steady state with large variance, as seen in Fig. 2 (b).
This is an example where the mean field theory does not predict correctly the steady state and the finite system size remains relevant even for large N , see Ref. [12] for a more detailed discussion.
In Fig. (3) we show an example of a steady state P-function of the model as a density on the sphere for different system sizes, in the case of large damping (|B x | > κ + ). As the system size N increases, the distribution gets more localized due to weaker diffusion. For large system sizes the steady state distribution is sharply peaked at the stable fixed point of the corresponding mean-field theory.

VI. DISCUSSION
In this article we have further developed a particular phase space picture for driven dissipative spin systems first introduced in the context of laser theory [24]. Our approach incorporates both collective and permutation invariant local decay processes into a unified framework, which allows for a very simple derivation of evolution equations for the corresponding P -function.
As is well known, the coefficient matrix for the second order terms in the equation of motion for phase space quasiprobability distributions is typically not positive semidefinite. We have examined conditions under which initially positive phase space densities remain positive for all times. Then the dynamics can be considered classical in the sense that no entanglement is generated between the two-level systems. In particular, in the case of purely collective dynamics, presence of sufficiently large dephasing ensures the positivity. Another interesting case is the 'infinite temperature' limit κ + = κ − , which results in positive diffusion even in presence of non-collective processes. In case the diffusion matrix is not positive, solving the problem using phase space techniques is challenging as this leads to regions where P is negative and the evolution equation is no longer parabolic.
The phase space approach provides immediately a consistent mean field picture just by neglecting the subleading second order terms in the evolution equation. By analyzing the stability of the fixed points of the mean field equations, phase diagrams for non-equilibrium phases can be obtained. However, scenarios are known where finite-size corrections become relevant even for large system sizes [12]. In the case of positive diffusion, such finite-size corrections can be incorporated exactly by adding noise to the mean-field equation. This remarkable feature allows to efficiently solve the problem with stochastic methods, where the numerical effort is independent of the dimension of the underlying Hilbert space.
The formalism presented in this paper can be straightforwardly generalized for SU(n)-type models [32,33], by using generalized Bloch-representations, for example for three-level systems, and defining the coherent states accordingly.
We believe that the presentation of the phase space methods in this article is easily accessible and can be applied effortlessly to problems of current experimental and theoretical interest. In this appendix we provide the full expressions for the diffusion matrix and the drift term. These can be derived by expressing the action of the generator of master equation (1) onto a spin coherent state Expressing the state via (13), this results in an evolution equation for the P -function of model (1) The drift term a reads and the symmetric diffusion matrix D = D T is given as D xx = 2 κ − z −x 2 + z + 1 + κ + z x 2 + z − 1 + κ z y 2 , D xy = −2xy(κ z + z(κ − − κ + )) , D xz = x κ − x 2 + y 2 − (z + 1) 2 − κ + x 2 + y 2 − (z − 1) 2 , D yy = 2 κ z x 2 + κ − z −y 2 + z + 1 + κ + z y 2 + z − 1 , D yz = y κ − x 2 + y 2 − (z + 1) 2 − κ + x 2 + y 2 − (z − 1) 2 .