Superradiant to subradiant phase transition in the open system Dicke model: dark state cascades

Collectivity in ensembles of atoms gives rise to effects like super- and subradiance. While superradiance is well studied and experimentally accessible, subradiance remains elusive since it is difficult to track experimentally as well as theoretically. Here we present a new type of phase transition in the resonantly driven, open Dicke model that leads to a deterministic generation of subradiant states. At the transition the system switches from a predominantly superradiant to a predominantly subradiant state. Counterintuitively, the cavity decay is the crucial parameter for subradiant state generation and not the individualizing process of spontaneous decay. The observed effect is thus a cavity assisted generation of subradiant quantum coherences. Clear experimental signatures for the effect are presented and entanglement properties are discussed. Letting the system relax into the ground state generates a cascade of dark Dicke states, with dark state populations up to unity. Furthermore we introduce a collectivity measure that allows to quantify collective behaviour.


I. INTRODUCTION
The open (and closed) system Dicke model has been a work horse in quantum optics and beyond for decades  . Current research on Dicke model based systems includes novel laserlike systems 18 , phase transitions 15,22 , quantum information and super/subradiance 10,13,19,20,23 .
In recent years superradiance has been investigated with respect to entanglement 19 and subradiance for its prospects to store quantum information 20,23 . The Dicke model assumes N identical two-level systems, interacting with a bosonic cavity mode.
Investigating subradiant effects in a consistent open system theory was not feasible for a long time since in a straight forward approach the master equation scales exponentially in the number N of two-level systems. This renders full simulations even for small N impossible, however subradiance is a few and many particle effect. Existing limits and approximations for both analytical and numerical treatments addressing this problem are not suited to study subradiance for moderate N [4][5][6][7][8]10,13 . Usually for superradiance total spin conservation (explained below) is assumed, entirely neglecting subradiant states. This reduces the numerical complexity to ∼ N 2 or sometimes even allows analytic solutions 7,8,13 .
However ubiquitous phenomena in real systems like decay processes and pure dephasing break this symmetry. Therefore, both realistic treatments and subradiant effects require a different methodology.
Symmetries in the associated master equations reduce the complexity from an exponential scaling in N to a polynomial scaling ∼ N 3 , even without total spin conservation 18,22,24-26 . This makes exact calculations for moderate emitter numbers feasible and removes constraints imposed by assumptions and approximations. Furthermore the method can be applied to any permutation symmetric multi-level system setup 24 .
In this work we investigate the population of subradiant states through decay and pure dephasing processes -both do not conserve the total spin. Counterintuitively, the cavity lifetime determines the population of the subradiant states: Increasing the external driving results in a nonequilibrium phase transition and for short cavity lifetimes subradiant states are always suppressed by quantum coherence. Contrary increasing the cavity lifetime results in an amplification of subradiant states due to quantum coherence. Experimentally accessible signatures of this effect and entanglement properties via spin squeezing are discussed.
Switching off the external driving, the subsequent relaxation into the ground state forms

