Table Of Content

Collective phases of strongly interacting cavity photons Ryan M. Wilson1, Khan W. Mahmud2, Anzi Hu3, Alexey V. Gorshkov2,4, Mohammad Hafezi2,5,6, and Michael Foss-Feig2,7 1Department of Physics, The United States Naval Academy, Annapolis, MD 21402, USA 2Joint Quantum Institute, NIST/University of Maryland, College Park, MD 20742, USA 3Department of Physics, American University, Washington, DC 20016, USA 4Joint Center for Quantum Information and Computer Science, NIST/University of Maryland, College Park, Maryland 20742, USA 5Department of Electrical Engineering and Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, MD 20742, USA 6Kavli Institute of Theoretical Physics, Santa Barbara, CA 93106, USA and 7United States Army Research Laboratory, Adelphi, MD 20783, USA 6 We study a coupled array of coherently driven photonic cavities, which maps onto a driven- 1 dissipative XY spin-1 model with ferromagnetic couplings in the limit of strong optical nonlineari- 2 0 ties. Using a site-decoupled mean-field approximation, we identify steady state phases with canted 2 antiferromagneticorder,inadditiontolimitcyclephases,whereoscillatorydynamicspersistindef- p initely. We also identify collective bistable phases, where the system supports two steady states e among spatially uniform, antiferromagnetic, and limit cycle phases. We compare these mean-field S resultstoexactquantumtrajectoriessimulationsforfiniteone-dimensionalarrays. Theexactresults exhibit short-range antiferromagnetic order for parameters that have significant overlap with the 6 mean-fieldphasediagram. Inthemean-fieldbistableregime,theexactquantumdynamicsexhibits 1 real-time collective switching between macroscopically distinguishable states. We present a clear physical picture for this dynamics, and establish a simple relationship between the switching times ] h and properties of the quantum Liouvillian. p - nt Despite numerous outstanding questions, the study of acting systems in mind. For example, it is not fully un- a quantum many-body systems in thermal equilibrium is derstoodhowthesteadystatesofthesesystemsrelateto u on relatively solid ground. In particular, very general the equilibrium states of their “closed” counterparts, or q guiding principles help to categorize the possible equi- how conventional optical phenomena, such as bistability, [ librium phases of matter, and predict in what situations manifests in the presence of strong optical nonlinearities 2 they can occur [1–3]. In comparison, quantum many- and spatial degrees of freedom. v bodysystemsthatarefarfromequilibriumarelessthor- Weconsideranarrayofcoupled,single-modephotonic 7 oughly understood, motivating a large scale effort to ex- cavities described by a driven-dissipative Bose-Hubbard 5 plore non-equilibrium dynamics experimentally, in par- model [31–35], which maps onto a driven-dissipative XY 8 ticular using atoms, molecules, and photons [4–9]. At spin-1 model in the limit of strong optical nonlinearity. 6 2 0 the same time, it has become clear that studying non- Weperformacomprehensivemean-field(MF)study,and . equilibrium physics in these systems is often more natu- identify a variety of interesting steady states including 1 ral than studying equilibrium physics; they are, in gen- spin density waves and limit cycles, which break the dis- 0 eral,intrinsically non-equilibrium. Forexample,thermal crete translational symmetry of the system. The spin 6 1 equilibrium is essentially never a reasonable assumption density waves possess canted antiferromagnetic order for : inphotonicsystems,wheredissipationmustbecountered a range of drive strengths, despite the ferromagnetic na- v by active pumping [10]. Indeed, the inadequacy of equi- ture of the spin couplings. Interestingly, the exact quan- i X librium descriptions for photonic systems has long been tum solutions exhibit short-range antiferromagnetic cor- r recognized [11], even though close analogies to thermal relations for parameters that have notable overlap with a systems sometimes exist [12–16]. the MF results. The system also supports collective bistable phases, which manifest in the exact quantum Until recently, photonic systems have been restricted dynamics as fluctuation-induced collective switching be- to a weakly interacting regime. With notable progress tween MF-like states. We present a simple relationship towards generating strong optical nonlinearities at the between this dynamics and properties of the quantum few-photon level, for example with atoms coupled to Liouvillian. small-mode-volume optical devices [17–22], Rydberg po- laritons [23, 24], and circuit-QED devices [25–28], this situationisrapidlychanging. Theproductionofstrongly interacting, driven and dissipative gases of photons ap- I. MODEL pearstobefeasible[29,30],andaffordsexcitingopportu- nitiestoexplorethepropertiesofopenquantumsystems For a system weakly coupled to a Markovian environ- in unique contexts, while studying the applicability of ment,thedynamicsofitsdensitymatrixρîsgovernedby theoretical treatments designed with more weakly inter- amasterequation∂ ρˆ=L [ρˆ],whereL [ρˆ]=−i[Hˆ,ρˆ]+ t 2 trivial steady states ρˆss, which satisfy L[ρˆss] = 0. Gen- erally, these non-equilibrium steady states are qualita- tively distinct from the equilibrium states of the closed (γ = Ω = 0) BH model described by Hˆ , which are BH characterizedbyasuperfluidorderparameterthatspon- taneously breaks the U(1) symmetry associated with particle number conservation. This U(1) symmetry is explicitly broken by the coherent laser drive, and the driven-dissipative BH (DDBH) model conserves neither energynorparticlenumber. Thereforesuperfluiditycan- not emerge in the DDBH model; this is in contrast to similar models with incoherent pumps, which can support superfluid phases [39–45]. The DDBH model does, however, possess a spatial symmetry generated by dis- crete translations along the Bravais vectors of the cavity FIG. 1. (color online). Mean-field phase diagram for array, which is broken spontaneously if spatial structure J/γ = 10 in the hard-core (U → ∞) limit. Represented are develops in the steady state. dark (U ) and bright (U ) uniform states, canted antiferro- 1 2 We study the steady states of the DDBH model us- magnetic states (AF), frustrated AF states (f-AF), and limit ing both a site-decoupled Gutzwiller mean-field (MF) cycles (LC). Regions with double labeling exhibit bistability method [46] and exact quantum trajectory simulations between the indicated states. Discontinuous transitions are indicated by thick black lines. The inset shows a limit cycle of finite systems [47–50]. In this MF approximation, forµ/γ =2.5andΩ/γ =6.5projectedontotheBlochsphere. the density matrix is decomposed as a matrix product (cid:78) Thered(blue)linesrepresentthedynamicsofcavitiesinthe ρˆ = ρˆ, where ρˆ is the local density matrix at cav- i i A (B) sublattice. ity i. Further, we restrict our study to the “hard-core” limit, where strong optical nonlinearities produce a per- fect photon blockade, by taking U → ∞ [33]. In this D[ρˆ] is the Liouvillian, Hˆ is the system Hamiltonian, limit, the photons can be mapped onto spins by an in- and D[ρˆ] is a dissipator in the Lindblad form [36, 37]. verseHolstein-Primakofftransformation[51],resultingin Here,weconsideranarrayofcoherentlycoupled,nonlin- an effective driven-dissipative XY spin-1 model [40, 52] ear, single-mode photonic cavities driven by a spatially with symmetric, ferromagnetic spin coup2lings, described uniformlaserfieldwithfrequencyωl,whichleakphotons by the Hamiltonian Hˆ = −J (cid:80) (cid:0)σˆxσˆx+σˆyσˆy(cid:1) + intotheenvironmentatarateγ. Intheframerotatingat Ω(cid:80) σˆx− µ(cid:80) σˆz, where σˆx,4yd,z ar(cid:104)ei,jt(cid:105)heiPajuli miatrjices. the driving frequency, the Hamiltonian for the coherent i i 2 i i i We derive equations of motion for the spin components drive is Hˆl =Ω(cid:80)i(aî+aˆ†i). The system Hamiltonian is σα = Tr[ρˆσˆα]; these are given by Eqs. (A4) in the ap- then Hˆ =HˆBH+Hˆl, where piendix. Foriclarity,wespecializetothecaseofJ/γ =10, though many of the qualitative features discussed below Hˆ =−J (cid:88)aˆ†aˆ −µ(cid:88)nˆ + U (cid:88)nˆ (nˆ −1) (1) are valid more generally. BH d i j i 2 i i (cid:104)i,j(cid:105) i i is the Bose-Hubbard (BH) Hamiltonian [38]. Here, aˆ† II. MEAN-FIELD PHASE DIAGRAM i (aˆ ) creates (annihilates) a photon in cavity i, nˆ =aˆ†aˆ i i i i isthelocalnumberoperator,andthesumsrunfromi=1 By solving the MF equations for a variety of parame- to N, the number of cavities. The notation (cid:104)i,j(cid:105) implies ters,wefindheuristicallythatallsteadystatesareeither a further summation over all cavities j that are nearest- spatially uniform, where all spins point in the same di- neighbors to cavity i, the number of which is given by rection (σα = σα for all i,j), or have antiferromagnetic i j the coordination number z. The nonlinearity of the cav- (AF) spin density wave order, where neighboring spins ities is quantified by effective two-photon interactions of point in different directions, but next nearest neighbors strength U. The first term in Eq. (1) describes the hop- point in the same direction (σα (cid:54)= σα and σα = σα i i±1 i i±2 ping of photons between cavities; J is the hopping rate, for all i). Because the neighboring spins are not antipar- and d is the dimensionality of the system (d = z/2 for allel, this AF order is “canted.” This motivates the use hypercubic arrays). The laser drives the system with of a two-sublattice ansatz, which we solve by evolving strength Ω, and is detuned from the cavity resonance the MF equations of motion on two sites. We find a frequency ωc by µ = ωl −ωc. The coupling of the sys- variety of interesting steady state phases, shown by the tem to the environment is described by the dissipator colored regions in Fig. 1. In the blue region, the sys- D[ρˆ]= γ (cid:80) (2aˆ ρâˆ†−ρˆnˆ −nˆ ρˆ). temexhibitsbistabilitybetweenspatiallyuniformdarker 2 i i i i i For Ω = 0, the system evolves into a trivial vac- (U ) and brighter (U ) steady states. In the red region, 1 2 uum at long times, with (cid:104)nˆ (cid:105) = 0 for all i. The laser there is a unique steady state with AF order. In the i drive provides a photon source, and can stabilize non- green region, the system is bistable between U and AF 1 3 steady states. All phase boundaries in Fig. 1 correspond tocontinuoustransitionsexceptthoseatthethresholdof bistability(darklines),wheretheadditionalsteadystate appears discontinuously. The bistability in this system is inherently collective, in that it does not exist for a single cavity in the hard- core limit [53]. We note that collective bistability exists in a variety of other driven and dissipative systems [54– 62], andwasrecentlyobservedinangaseousensembleof Rydberg atoms [63]. Gases of Rydberg atoms are also predicted to exhibit AF order [64–66], though unlike the modelweconsiderhere,theirinteractions(duetotheRy- dberg blockade) are effectively antiferromagnetic in nature. Other works studying the hard-core DDBH model also predict AF order, though they consider variants of the model that include spatially varying drive fields [67], two-cavity pumping [52], and cross-Kerr terms [68, 69]. OursystemexhibitsAForderintheabsenceofthesefea- tures, despite the ferromagnetic nature of the couplings. IntheregionenclosedbythedashedlineinFig.1, the AFstatesacquirealimitcycle(LC)character[64,70,71], exhibiting periodic oscillatory dynamics at long times, andthusbreakingthecontinuoustime-translationalsym- metry of the system. The inset in Fig. 1 shows an example limit cycle projected onto a Bloch sphere. Interest- FIG. 2. (color online). Both panels correspond to N = ingly,thelimitcycleexistsbothasauniquesteadystate, 12 cavities with periodic boundary conditions and nearest- neighbor couplings for J/γ = 10. The black lines show the and as part of a U /LC bistable pair (green region). We 1 mean-field phase diagram boundaries. (a) Nearest-neighbor note that similar limit cycles were recently predicted for part of the σˆy correlation function, Σy. The red region a driven-dissipative XYZ spin-1 model [72]. i 1 2 indicates antiferromagnetic correlations. (b) δN2/N (see To determine the validity of the two-sublattice ansatz, text),whichisstronglyenhancedinthepresenceofcollective we perform a linear stability analysis on the spatially bistable switching. The inset in (a) shows Σy for µ/γ = 10 i uniform steady states [34, 35, 52] (see Sec. A2 in the and Ω/γ = 6 (12), exhibiting short-range incommensurate appendix). In the light blue and light green regions spin density wave (antiferromagnetic) order. in Fig. 1, the U steady state is dynamically unstable 1 to the formation of incommensurate (k < π) spin density waves, which cannot be captured within the two- ence of quantum and classical fluctuations, which exist sublattice ansatz. To understand the effects of this in- in the true steady state of the DDBH model, and play stability, we solve the inhomogeneous MF equations for particularly important roles in low dimensions. Toward a 1D chain, with randomly seeded initial states. In the this end, we employ a quantum trajectories algorithm to light blue region, the instability results in the disappear- study finite 1D systems [47–50] (see Sec. B in the ap- ance of the U1 steady state. In the light green region, pendix). This method provides exact results for physical we find steady states that exhibit AF order over finite observables in the steady state under the ensemble aver- size domains, which are interrupted by domains of the aging of trajectories. U state; we refer to this as frustrated AF order (f-AF). 1 We present results for a system with N = 12 cavi- WhiletheAFphaseistheonlytrueperiodicallyordered tiesandperiodicboundaryconditionsinFig.2. Inpanel steady state in this region, U1 domains remain stable if (a), we show the nearest-neighbor part of the σˆy con- they are sufficiently small, and do not sample unstable nected correlation function Σy, where Σy = (cid:104)σˆyσˆy (cid:105)− wave vectors that often exist over a very narrow range (cid:104)σˆy(cid:105)(cid:104)σˆy (cid:105), which is indepen1dent of j. i We cjhojo+sei to of k. The limit cycles remain mostly AF ordered, with j j+i study σˆy correlations because the y-components of the some low amplitude, small-k features, and the U state 1 spins exhibit the strongest AF order in the MF re- of the U1/LC pair ceases to exist in the light green LC sults, and Σy can be measured via correlated homo- region. i dyne detection in an experimental setup. The blue region shows where Σy is positive, and the red region 1 showswhereΣy isnegative,correspondingtoAFnearest- 1 III. BEYOND MEAN-FIELD neighbor correlations. The inset shows Σy for µ/γ = i 10 and Ω/γ = 6 (12), represented by green diamonds Itisimportanttounderstandhowtherichphysicspre- (red squares). While both correlation functions exhibit dicted by the Gutzwiller MF theory survives in the pres- nearest-neighbor AF order, the correlations for Ω/γ =6 4 FIG.3. (coloronline). AllpanelscorrespondtoN =12cav- FIG. 4. (color online). (a) Liouvillian gap Γ/γ for N = 20 ities and µ/γ = −2.5. (a) Example quantum trajectory for cavitieswithJ/γ =10andinfinite-rangecouplinginthehard- Ω/γ =2.5, in the bistable regime, showing collective switch- core (U → ∞) limit. The solid black line shows the mean- ingbetweendark(white)andbright(blue)states. (b)Mean- fieldbistablephaseboundaryforspatiallyuniformstates. (b) field calculation of total photon number, showing collective Liouvillian gap for µ/γ = −5. The solid dark (light) blue bistability. The solid (dashed) black lines correspond to sta- lineshowsthegapcalculatedbydiagonalizingtheLiouvillian ble(unstable)solutionsofthemean-fieldequations. (c)Total for N = 12 (8) cavities. The circles (squares) show the gap photon number for the quantum trajectory shown in (a). extracted from quantum trajectories simulations (see text); the error bars represent standard error. have incommensurate spin density wave character while cated by the grey markers. We show N for the quantum the correlations for Ω/γ = 12 have a true AF character, trajectory as a function of time in panel (c). Here, N switchingfromnegativetopositivevalueswitheachsuc- fluctuates about two mean values for extended periods cessive cavity. The incommensurate spin density wave of time, which are interrupted by switching events that correlations are present where the MF linear stability drive the system from one MF-like state to the other. analysis predicts an incommensurate spin density wave instability of the U steady state, corresponding to the 1 f-AF region in Fig. 1. IV. COLLECTIVE BISTABILITY & We performed finite size scaling calculations for these SWITCHING and a number of other parameters, and find that N = 12 accurately captures AF order of the DDBH model in A recent work [73] showed that collective bistabil- the thermodynamic limit for a large region of the phase ity in the driven-dissipative XY spin-1 model vanishes diagram; this is due to the short-range, exponentially 2 as the Gutzwiller approximation is systematically im- decaying nature of Σy. The black lines overlaid in Fig. 2 i proved. Our results demonstrate that while ρˆss is indeed show the MF phase diagram boundaries. Interestingly, unique[74,75],itretainsclearsignaturesoftheMFbista- thereisreasonablygoodagreementbetweentheregionof bility, namely that the system dynamically switches be- short-range AF correlations in the 1D quantum system tweentwomacroscopicallydistinguishableconfigurations and the MF results. thatresembletheMFsteadystatesρˆ1andρˆ2;infact,this In Fig. 2(b), we show the normalized fluctuations in isagenericfeatureofbistablesystems[60,61,65,76,77]. the total photon number Nˆ =(cid:80)Ni=1nî, given by δN2/N The approach to ρˆss in the bistable regime is then char- where δN2 = (cid:104)Nˆ2(cid:105)−(cid:104)Nˆ(cid:105)2 and N = (cid:104)Nˆ(cid:105). Interestingly acterized by switching between these two states, and the δN2/N becomesanomalouslylargeinthedarkerbluere- rateofconvergenceisdirectlyrelatedtotheaveragetime gion, which has significant overlap with the MF bistabil- spent in ρˆ1 and ρˆ2 between switching events; we refer to ity, indicating that photon number fluctuations become thesetimesasτ andτ ,respectively. Inthemasterequa- 1 2 strongly correlated when the MF theory predicts collec- tion formalism, the asymptotic rate of convergence is set tivebistability. TheoriginoftheδN2/N enhancementis by the Liouvillian gap Γ = −Re[ε ], where ε is the ex ex revealed upon inspection of the trajectories themselves. eigenvalue of L with the smallest magnitude non-zero InFig.3(a),weshowthephotonnumberasafunctionof real part [78, 79]; this suggests that Γ and τ are inti- 1,2 timeforN =12cavitieswithµ/γ =−2.5andΩ/γ =2.5, mately related [80]. The presence of true bistability in in the region of enhanced δN2/N (shown by the white the MF theory reflects the fact that the Liouvillian gap circle in Fig. 2(b)). Interestingly, the trajectory exhibits vanishes at this level of approximation. The collective collectiveswitchingbetweenmacroscopicallydistinguish- switching in the exact dynamics, and the concomitant ablestates,whichresembletheMFsteadystatesforthese opening of the Liouvillian gap, can be understood to re- parameters. We plot the total photon number obtained sult from both quantum and dissipation-induced (classi- via MF calculations as a function of Ω/γ for µ/γ =−2.5 cal) fluctuations that are not included in the MF theory. in panel (b) of Fig. 3; the values at Ω/γ = 2.5 are indi- It is natural to expect that MF behavior can be re- 5 covered in the quantum system as its coordination num- ACKNOWLEDGMENTS ber z is increased. We explore this possibility by taking the spin couplings (or photon hopping) to be infinite- WethankMohammadMaghrebi,SarangGopalakrish- range. This corresponds to modifying the hopping term in Eq. (1) to −J (cid:80) aˆ†aˆ , where J ≡ 2J/(N − 1). nan, and Jeremy Young for insightful discussions. RW, i(cid:54)=j i j KM, AG, and MH thank the KITP for hospitality. We While this limit may seem unnatural for arrays of pho- acknowledge partial support from the NSF, ONR, ARO, toniccavities,itcouldbeachievedusinganexternalmir- ARL, AFOSR, NSF PIF, NSF PFC at the JQI, and the ror,anditisinfactquitenaturalforotheropenquantum Sloan Foundation. This research was supported in part systems,suchasensemblesofRydbergatoms[64,66,81] by the NSF under Grant No. NSF PHY11-25915. or trapped ions [82, 83]. We calculate the Liouvillian gap exactly by taking advantage of an efficient param- eterization of the accessible space of density matrices (see [57, 84] and Sec. C in the appendix); we show the Liouvillian gap for N =20 as a function of µ/γ and Ω/γ Appendix A: Equations of motion in the U →∞ limit in Fig. 4(a). There is striking agreement between the exact quantum results and the MF results; where MF In this appendix, we provide technical details to sup- theory predicts collective bistability, the Liouvillian gap port the theory and numerics in the main text. Though decreases to Γ(cid:28)γ [57]. many of our results are valid more generally for soft- We plot Γ/γ for N =8 (12) cavities and µ/γ =−5 as core bosons, here we specialize to the limit of hard-core a function of Ω/γ in Fig. 4(b), shown by the light (dark) bosons, valid when U → ∞, where U is the local inter- blue solid lines. Savage and Carmichael proposed a two- action strength in Eq. (2) of the main text. Hard-core statetoymodeltodescribeabistablesystemwithasmall bosons can be conveniently mapped onto spins via an Liouvillian gap [85], which has a gap Γtoy = τ1−1+τ2−2. inverse Holstein-Primakoff transformation, We extract values for τ and τ heuristically by measur- 1 2 ingthetimespentinthedark(1)andbright(2)statesof the quantum trajectory simulations, and plot τ1−1+τ2−1 aˆ†i →σî+ and aî →σî−, (A1) in Fig. 4(b), shown by the squares (circles) for N = 8 (12). Already at N = 12, there is excellent quantita- where σˆ± = (σˆx±iσˆy)/2 and σˆ+σˆ− = (σˆz +1)/2, 1 is tive agreement between the exact Liouvillian gap and i i i i i i theSU(2)identityoperator,andσˆx,y,z arePaulimatrices the results of the simple two state model as extracted i that act on cavity i. Following this transformation, the from the quantum trajectories with infinite-range cou- master equation becomes plings. ThisprovidesaclearconnectionbetweentheMF and quantum solutions in the bistable regime. We note thatKinslerandDrummondperformedasimilaranalysis ∂ ρˆ=−i(cid:2)Hˆ,ρˆ(cid:3)+ γ (cid:88)(cid:0)2σˆ−ρˆσˆ+−σˆ+σˆ−ρˆ−ρˆσˆ+σˆ−(cid:1), for the single-mode quantum parametric oscillator, and t 2 i i i i i i i also found good quantitative agreement for large photon (A2) numbers [80]. where Hˆ is the system Hamiltonian, Hˆ =−J (cid:88)(cid:0)σˆxσˆx+σˆyσˆy(cid:1)+Ω(cid:88)σˆx− µ(cid:88)σˆz. V. DISCUSSION 2z i j i j i 2 i (cid:104)i,j(cid:105) i i (A3) Thethermodynamiclimitofthe1DDDBHmodelwith nearest-neighbor interactions is challenging to study nu- Here the summations run over i=1,...,N, and the no- merically, but we expect its dynamics to exhibit collec- tation (cid:104)i,j(cid:105) indicates an additional sum over all cavities tive switching over finite-size domains. This behavior is j that are coupled to cavity i; the number of cavities reminiscent of equilibrium systems that, while exhibit- j are quantified by the coordination number z. Equa- ing a first-order phase transition in higher dimensions, tions (A2) and (A3) describe a driven-dissipative spin-1 failtodosoin1D.Whetherornotmean-fieldbistability 2 XYmodel,withisotropicinteractionsandinthepresence is associated with a true first-order phase transition in of a homogeneous applied field. higher spatial dimensions is an interesting question that warrants further study. Finally, we note that we have To study the steady states of Eq. (A2), it is conve- studied the soft-core DDBH model (with finite U), and nienttouseequationsofmotionforthespincomponents identified features that are directly analogous to those σx,y,z = (cid:104)σˆx,y,z(cid:105); these equations are readily derived by i i discussed here. A comprehensive study of the soft-core taking ∂tσix,y,z = Tr[σîx,y,z∂tρˆ]. Using the Pauli ma- DDBH model is the subject of future work. trix commutation relations [σˆα,σˆβ]=2iδ ε σˆγ where i j ij αβγ i 6 ε is the Levi-Civita symbol, we find form steady states within the MF approximation, which αβγ is carried out as follows. Spatially uniform states have ∂ σx =−2J (cid:88)(cid:104)σˆzσˆy(cid:105)+µσy− γσx the property σiα = σjα ≡ σα for all i,j. Using this as an t i z i j i 2 i ansatz for Eqs. (A5), we find the following equations for (cid:104)i,j(cid:105) the uniform steady states, 2J (cid:88) γ ∂ σy = (cid:104)σˆzσˆx(cid:105)−µσx−2Ωσz− σy t i z i j i i 2 i γ (cid:104)i,j(cid:105) 0=−2Jσzσy+µσy− σx 2 2J (cid:88) ∂tσiz = z (cid:104)σîxσˆjy−σîyσˆjx(cid:105)+2Ωσiy−γ(σiz+1). 0=2Jσzσx−µσx−2Ωσz− γσy 2 (cid:104)i,j(cid:105) (A4) 0=2Ωσy−γ(σz+1), (A6) which can be solved analytically to find the spin config- 1. Gutzwiller mean-field approximation urations of the uniform steady states. For example, we plot solutions for J/γ = 10 and µ/γ = −5 as a func- In the Gutzwiller mean-field (MF) approximation, the tionofΩ/γ inFig.5(a);thesesolutionsexhibitcollective density matrix is assumed to factorize over all cavities, bistability, wherethedynamicallystablesteadystateso- (cid:78) ρˆ= iρî. Inthespinformalism,thiscorrespondstothe lutions are shown by the solid black lines. The dashed factorizationofallnon-localtwo-spinexpectationvalues. black line in this figure indicates a dynamically unstable Thus, in the MF approximation, Eqs. (A4) become solution, which we discuss in more detail below. The ansatz of spatially uniform steady states is not 2J (cid:88) γ ∂ σx =− σz σy+µσy− σx appropriate under all circumstances, for example, if the t i z i j i 2 i (cid:104)i,j(cid:105) steady state develops a spatial order that breaks the dis- 2J (cid:88) γ crete translational symmetry of the DDBH model. The ∂tσiy = z σiz σjx−µσix−2Ωσiz− 2σiy presence of such spatially ordered states can sometimes (cid:104)i,j(cid:105) becapturedbyaninstabilityoftheuniformsteadystate. 2J (cid:88) 2J (cid:88) To explore this possibility, we perform a linear stabil- ∂ σz = σx σy− σy σx+2Ωσy−γ(σz+1). t i z i j z i j i i ity analysis on the uniform steady state solutions of (cid:104)i,j(cid:105) (cid:104)i,j(cid:105) Eqs. (A6), specializing to one-dimensional systems with (A5) z =2. Inasystemwithinfinitespatialextent, thisisac- complishedbyaddingasmallplanewaveperturbationto InthemajorityofthephasediagraminFig.1ofthemain the uniform steady state of the form σ = σ +δeikm, m 0 text, the steady states are obtained by evolving these where σ = (σx,σy,σz)T and σα are the solutions of equations numerically using a 4th order Runge-Kutta al- 0 0 0 0 0 Eqs. (A6), m are cavity indices, and k is the wave num- gorithm. ber of the perturbation. Linearizing in δ, we find the equation 2. Linear stability analysis ∂ δ =Mδ, (A7) t ThephasediagraminFig.1ofthemaintextalsoshows results from a linear stability analysis on spatially uni- where  −γ µ−2Jσzcos(k) −2Jσy  2 0 0 M= 2Jσ0zcos(k)−µ −γ2 2Jσ0x−2Ω. (A8) 2Jσy(1−cos(k)) −2Jσx(1−cos(k))+2Ω −γ 0 0 ThematrixMhaseigenvaluesωthatdependonthewave ics of the instability. For example, in Fig. 5(b) we plot number k. When the real part of ω is negative for all k, Re[ω/γ] as a function of k for three sets of parameters, theuniformsteadystateσ isdynamicallystable. When all with J/γ = 10. The blue lines correspond to one of 0 some ω acquires a positive real part, this signifies an in- three solutions with Ω/γ = 2, shown by the blue circle stability of the uniform steady state, and the wave num- in panel (b) of this figure. These modes are unstable beratwhichRe[ω]ismaximumcorrespondstothemode at k = 0, indicating a global instability of the uniform that is maximally unstable, and dominates the dynam- steadystatesolution. Thegreenandredlinescorrespond 7 0.5 a) b) µ/γ=−5,Ω/γ=2 10 µ/γ=10,Ω/γ=12 0.4 µ/γ=10,Ω/γ=6 5 ×10 +1)/20.3 γ[ω/] 0 zσ00.2 Re ( -5 0.1 -10 0 0 1 2 3 4 5 6 7 8 0 1/4 1/2 3/4 1 Ω/γ k/π FIG. 5. a) Gutzwiller mean-field result for (σz+1)/2 (equivalent to local photon number) for J/γ =10 and µ/γ =−5. The 0 solid black lines represent steady state solutions that are stable at k = 0; the dashed black line represents a solution of the mean-fieldequations thatis unstableatk=0. b) Realpart ofthe linearstabilityeigenvalues asa functionof wave numberk. The values of k for which Re[ω/γ]>0 indicates an instability of the uniform steady state with these values of k. to parameters exhibiting spin density wave instabilities. The norm of |ψ(cid:105) is not conserved in real-time evolution The green line exhibits an incommensurate spin density under Hˆ due to its being non-Hermetian, so in general eff wave instability, as the modes are stable at k = π but we have unstableinasmallregionofk <π. Theredlineexhibits anantiferromagnetic(AF)spindensityinstability,asthe (cid:104)ψ(cid:48)(t+dt)|ψ(cid:48)(t+dt)(cid:105) =1−δp. (B4) most unstable eigenvalue occurs for k =π. (cid:104)ψ(t)|ψ(t)(cid:105) and antiferromagnetic (AF) instabilities, respectively. In the quantum trajectories formalism, we choose the wave function at t+δt stochastically. With probability 1−δp, we choose Appendix B: Exact numerical solution of Eq. (A2) |ψ(cid:48)(t+δt)(cid:105) |ψ(t+δt)(cid:105)= √ , (B5) We employ a quantum trajectories algorithm to solve 1−δp the master equation (A2), which is a powerful, exact method for studying open quantum systems that relies and with probability δp we take our next state to be one on treating the system-environment coupling (the dissi- that emitted a photon from cavity i, pator in the master equation) stochastically. Here, we aˆ |ψ(t)(cid:105) present the algorithm for generating a quantum trajec- |ψ(t+δt)(cid:105)= i , (B6) (cid:112) tory using quantum jumps. We begin by defining the δp /δt i effective non-Hermetian Hamiltonian, where the jump-site i is chosen with probability Hˆ =Hˆ−iγ (cid:88)nˆ , (B1) eff 2 i i Πi = δδppi = δt(cid:104)ψ(t)δ|npî|ψ(t)(cid:105), (B7) so Eq. (A2) can be written as (cid:80) and δp =δp. i i ∂ ρˆ=−i(cid:16)Hˆ ρˆ−ρˆHˆ† (cid:17)+γ(cid:88)aˆ ρâˆ†. (B2) t eff eff i i i Appendix C: Infinite range interactions The quantum trajectories method amounts to evolving In general, long-range interactions introduce consider- a wave function, as opposed to a density matrix, under Hˆ , while treating the right-most “recycling” term in able difficulties when numerically studying the dynam- eff ics of many-body systems. However, when all pairs of Eq. (B2) stochastically. spinsinteractwiththesamestrength,themodelbecomes Consider a wave function |ψ(t)(cid:105) at time t that evolves symmetric under the exchange of any two spins, which under Hˆ . For a small time step δt, the wave function eff allows for an extremely efficient parametrization of the at t+δt can be written as (using Euler integration) accessible Hilbert space. In the absence of spontaneous emission, permutation symmetry of the Hamiltonian re- |ψ(cid:48)(t+dt)(cid:105)=(1−iHˆ dt)|ψ(t)(cid:105). (B3) eff stricts the dynamics of initially permutation symmetric 8 pure states to the subspace of collective spin-states, of- permutation symmetric density matrices, which can be ten referred to as Dicke states, and the dimension of parametrized in terms of symmetrized direct products of the accessible Hilbert space is therefore O(N). For an Pauli matrices. The space of permutation symmetric, open quantum system governed by a permutation sym- Hermitian matrices over the product Liouville space of metric Liouvillian, dynamics is restricted to the space of N spins is spanned by basis states of the following form, Mˆ(n)= 1 (cid:88)(cid:16)σˆx ⊗···⊗σˆx (cid:17)⊗(cid:16)σˆy ⊗···⊗σˆy (cid:17) 2N χ(1) χ(nx) χ(nx+1) χ(nx+ny) χ (cid:16) (cid:17) (cid:16) (cid:17) ⊗ σˆz ⊗···⊗σˆz ⊗ σˆ0 ⊗···⊗σˆ0 . (C1) χ(nx+ny+1) χ(nx+ny+nz) χ(nx+ny+nz+1) χ(N) Here the parameters n = {n ,n ,n } are positive inte- a valid density matrix, as it may have negative eigenval- x y z gersthatarearbitraryuptotheconstraintn +n +n ≤ ues (equivalently, it need not satisfy Tr[ρˆ2]≤Tr[ρˆ]=1). x y z N, σˆ0 is the identity matrix on the Hilbert space of spin However, the master equation provides a positive map, i i, and χ is a permutation of integers {1,...,N}. The and preserves the positivity of the density matrix eigen- most general permutation symmetric Hermitian matrix values dynamically. Therefore if we start with an initial can be written state that is a valid density matrix, it will remain restricted to the space of valid density matrices during the ρˆ=(cid:88)c(n)Mˆ(n). (C2) time evolution. n The permutation symmetry of the Liouvillian endows itwithablockdiagonalstructure, withoneoftheblocks Note that Mˆ({0,0,0}) has unity trace, while all the spanned by the states ρˆ in Eq. (C2). Dynamics of an other basis states are traceless; hence ρˆ can only be a initial density matrix within this block can therefore be valid density matrix if c({0,0,0}) = 1. The number calculatedbydeterminingtheactiononallstatesMˆ(n): of unconstrained coefficients, and hence the dimensionality of the space that must be considered in the dy- L(Mˆ(n))=(cid:88)L Mˆ(m). (C4) m,n namics, is the number of ways of choosing n such that m 0 < n ,n ,n ≤ N. It is straightforward to check that x y z there are D such n, where The coefficients Lm,n are straightforward to compute, andusingthisEq.(C4)themasterequationcanbewrit- (N +3)(N +2)(N +1) ten D = −1, (C3) 3! ∂ c(n)=L c(m). (C5) t n,m hencethedimensionoftheaccessiblespaceofdensityma- trices is O(N3). Even after the restriction c({0,0,0})= ThissetofO(N3)first-order,linear,ordinarydifferential 1,amatrixoftheforminEq.(C2)isnotguaranteedtobe equations is straightforward to integrate for N (cid:46)100. [1] M. E. Fisher, Rev. Mod. Phys. 46, 597 (1974). [7] A. Amo, J. Lefrére, S. Pigeon, C. Adrados, C. Ciuti, [2] M. E. Fisher, Rev. Mod. Phys. 70, 653 (1998). I.Carusotto,R.Houdré,E.Giacobino, andA.Bramati, [3] S. Sachdev, Quantum phase transitions (Wiley Online Nat. Phys. 5, 805 (2009). Library, 2007). [8] H.Deng,H.Haug, andY.Yamamoto,Rev.Mod.Phys. [4] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, 82, 1489 (2010). P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. [9] A.Polkovnikov,K.Sengupta,A.Silva, andM.Vengalat- Szyman´ska,R.André,J.L.Staehli,V.Savona,P.B.Lit- tore, Rev. Mod. Phys. 83, 863 (2011). tlewood, B. Deveaud, and L. S. Dang, Nature 443, 409 [10] I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2006). (2013). [5] S. Hofferberth, I. Lesanovsky, B. Fischer, T. Schumm, [11] M. O. Scully and W. E. Lamb, Phys. Rev. 159, 208 and J. Schmiedmayer, Nature 449, 324 (2007). (1967). [6] A. Amo, D. Sanvitto, F. P. Laussy, D. Ballarini, E. del [12] V. DeGiorgio and M. O. Scully, Phys. Rev. A 2, 1170 Valle, M. D. Martin, A. Lemaˆıtre, J. Bloch, D. N. (1970). Krizhanovskii, M. S. Skolnick, C. Tejedor, and L. Vina, [13] E. G. D. Torre, S. Diehl, M. D. Lukin, S. Sachdev, and Nature 457, 291 (2008). P. Strack, Phys. Rev. A 87, 023831 (2013). 9 [14] P. Kirton and J. Keeling, Phys. Rev. Lett. 111, 100404 [40] D.G.Angelakis, M.F.Santos, andS.Bose,Phys.Rev. (2013). A 76, 031805(R) (2007). [15] M. Hafezi, P. Adhikari, and J. M. Taylor, Phys. Rev. B [41] M. J. Hartmann and M. B. Plenio, Phys. Rev. Lett. 99, 92, 174305 (2015). 103601 (2007). [16] M. F. Maghrebi and A. V. Gorshkov, Phys. Rev. B 93, [42] D. Rossini and R. Fazio, Phys. Rev. Lett. 99, 186401 014307 (2016). (2007). [17] S. M. Spillane, G. S. Pati, K. Salit, M. Hall, P. Kumar, [43] J.KochandK.L.Hur,Phys.Rev.A80,023811(2009). R. G. Beausoleil, and M. S. Shahriar, Phys. Rev. Lett. [44] S.Diehl,A.Tomadin,A.Micheli,R.Fazio, andP.Zoller, 100, 233602 (2008). Phys. Rev. Lett. 105, 015702 (2010). [18] E. Vetsch, D. Reitz, G. Sagué, R. Schmidt, S. T. [45] A.-W. de Leeuw, O. Onishchenko, R. A. Duine, , and Dawkins, andA.Rauschenbeutel,Phys.Rev.Lett.104, H. T. C. Stoof, Phys. Rev. A 91, 033609 (2015). 203603 (2010). [46] D.S.RokhsarandB.G.Kotliar,Phys.Rev.B44,10328 [19] K. Hennessy, A. Badolato, M. Winger, D. Gerace, (1991). M. Atatu¨re, S. Gulde, S. Fa¨lt, E. L. Hu, and [47] J.Dalibard,Y.Castin, andK.Mølmer,Phys.Rev.Lett. A. Imamo˘glu, Nature 445, 896 (2007). 68, 580 (1992). [20] T. G. Tiecke, J. D. Thompson, N. P. d. Leon, L. R. Liu, [48] R.Dum,P.Zoller, andH.Ritsch,Phys.Rev.A45,4879 V. Vuletic, and M. D. Lukin, Nature 508, 241 (2014). (1992). [21] M. Hafezi, D. E. Chang, V. Gritsev, E. Demler, and [49] M.B.PlenioandP.L.Knight,Rev.Mod.Phys.70,101 M. D. Lukin, Phys. Rev. A 85, 013822 (2012). (1998). [22] J. D. Thompson, T. G. Tiecke, N. P. de Leon, J. Feist, [50] A. J. Daley, Adv. Phys. 63, 77 (2014). A.V.Akimov,M.Gullans,A.S.Zibrov,V.Vuletic, and [51] T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 M. D. Lukin, Science 340, 1202 (2013). (1940). [23] A.V.Gorshkov,J.Otterbach,M.Fleischhauer,T.Pohl, [52] C. Joshi, F. Nissen, and J. Keeling, Phys. Rev. A 88, and M. D. Lukin, Phys. Rev. Lett. 107, 133602 (2011). 063835 (2013). [24] T.Peyronel,O.Firstenberg,Q.-Y.Liang,S.Hofferberth, [53] P.D.DrummondandD.F.Walls,J.ofPhys.A13,725 A. V. Gorshkov, T. Pohl, M. D. Lukin, and V. Vuletic, (1980). Nature (London) 488, 57 (2012). [54] R. Bonfacio and L. A. Lugiato, Phys. Rev. A 18, 1129 [25] A. A. Houck, H. E. Tu¨reci, and J. Koch, Nat. Phys. 8, (1978). 292 (2012). [55] C. M. Bowden and C. C. Sung, Phys. Rev. A 19, 2392 [26] F. Nissen, S. Schmidt, M. Biondi, G. Blatter, H. E. (1979). Tu¨reci, and J. Keeling, Phys. Rev. Lett. 108, 233603 [56] H. J. Carmichael, J. S. Satchell, and S. Sarkar, Phys. (2012). Rev. A 34, 3166 (1986). [27] J.Reftery,D.Sadri,S.Schmidt,H.E.Tu¨reci, andA.A. [57] S. Sarkar and J. S. Satchell, Europhys. Lett. 3, 797 Houck, Phys. Rev. X 4, 031043 (2014). (1987). [28] R.Barends,L.Lamata,J.Kelly,L.Garcia-Alvarez,A.G. [58] B.Nagorny,T.Elsa¨sser, andA.Hemmerich,Phys.Rev. Fowler, A. Megrant, E. Jeffrey, T. C. White, D. Sank, Lett. 91, 153003 (2003). J.Y.Mutus,B.Campbell,Y.Chen,Z.Chen,B.Chiaro, [59] T.Elsa¨sser,B.Nagorny, andA.Hemmerich,Phys.Rev. A. Dunsworth, I. C. Hoi, C. Neill, P. J. J. O’Malley, A 69, 033403 (2004). C. Quintana, P. Roushan, A. Vainsencher, J. Wenner, [60] M.A.ArmenandH.Mabuchi,Phys.Rev.A73,063801 E.Solano, andJ.M.Martinis,Nat.Comm.6,1(2015). (2006). [29] D. E. Chang, V. Gritsev, G. Morigi, V. Vuletić, M. D. [61] T.E.Lee,H.Ha¨ffner, andM.C.Cross,Phys.Rev.Lett. Lukin, and E. A. Demler, Nat. Phys. 4, 884 (2008). 108, 023602 (2012). [30] D. E. Chang, V. Vuletić, and M. D. Lukin, Nat. Phot. [62] M. Marcuzzi, E. Levi, S. Diehl, J. P. Garrahan, and 8, 685 (2014). I. Lesanovsky, Phys. Rev. Lett. 113, 210401 (2014). [31] M.J.Hartmann,F.G.S.L.Brandao, andM.B.Plenio, [63] C.Carr,R.Ritter,C.G.Wade,C.S.Adams, andK.J. Nat. Phys. 2, 849 (2006). Weatherill, Phys. Rev. Lett. 111, 113901 (2013). [32] M.J.Hartmann,F.G.S.L.Brandao, andM.B.Plenio, [64] T. E. Lee, H. Ha¨ffner, and M. C. Cross, Phys. Rev. A Las. & Photon. Rev. 2, 527 (2008). 84, 031402(R) (2011). [33] I. Carusotto, D. Gerace, H. E. Tureci, S. D. Liberato, [65] A. Hu, T. E. Lee, and C. W. Clark, Phys. Rev. A 88, C.Ciuti, andA.Imamog˘lu,Phys.Rev.Lett.103,033601 053627 (2013). (2009). [66] M. Hoening, W. Abdussalam, M. Fleischhauer, and [34] A. Le Boité, G. Orso, and C. Ciuti, Phys. Rev. Lett. T. Pohl, Phys. Rev. A 90, 021603(R) (2014). 110, 233601 (2013). [67] M. J. Hartmann, Phys. Rev. Lett. 104, 113601 (2010). [35] A. Le Boité, G. Orso, and C. Ciuti, Phys. Rev. A 90, [68] J. Jin, D. Rossini, R. Fazio, M. Leib, and M. J. Hart- 063821 (2014). mann, Phys. Rev. Lett. 110, 163605 (2013). [36] H.-P. Breuer and F. Petruccione, The theory of open [69] J. Jin, D. Rossini, M. Leib, and M. J. H. and. R. Fazio, quantum systems (Oxford University Press, 2002). Phys. Rev. A 90, 023827 (2014). [37] H. J. Carmichael, An Open Systems Approach to Quan- [70] G.NicolisandI.Prigogine,Self-OrganizationinNonequi- tum Optics (Springer, Berlin, 1993). librium Systems (Wiley, 1977). [38] M.P.A.Fisher,P.B.Weichman,G.Grinstein, andD.S. [71] J. Qian, G. Dong, L. Zhou, and W. Zhang, Phys. Rev. Fisher, Phys. Rev. B 40, 546 (1989). A 85, 065401 (2012). [39] A. D. Greentree, C. Tahan, J. H. Cole, and L. C. L. [72] C.-K. Chan, T. E. Lee, and S. Gopalakrishnan, Phys. Hollenberg, Nat. Phys. 2, 856 (2006). Rev. A 91, 051601(R) (2015). 10 [73] J.-J. Mendoza-Arenas, S. R. Clark, S. Felicetti, [79] Z. Cai and T. Barthel, Phys. Rev. Lett. 111, 150403 G. Romero, E. Solano, D. G. Angelakis, and D. Jaksch, (2013). Phys. Rev. A 93, 023821 (2015). [80] P.KinslerandP.D.Drummond,Phys.Rev.A43,6194 [74] H. Spohn, Lett. Math. Phys. 2, 33 (1977). (1991). [75] S. G. Schirmer and X. Wang, Phys. Rev. A 81, 062306 [81] H. Weimer, Phys. Rev. Lett. 114, 040402 (2015). (2010). [82] D. Porras and J. I. Cirac, Phys. Rev. Lett. 92, 207901 [76] J. Kerckhoff, M. A. Armen, and H. Mabuchi, Optics (2004). Express 19, 24468 (2011). [83] K. Kim, M.-S. Chang, S. Korenblit, R. Islam, E. E. Ed- [77] A. Dombi, A. Vukics, and P. Domokos, J. Phys. B 46, wards, J. K. Freericks, G.-D. Lin, L.-M. Duan, and 224010 (2013). C. Monroe, Nature 465, 590 (2010). [78] F. Verstraete, M. M. Wolf, and J. I. Cirac, Nat. Phys. [84] M. Xu, D. A. Tieri, and M. J. Holland, Phys. Rev. A 5, 633 (2009). 87, 062101 (2013). [85] C. M. Savage and H. J. Carmichael, IEEE J. Quantum Electron. 24, 1495 (1988).