Atom-only theories for U(1) symmetric cavity-QED models

We consider a generalized Dicke model with U(1) symmetry, which can undergo a transition to a superradiant state that spontaneously breaks this symmetry. By exploiting the difference in timescale between atomic and cavity dynamics, one may eliminate the cavity dynamics, providing an atom-only theory. We show that the standard Redfield theory cannot describe the transition to the superradiant state, but including higher-order corrections does recover the transition. Our work reveals how the forms of effective theories must vary for models with continuous symmetry, and provides a template to develop effective theories of more complex models.

Introduction-While phase transitions in equilibrium many-body systems have been extensively studied and are broadly understood [1, 2], the critical behavior of driven-dissipative out-of-equilibrium systems [3][4][5][6] poses open questions. A central challenge in numerical exploration of such systems is the exponential growth of Hilbert space dimension with problem size. As such, any ability to reduce Hilbert space size, e.g. by a priori identifying a low-energy (slow) subspace can be very powerful [7][8][9][10]. A widely-used approach to derive a reduced model is Redfield theory [11,12]. Previous work [13] has shown this works well for the steady states and collective modes of the Dicke model. However, as we show below, Redfield theory can fail for models with U(1) symmetry, requiring higher-order approaches [14].
In this Letter, we consider a two-mode generalized Dicke model which has U(1) symmetry [38], and discuss how an effective theory can be developed. Given the mo- tivation above, there are certain conditions required of an effective theory: it must describe the transition to a symmetry broken state, and must correctly describe the frequencies and damping rates of the low-energy (soft modes) associated with this symmetry breaking. In the following, we first show why standard Redfield theory fails, and then present an alternative method which succeeds. We also discuss how, for such an effective model, we can derive semiclassical equations of motion, applicable in the large N limit.
U(1) symmetric model -We begin by introducing the U (1) symmetric model that we will consider, and summarizing its mean-field behavior. The model, introduced in Ref. [38], describes N two-level systems-described via a collective spin degree of freedom-interacting with two cavity modes. As shown in Fig. 1(a,b), this can be realized by a Raman driving scheme which couples two low-energy atomic states [17]. As shown, two transition pathways exist. Each pathway involves a different cavity mode-this is key to realizing a U (1) symmetry. Including cavity losses, we find an equation of motion for density matrix of the total system, ρ t : Here a, b are the two cavity mode annihilation operators, while S is a collective spin of the atoms, with modulus N/2. Cavity loss at rate κ is described by Lindblad terms L[X] = 2XρX † − {X † X, ρ}. In the Hamiltonian, the energies ω A,B describe the cost of scattering a photon from the pump into each cavity mode, while ω 0 is the splitting of the two hyperfine atomic levels. The effective coupling between light and matter is given by g = g 0 Ω p /∆ 2 a , where g 0 is the bare atom-cavity coupling, Ω p the Rabi frequency of the pump field, and ∆ a is the detuning between the pump and the atomic resonance, see Fig. 1(b). As shown in Ref. [38], this . This corresponds to a U (1) symmetry.
For large N , the composite atom-cavity system is well described by semiclassical equations [50] for a , b , S , which show two distinct types of steady state. At small g, there is a normal phase with S ± = a = b = 0, respecting the U (1) symmetry. At large g, there is a superradiant phase with a non-zero photonic field, a , b = 0. Thus, at a critical value g = g c -which, for ω A = ω B = ω obeys 2g 2 c N = ω 0 (ω 2 + κ 2 /4)/ω-the system undergoes a continuous phase transition to a state which spontaneously breaks the U(1) symmetry as shown in Fig. 1(c). Because of the cavity loss, these steady states are attractors of the dynamics, and the system undergoes damped relaxation towards these states.
In typical experiments, there is a separation of timescales between the atomic and cavity degrees of freedom, κ ω 0 , g √ N . We thus next consider how to adiabatically eliminate the cavity degrees of freedom and still describe the same behavior as discussed above.
Adiabatic elimination at semiclassical level -We first consider adiabatic elimination of a , b from the semiclassical equations given in [38]. This yields an equation of motion for S = S . The resulting equation is conservative, defined by a Poisson bracketṠ = {S, H sc } where, for ω A = ω B = ω, the classical Hamiltonian takes the form H sc = ω 0 S z + 2g 2 ω ω 2 +κ 2 /4 (S z ) 2 . While this Hamiltonian has a ground state phase transition at g = g c as expected, the purely conservative dynamics is in contrast with the dissipative evolution expected for an open system. Similar behavior was found for the single-mode Dicke model [13,51]. The results of [13] suggest that to recover the correct dissipative dynamics one should instead eliminate the cavity degrees of freedom in a quantum model and then derive the semiclassical limit.
Redfield Theory-To eliminate the cavity modes in a full quantum model, we use the standard Redfield ap-proach [11,12]. Specifically, we take the collective spin as the system, and all other modes form the bath. The system-bath coupling in the interaction picture is H I = g[S + (t)X(t) + S − (t)X † (t)], where S ± (t) = S ± e ±iω0t and X(t) = a(t) + b † (t). The time evolution of X(t) is discussed below. Redfield theory states: Evaluating this requires two-time correlations of the bath operators X(t).
These are found by solving Heisenberg-Langevin equations for the cavity modes, including loss, giving In the Schrödinger picture, the 2 nd -order Redfield equation (2RE) is thus: While this equation includes dissipative effects, as we discuss next, it does not show a phase transition with increasing g. To understand why, we first rewrite this equation in Linblad formρ The U(1) symmetry means the steady state must commute with S z , which implies ρ = Since this ratio is independent of g, no transition occurs at g = g c . At large N we find , the system is always in a normal or inverted state. This absence of a phase transition depends on two features of the equation. First, symmetry ensured that both the effective Hamiltonian and steady-state density matrix are diagonal in the S z basis, so they commute. This means the steady state depends only on the dissipative terms. Such a statement will always be true for U(1) symmetric models. Secondly, the dissipative terms in the Redfield equation are all proportional to g 2 , so no g dependence occurs in their ratio. This statement is not generic, so we next consider how contributions of higher-order in g change the equation.
4 th -order Keldysh-Redfield Theory-A systematic method to derive higher-order density matrix equations was introduced by Müller and Stace [14], making use of Keldysh diagrammatic perturbation theory. This technique allows one to take into account all contributions at each order while avoiding double counting. We write the density matrix equation in the formρ = D 0 ρ+D 2 ρ+D 4 ρ, where D 0 ρ + D 2 ρ is given in Eq. (4), and D 4 ρ is 4 th -order in H I . Crucially, the diagrammatic expansion ensures the terms in D 4 ρ are not separable-i.e. they correspond to genuine 4 th -order processes, not products of 2 nd -order terms.
FIG. 2. Example diagrams at 4 th order. The solid black lines are periods of free evolution of the system, interrupted by the action of the interaction Hamiltonian at times t1, t2, t3 (purple dots). These vertices are connected by purple dashed lines representing the cavity mode correlation functions. Figure 2 shows two examples of 4 th -order diagrams. As described in Ref. [14], the solid horizontal lines indicate the branches of the Keldysh contour. The system undergoes free evolution along these branches, interrupted by the action of H I at times t > t 1 > t 2 > t 3 (purple dots). Points on the lower branch correspond to operators to the left of the density matrix, while those on the upper branch are to the right. Because the bath is quadratic, expectations of bath operators factorize into pairwise correlations (purple dashed lines). The form of H I means each dashed line must connect opposite S ± operators. Following these rules, after integrating over t 1 , t 2 , t 3 the diagrams in Fig. 2 correspond to: Overall, at 4 th order, there are 32 diagrams. Considering the patterns of S ± operators, each diagram contributes 4 terms. These are written in full in the supplemental material [52].
The resulting equation is not of Lindblad form [53] and so does not necessarily preserve positivity-this is shown in [52], where we diagonalize the Lindblad-Kossakowski matrix. As discussed in many other contexts, such nonpositive equations can nonetheless predict correct behavior [13,[54][55][56][57][58]. In 2 nd -order theories, a Lindblad form equation is known to arise from secularizating the Redfield equation [59]-i.e., deleting terms which are time dependent in the interaction picture. Since our equation is U (1) symmetric, it already has a secularized form.
While U (1) symmetry still means the density matrix is diagonal in the S z basis, the presence of the 4 th order contribution gives a non-trivial dependence on g. In particular, as shown in Fig. 3(a), the steady state of this 4 thorder Keldysh-Redfield equation (4KRE) shows a transition to a superradiant state. As N increases, the results converge to the mean-field predictions of the full atomcavity model, with a discontinuity in S z at g = g c .
Semiclassical equations-We now consider whether we can use the 4KRE to derive semiclassical equations which capture the dissipative phase transition. A naive approach to this is to assume that terms S α S β S γ can be factorized. Applying such an approach directly to 4KRE fails, in the sense that after factorization, this yields equations for which S x,y = 0, S z = −N/2 is always a stable solution, in contrast to the clear instability seen in the full quantum solution. The structure of the full quantum solution (i.e., a density matrix diagonal in the S z basis) suggests the origin of this failure: while U (1) symmetry ensures S x,y is zero for the full solution, products such as (S x ) 2 can be non-zero. The diagonal structure of the steady-state density matrix suggests a better way to approach a semiclassical limit is by considering the Gaussian form of the probabilities P M (see [52]), and characterising P M by its first two moments. This yields a form of cumulant equation (CE), as used elsewhere [50,60]. Considering coupled equations for S z S z , S z (given in [52]), one finds results (dotted lines in Fig. 3(a)) that match 4KRE well for a range of N . At N → ∞, P M is sharply peaked so S z S z → S z 2 , giving a single equation for S z : Liouvillian spectrum-To explore whether the atomonly equations we have derived correctly capture the slow dynamics, we next consider the Liouvillian spectrum. Due to the U (1) symmetry, the quantum dynamics decouples into sectors labeled by index k, ρ = The matrices L (k) in each sector can be numerically diagonalized to obtain eigenvalues λ (k) i . The real part of these eigenvalues describe the relaxation rate toward the steady state. For our 4 th -order approach, the matrices L (k) take a simple structure-only terms with M = M, M ± 1, M ± 2 are non-zero, yielding a pentadiagonal matrix [52]. Together, the separation into sectors and this banded structure mean we can numerically find the eigenspectrum for relatively large N .
The k = 0 sector describes the dynamics of the populations (in the S z basis). As such, the eigenvalues in this sector provide information about the evolution of S z , and in particular, the damping rate towards steady state. The first non-zero eigenvalue [61] in the k = 0 sector, λ 1 , is shown in Fig. 3(b). We can also compare this eigenvalue to the decay rate found by linearizing the semiclassical Eq. (6): shown as the black solid line in Fig. 3(b). While the first non-zero eigenvalue for k = 0 tells us about relaxation to the steady state, it does not describe the slowest dynamics of this U (1) symmetric systemi.e., the smallest non-zero eigenvalue, also known as the Liouvillian gap. As explained in Refs. [49], when an open system spontaneously breaks a symmetry, the Liouvillian gap should vanish as N → ∞ throughout the symmetrybroken phase. This occurs because spontaneous symmetry breaking means more than one steady state is possible; since any mixture of two steady-state density matrices is also a steady state, an extra zero mode must arise. When the spontaneously broken symmetry is continuous, this also relates to the Goldstone mode.
For our model, the gapless mode is associated with how the U (1) symmetry is broken. As such, it must involve terms which are off-diagonal in the S z basisspecifically terms which correspond to the long-time evolution of S + (t)S − (0) . This means we should consider the smallest eigenvalue (by real part) of the |k| = 1 sector, λ (1) 0 . This eigenvalue is shown in Fig. 4(a), as a function of coupling g. One may see that the gap reduces with increasing N . However, a detailed analysis of the Liouvillian gap as a function of system size, Fig. 4(b), shows a non-zero gap remains at N → ∞. Specifically, considering g > g c , we see this eigenvalue matches well to λ (1) 0 = A + B/N , with a finite intercept A. As noted above, the 4KRE is not of Lindblad form (i.e. is not completelely positive), and it is possible this may be associated with the non-vanishing Liouvillian gap. However, the equation is Hermitian, trace-preserving, respects the U (1) symmetry, and predicts the correct steady states. Moreover, other examples of non-positive Redfield equations do show the expected vanishing Liouvillian gap, as shown in [52].
Conclusion-We have shown that constructing an FIG. 4. Liouvillian gap in the |k| = 1 sector. a) Gap vs g √ N . The colored vertical lines indicate g √ N values above threshold used to evaluate the trend of the gap closure with the system size shown in b). This is computed for N in the range 5 × 10 4 − 2 × 10 5 . The intercept shows the extrapolated gap as N → ∞. All parameters as for Fig. 3 atom-only effective theory for a model with U(1) symmetry presents surprising challenges. We have shown why standard 2 nd order Redfield theory fails to describe the superradiant transition, and we have shown how this can be rectified by the 4 th order terms of a diagrammatic expansion. We see this correctly describes the steady state and relaxation to that state in the k = 0 sector. Moreover, we may see how a semiclassical approximation becomes valid as N → ∞, through the emergence of an increasingly sharp distribution-such an approach may provide alternate ways to understand models where it is found that the semiclassical (mean-field) approximation fails [62, 63]. When considering spin coherences (i.e. the k = 1 sector), we surprisingly find that the 4KRE does not predict a vanishing Liouvillian gap in the symmetrybroken state. Moreover, we note that the 4KRE provides a Hermitian, trace-preserving, and secularized density matrix equation which is nevertheless not of Lindblad form. A key task for future work is to extend the methods developed here to a full multimode (e.g. confocal cavity) and spatially extended system. This would provide a powerful tool to theoretically explore non-equilibrium phase transitions in these driven-dissipative multimode systems.
Acknowledgments-We are grateful for helpful discussions with A. Daley, F. Damanet. R.P. wishes to thank G. Baio for stimulating discussions. We are grateful to B. Lev for comments on an earlier version of the manuscript. R.P. was supported by the EPSRC Scottish Doctoral Training Centre in Condensed Matter Physics (CM-CDT), Grant number: EP/L015110/1.  This section provides details of the fourth order Keldysh-Redfield equation (4KRE), using Keldysh diagrammatic perturbation theory introduced by Müller and Stace [14].
The quantum dynamics of the system density matrix, ρ(t), written in the interaction picture, is given by: where Σ(t 1 , t) is a self-energy superoperator. Müller and Stace [14] provide a diagrammatic recipe to derive this self energy perturbatively in the system bath coupling (which here means g). We use this approach up to fourth order to find Σ(t 1 , t), and then make a Markovian approximation ρ(t 1 ) ρ(t) to provide a time-local equation of motion for ρ. Because we consider evolution of a density matrix, the diagrammatic representation of the super-operators contributing to the self-energy consists of two horizontal solid lines representing free evolution of the system to the left and right of the density matrix (bottom and top branches respectively). An nth order diagram then contains n vertices corresponding to the action of H I at the times t 1 , ...t n−1 and at t. Dashed lines connecting these vertices represent pairwise correlations of the (Gaussian) bath, following Wick's theorem. Crucially, only irreducible diagrams are drawn, describing the contributions of nth order terms which cannot be factorised into products of lower order terms.
The irreducible diagrams shown in Fig. S1 are the fourth order contributions. (The second order contributions are the standard Redfield theory, and their contribution D 0 ρ + D 2 ρ is given explicitly in the Letter). In total there are 32 fourth order diagrams, 16 shown in Fig. S1, and an additional 16 diagrams obtained by swapping vertices between the lower and upper branches. The diagrams can be divided in three categories: (a) terms where all operators act to the left of ρ, (b) terms where operators appear in pairs to the left and the right, and (c) terms with three operators to the left and one to the right. As stated in the Letter, the form of our interaction Hamiltonian means that vertices which are connected must have opposite spin flip operators S ± .
Evaluating these diagrams, the 4th order contribution to the density matrix equation takes the form: where we have used bers of raising and lowering operators. In this section we explore whether this structure also preserves positivity, equivalent to requiring that the density matrix equation of motion takes Lindblad form [53]. We note however that, as discussed elsewhere [13,[55][56][57]64], there are many situations where a Redfield equation which does not preserve positivity may yet give an accurate description of the reduced system dynamics. To check for the Lindblad form we follow the method in Ref. [12] to extract the Lindblad-Kossakowski matrix. We describe here the approach we use to perform this numerically.
To construct the Lindblad-Kossakowski matrix we start by writing our equation of motion,ρ = Mρ, in matrix form:ρ The next step is to rewrite Eq. (S3) in terms of a complete set of linearly independent N × N matrices spanning the Hilbert space. It will be convenient to make use of two sets of basis matrices. The first is the element basis O i , with i = 0, ...N 2 − 1. These matrices are defined by [O i ] j,k = δ j,Q i,N δ k,R i,N where Q i,N and R i,N are respectively the quotient and the remainder of i modulo N . The second set is the normalized generalized Gell-Mann (gGM) basis, γ i , with i = 0, ...N 2 − 1.
These matrices are defined by T r(γ i γ j ) = δ ij . We include the identity matrix in this set as γ 0 = N / √ N . The next N matrices, p = 1 . . . N are diagonal matrices, γ p = diag(1, 1, . . . , −p, 0, 0, . . .)/ p(1 + p) (i.e. matrices with first p diagonal elements being one, the next element being −p, and all other elements zero). The remaining matrices are off-diagonal, of the form: labeled by ordered pairs of integers I > J. We use two bases as it is simplest to translate M to the element basis, but to find the Lindblad-Kossakowski matrix we need a basis with the identity matrix as one of the elements.
Using the above, we first write the density matrix equation in terms of the element basis aṡ Summing explicitly over the quotient and remainders of i, j, and defining the function I(Q, R) = QN + R, yields: Therefore, by comparison with Eq. (S3), we find: Now we have L O , the next step is to make a basis transformation to the gGM basis, . We thus have: where we have used γ † i = γ i . The matrix L γ thus has the form: Because the gGM basis contains the identity matrix as element 0, we may now use standard results [12] to write: where we have defined: We thus see that the submatrix of L γ excluding row and column 0 is the Lindblad-Kossakowski matrix. Extracting this numerically from our density matrix equation of motion, we find that the 4KRE is not of Lindblad form since not all eigenvalues are positive. We also note that independent of N , we find only 8 non-zero eigenvalues.
To conclude this section, we note that for second-order Redfield equations, one way to put the equation into Lindblad form is secularization [59]: eliminating from the Redfield equation those terms which are time dependent in the interaction picture. For our equation, the U (1) symmetry (and thus matching numbers of raising and lowering operators) automatically means that the equation is time-independent in the interaction picture. As such, it is notable that the 4KRE yields an equation for which secularization does not ensure positivity.