II. MODEL SYSTEM
We consider the usual Dicke model with an additional classical optical, cw field E driving all TLS identically. Driving is necessary since subradiant states are excited states. In a frame rotating at the external laser frequency, using the rotating wave approximation the system Hamiltonian reads where ∆ 0 , ∆ 1 are the mode and TLS detuning, g is the TLS-mode coupling, E is the optical driving, b, b † are photonic operators and J k = i σ i k , k = 11, 10, 01, 00 are the collective spin operators. Excited and ground state of the individual TLS i are |1 i , |0 i , the lower indices of the spin operators represent the Ket and Bra notation: We assume resonant excitation field, cavity and TLS. Both cavity and TLS are subject to loss and dephasing, using Lindblad formalism 27 . The master equation reads The Lindblad dissipators describe decay processes like individual radiative and non-radiative Fig. 1 (a). We use σ i z = σ i 11 − σ i 00 . All contributions to the master equation except D de and D pd are total spin preserving, ( Fig.   1 (b)). The total spin l(l + 1) is the eigenvalue of the J 2 = (J 10 J 01 + J 01 J 10 )/2 + J z operator, The value of l varies between l max = N/2 for the superradiant subspace and l min = 0, 1/2 for the (most) subradiant subspace. The J 2 and J z eigenvalues determine the coupling strength of the multi TLS (Dicke) state to an optical mode, the collective dipole transition element. This coupling strength distinguishes between superradiance and subradiance. For superradiant states the dipole element scales superlinear in N , while for subradiant states the scaling is sublinear in N and some subradiant states are dark 28 . Dark means that the dipole transition element of the collective excitation vanishes, meaning these states cannot decay e.g. creating a cavity photon. However these states still decay into other states via the decay and dephasing processes D de and D pd acting individually on the emitters, c.f. Fig. 1 (b). Generally the spin preserving contributions in the master equation In the bad cavity limit (κ g) equation (2) corresponds to the cooperative resonance fluorescence setup 4,5 . The system exhibits a non-equilibrium phase transition for increasing E for both total spin preserving and nonpreserving setups, where the nonpreserving setup was studied using mean field theory 4 . For longer cavity lifetimes κ the system more and more resembles the absorptive optical bistability setup 29 (instead of driving the TLS, in optical bistability the cavity is driven, opposed to Fig. 1 (a)). In the range investigated in this work (κ ∼ g) the clear distinction between cooperative resonance fluorescence and optical bistability breaks down, thus combining these distinct fields of quantum optics. Besides the steady state, density matrix states with very long lifetimes can exist in these systems, which lead to the observation of bistabilities in experiments with finite measurement time 30 .
In some limits these lifetimes go to infinity, resulting in a second steady state. For optical bistability these long lifetimes are called tunneling times 31,32 , more general this phenomenon is called dissipative phase transition 33 .

III. PERMUTATION SYMMETRIC METHOD
The permutation symmetry allows the incorporation of the individual TLS decay and  To quantify the effect of collectivity and distinguish between collective (super-and subradiance) and individual (dipole moment scales linear in N ) behavior we introduce the ratio between the full Dicke subspace population and its incoherent part We use γ = 1.0 ns −1 and g = 3.3 meV throughout this work. Please note that ultrastrong coupling effects are not present in the investigated parameter range. There are two types of dephasing/individualization processes: spontaneous decay and pure dephasing. We first investigate the spontaneous decay and investigate the effects of pure dephasing later.
Including small pure dephasing preserves all effects (see Section V B for a discussion).

A. Nature of the phase transition
In the steady state the most basic feature of the nonequilibrium phase transition is the change from the ground state to a half excited TLS state with increasing external driving field ( Fig. 2 (a)). Increasing the cavity quality (decreasing the ratio between cavity decay rate and TLS-cavity coupling strength κ/g) makes the transition sharper but the overall effect does not change much. Contrary a drastic change is seen in the behavior of the collectivity measure for the superradiant subspace R(l max = N/2), Fig. 2 Fig. 3 (c), we see that the increase due to The total occupation in the superradiant subspace goes to zero above the phase transition for N → ∞, Fig. 3  However by looking at both the absolute and relative populations we conclude that in the good cavity and large N limit the system changes from a predominantly superradiant to a predominantly subradiant state at the phase transition. This constitutes the main result of this work.
In Fig. 4 the scaling of experimentally more accessible quantities with the number of individual TLS N is presented: The normalized TLS excitation develops a kink for increasing  N , indicating a second-order transition, Fig. 4 (a). The smallest magnitude nonzero eigenvalue λ 1 of the Liouville operator L (c.f. equation (2)), which corresponds to the slowest time scale in the system to reach steady state, decreases around the phase transition for increasing N , Fig. 4 (b). It might even vanish for N → ∞, creating a second steady state. This could be measured for instance in a hysteresis cycle typical for optical bistability experiments 15,30,41 . The intracavity mean photon number shows the formation of a local minimum at the transition and an increase in the peak intensity, Fig. 4 (c). Also bunching (g (2) (0) > 1) increases for increasing N , Fig. 4 (d). Overall the transition becomes sharper and more pronounced for increasing N and decreasing κ/g, since these parameters increase the system size. This displays a typical property of phase transitions, which are well defined only in the thermodynamic limit (infinite system size) and blur for small system sizes 4,42,43 .

B. Robustness test and entanglement properties
So far all results were presented without including pure dephasing. Now we investigate the robustness of the collective effects at the phase transition against pure dephasing: In and subradiant state above phase transition is preserved for δ ∼ γ. The general trend of total Dicke subspace occupation is not affected by pure dephasing, as in Fig. 3 (d).
In the spin preserving setup the TLS are entangled via spin squeezing below the phase transition 13 . Spin squeezing is a concept originating from quantum metrology, where it was developed around the idea that squeezed atomic coherent states could be used for measurement precision below the shot noise limit, but also has attracted a lot of attention as an entanglement witness [44][45][46][47] . Here we employ the spin squeezing inequalities introduced  5 (b)). Hence the entanglement detected in the spin preserving setup is still present for spin nonpreserving setups and even for moderate pure dephasing times.

C. Dark state cascades
Super-and subradiance are concepts related to time evolution and so far we have only discussed the steady state: Now, we drive the system to the steady state with maximum R(l min ) (see Fig. 3 (b) and (c)) and then, afterwards, we switch off the driving field.
The system relaxes into the ground state and we observe that a cascade of dark states is generated, Fig. 5 Fig. 1  The parameters used in this study are realistic for NV centers, quantum dots and Rydberg atoms and the behavior is stable over a wide parameter range.
In summary we have shown that the nonequilibrium phase transition of cooperative resonance fluorescence changes drastically when leaving the bad cavity limit: Subradiant Dicke states are amplified and clear experimental signatures of this effect emerge. Letting the system relax into the ground state generates a dark state cascade that can be utilized to store quantum information.

VII. ACKNOWLEDGEMENTS
We thank Nicolas Naumann for useful discussions, and gratefully acknowledge support

Appendix B: Spin sqeezing inequalities
We employ the spin squeezing inequalities (SSI) introduced by Tóth et al. 48,49 as entanglement measure. Tóth et al. derived seven inequalities that are satisfied by any separable N -qubit state, hence the violation of any of these inequalities implies entanglement. Four of the seven inequalities detect entanglement in our setup, but the violation of two equations is equivalent: the coherent driving field introduces a time dependent phase factor caused by local unitary transformations which do not affect entanglement 56 but cause the violation of the SSI to oscillate back and forth between the two associated inequalities (between (B1), where the variances are defined as (∆A) 2 = A 2 − A 2 . In order to simplify the discussion we only show one SSI in our plot: hence A is the quantity plotted in Fig. 5 (b). Since strictly speaking the quantities J 2 y and J 2 x do not have a defined steady state, but oscillate with the phase factor mentioned above, we set t = 0 and thus set the phase factor to unity throughout the plot in Fig. 5 (d).
Since, as stated above, the local unitary transformations causing the oscillation do not affect the entanglement, this is a valid approach. In the following the local unitary transformation is explained: On resonance the Hamiltonian of the system in a frame rotating at the external laser frequency ω l reads H = g(J 10 b + J 01 b † ) + E(J 10 + J 01 ). (B6) The corresponding master equation for the setup considered in this work is where ρ is the rotating frame density matrix. The transformation between normal frame and rotating frame is given by with the normal frame density matrix ρ n and the Hamiltonian The Hamiltonian acts locally on the density matrix, in the sense that each TLS experiences an individual unitary transformation, i.e.
tities arising in the SSI experience a time dependency through this transformation. In fact only the rotating frame density matrix has a stationary steady state, the normal frame density matrix ρ n exhibits an oscillating steady state, where diagonal entries are stationary and offdiagonal entries oscillate with a phase of multiples of ω l .
The quantities J 2 x,y and (∆J x,y ) 2 are explicitly time dependent in the normal frame. By adding Eqs. (B1), (B2) and (B3), (B4) respectively, one can derive time independent inequalities, which however do not detect entanglement in our setup.