Block structure of density matrix equation
In this section, we analyze the block structure of Eq. (S2). As noted in the Letter the U(1) symmetry means the equation of motion breaks up into separate blocks. By writing we find that each k block evolves independently: We may thus also label eigenvalues and eigenvectors by blocks, L M,M takes a pentadiagonal form, which we may write as follows: Equivalently, in explicit matrix form, this means: (S14) The k = 0 sector addresses the behavior of the symmetry-preserving variables as S z and S z S z , which may have a non-zero steady state value. In contrast k = 0 describes coherences, which (for finite N ) must vanish at long times due to the U (1) symmetry. The spectrum of the k = 0 block thus provides details about the relaxation dynamics to the stationary states S z SS and S z S z SS .
The five functions A In the Letter we note that the Liouvillian gap, evaluated from the k = ±1 sector, does not vanish as expected at large N . To verify there is no similar issue in the k = 0 sector, Fig. S2 shows how, for g < g c the value of S z /N converges to −1/2 as N → ∞.

Probability distribution at large N
In the Letter we discussed how the form of the probability P M suggested the approach of writing coupled equations for the first two moments of the probability distribution. Here we present the results that support this statement. The probability distribution in the steady state is given by P M = R  Figure S3 shows the probability distribution for four different values of g √ N , and a range of values of N . Below threshold (panel a), for g < g c , the distribution is peaked at M = −N/2, with a finite width that vanishes as N increases, consistent with Fig. S2. In contrast, for g > g c the peak of the distribution moves away from M = −N/2, and asymptotically approaches a peak at M = 0 at large g √ N (panels c,d).
One clearly sees that for g > g c the probability distribution has a quasi-Gaussian form, with a width that shrinks as N increases. For g < g c , the probability distribution is one-sided, but can still be fit effectively by a Gaussian (peaked at S z < −N/2). These observations explain why the cumulant equation adequately describes the steady state at finite N , and how mean-field theory becomes valid in the thermodynamic limit.

Cumulant equations of motion
As we have seen, the steady-state probability distribution is approximately Gaussian, and so it should be possible to accurately model using its first two moments. We thus derive equations of motion for S z and S z S z . These equations depend on third order correlators of spin operators (from the second order terms in the density matrix equation) and fourth order correlators (from the fourth order terms). Such terms can be decoupled to second order terms using the cumulant expansion [60].
Because the cumulant expansion here is based on the approximate Gaussian distribution of P M , such an expansion applies exclusively to decoupling products of S z operators. As such, before decoupling the equations of motion, we first rewrite all expressions in terms of S z operators by using the identity S + S − = S (S + 1) − S z S z + S z . The U (1) symmetry of our model ensures that this will always be possible. We also note that we must keep track of non-commuting operators S ± , S z in this process. Once all terms are written in terms of only S z , these operators commute, and so can be treated analogously to classical cumulant expansions.
In the single-mode Dicke model, because there is no U(1) symmetry, there is no separation of the density matrix equation of motion into separate sectors. Therefore, we must numerically diagonalize the complete Liouvillian. Figure. S4 shows the corresponding Liouvillian gap-Liouvillian eigenvalue with the smallest (non-zero) real part. As shown in Fig. S4(a), the gap appears to close over as one enters the superradiant phase (pink region), even for 500 spins. We also show the evolution of gap with system size. It is noteworthy that the gap found this way scales differently to the U(1) problem discussed in the Letter. For the U(1) model, we found a scaling λ (1) 0 = A + B/N while for the Z 2 Dicke model, the behavior is λ 1 ∝ exp(−CN ), as shown in Fig. S4(b).