Breaking strong symmetries in dissipative quantum systems: the (non-)interacting bosonic chain coupled to a cavity
BBreaking strong symmetries in dissipative quantum systems: the (non-)interactingbosonic chain coupled to a cavity
Catalin-Mihai Halati, Ameneh Sheikhan, and Corinna Kollath Physikalisches Institut, University of Bonn, Nussallee 12, 53115 Bonn, Germany (Dated: February 5, 2021)In dissipative quantum systems, strong symmetries can lead to the existence of conservation lawsand multiple steady states. The investigation of such strong symmetries and their consequenceson the dynamics of the dissipative systems is still in its infancy. In this work we investigate astrong symmetry for bosonic atoms coupled to an optical cavity, an experimentally relevant system,using adiabatic elimination techniques and numerically exact matrix product state methods. Weshow the existence of multiple steady states for ideal bosons coupled to the cavity. We find that theintroduction of a weak breaking of the strong symmetry by a small interaction term leads to a directtransition from multiple steady states to a unique steady state. We point out the phenomenon ofdissipative freezing, the breaking of the conservation law at the level of individual realizations inthe presence of the strong symmetry. For a weak breaking of the strong symmetry we see that thebehavior of the individual trajectories still shows some signs of this dissipative freezing before itfades out for a larger symmetry breaking terms.
Symmetries play a key role in order to classify andunify the physics occurring in different microscopic sys-tems. A famous example is the universal behavior arisingat (quantum) phase transitions. This universal behavioris independent of the microscopic details of the systemand can be classified by the symmetries which are spon-taneously broken at the transition. Typically, each sym-metry of a Hamiltonian is connected to a conservationlaw. This has the crucial consequence that the long timestate remembers the initial conditions and that the con-servation laws need to be considered constructing thermalensembles [1]. The ensemble taking these conservationlaws into account has been called the generalized Gibbsensemble [2].Surprisingly, in contrast to the Hamiltonian case, forsystems described by the dissipative Lindblad masterequations ∂∂t ρ = L ( ρ ), where ρ is the density matrixand L the Liouvillian, a symmetry of the Liouvilliandoes not always imply a conserved quantity and mul-tiple steady states [3, 4]. Let the Hermitian operator O be the generator of the symmetry U = exp( iφ O ),with real φ . If the symmetry operator satisfies the con-dition L (cid:0) U ρ U † (cid:1) = UL ( ρ ) U † , we have only a so-calledweak symmetry. In contrast to the Hamiltonian coun-terpart, this weak symmetry condition is not sufficientto imply the existence of a conserved quantity or multi-ple steady states. Only if additionally the generator ofthe symmetry O is commuting with both the Hamilto-nian and all jump operators J m , [ O , H ] = [ O , J m ] = 0, aso-called strong symmetry exists which implies conserva-tion laws. In the presence of a symmetry the Liouvilliancan be block diagonalized in different symmetry sectors.For a strong symmetry, in each symmetry sector at leastone steady state exists, thus giving a lower bound onthe number of steady states. In this case the quantity (cid:104)O(cid:105) = tr( ρO ) is conserved, ∂∂t (cid:104)O(cid:105) = 0.Recently, the consequences of the weak and strongsymmetries in open systems were discussed in the con-text of error correction for quantum information theory [5]. Experimental systems that can be described by aLindblad master equation are very frequent in the areaof quantum optics and solid state systems coupled tolight. In many situations, a unique steady state arises.However, in recent years a significant amount of workhas been devoted to go beyond this typical situation andto study the coexistence of several states in such Lind-blad systems [6–18]. For example, the phenomena of bi-/metastability have been investigated in various systems[15, 19–28]. Additionally, experiments [13–15, 25, 29]tried to identify the coexistence of phases by the phe-nomenon of intermittency, the random switching betweentwo distinct experimental measurement outputs in a sin-gle experimental run. However, care has to be taken con-cerning the interpretation of the occurrence or absence ofintermittency. A system with a unique steady state canshow intermittency if the measurement targets differentdiagonal components as e.g. shown for the stochasticunravelling [30], or intermittency can be absent if thetimescale for the transition between two states diverges,or in the case of disconnected multiple steady states. Thelater can occur due to the presence of a strong symmetry[31]. Additionally, steady states with sought-after prop-erties have been constructed employing symmetries of thesystem, such as steady states with η -pairing correlations[32] or state with enhanced currents [33, 34] or in weaklydriven systems [35].We show in this work for a realistic experimental sys-tem, a cavity coupled to a quantum gas, how the presenceof a strong symmetry can lead to the occurrence of mul-tiple dissipative phase transitions in different symmetrysectors.We further investigate how in the situation of the weakbreaking of the strong symmetry by an additional termin the Liouvillian the unique steady state is recovered.Thus, the weak breaking causes a drastic response of thesystem. Beside identifying the different possible phasetransitions and the corresponding steady states, we inves-tigate the timescales associated to the process of reaching a r X i v : . [ c ond - m a t . qu a n t - g a s ] F e b the steady state. Additionally, we show the absence of in-termittency, the dissipative freezing [31], in the presenceof the strong symmetry in single trajectories obtained bythe stochastic unravelling of the master equation. Wefind that this behavior of the absence of intermittencycan approximately survive for an intermediate time whenadding a small symmetry breaking term. Whereas forlarger symmetry breaking term the different symmetrysectors are no longer a good description of the system.We consider ultracold bosons confined to a one-dimensional chain coupled to a single cavity mode andtransversely pumped with a standing-wave laser beam[30]. The dynamics of the coupled cavity-atom system isdescribed by the Lindblad equation in which the excitedinternal state of the atoms is adiabatically eliminated[30, 36–39] ∂∂t ρ = L ( ρ ) = − i (cid:126) [ H, ρ ] + Γ2 (cid:0) aρa † − a † aρ − ρa † a (cid:1) , (1)where L ( ρ ) is the Liouvillian. The bosonic operators a and a † are the annihilation and creation operators for thephoton mode of the cavity and Γ the dissipation strength.The dissipator represents the losses from the cavity dueto the imperfections of the mirrors. The Hamiltonian, H = H + H int , is given by [38–40] H = (cid:126) δa † a − (cid:126) Ω( a + a † ) (cid:88) k b † k b k + π (mod 2 π ) (2) − J (cid:88) k cos( k ) b † k b k H int = U (cid:88) j n j ( n j − . The cavity mode is described by the first term in H , inthe rotating frame of the pump beam, where δ = ω c − ω p is the detuning between the cavity mode and the trans-verse pump beam. The operators b k and b † k are thebosonic annihilation and creation operators of the atomswith the unitless momentum k j = πjL and j = 1 , . . . , L ,assuming periodic boundary conditions. We note that inthe numerical results we considered open boundary con-ditions [41]. L denotes the number of sites of the bosonicchain and the total number of bosons is N . The cou-pling between the atoms and the cavity field is describedby the second term and introduces a coupling betweenthe momentum k and k + π (mod 2 π ). This is due tothe periodicity of cavity mode which has twice the peri-odicity of the lattice spacing within the chain. The lastterm describes the tunneling processes of the atoms alongthe chain with the tunneling amplitude J . The repulsiveon-site interactions of strength U ≥ H int , where j denotes the site ofthe chain and n j the atomic density.For U = 0, only transitions between the occupationof the momenta k and k + π (mod 2 π ) of the atoms arepossible. In the single particle case, L/ FIG. 1. The scaled photon number, (cid:104) a † a (cid:105) /N , as a functionof (cid:126) Ω √ N/J , for different symmetry sectors. We compare thevalues obtained with tMPS at tJ/ (cid:126) = 200 with the general-ized Gibbs ensemble within the mean field approach includ-ing thermal fluctuations. The parameters used are L = 10, N = 5, (cid:126) δ/J = 2, U/J = 0, and (cid:126) Γ /J = 1. At the dashedvertical line, the nature of the obtained steady states reachesfrom states with an empty cavity (red symbol/line) to theself-organized states (remaining ones). symmetry sectors exist, each spanned by the momentumstates | k j (cid:105) and | k j + π (mod 2 π ) (cid:105) , j = 1 , . . . , L/
2. Theseare corresponding to the strong symmetry, having asgenerator the atomic number operators in each symme-try sector, O k j = b † k j b k j + b † k j + π (mod 2 π ) b k j + π (mod 2 π ) ,and their average values are conserved quantities, m k j = (cid:104)O k j (cid:105) . Due to the strong symmetry, already for a singleparticle multiple steady states exist.For N atoms, the symmetry sectors can be con-structed from the different combinations in which onecan arrange the atoms in the single particle sectors.Thus, each symmetry sector will be labeled by K ≡ (cid:0) m k , ..., m k i , ..., m k L/ (cid:1) , with (cid:80) L/ i =1 m k i = N . However,even though the atoms can be arranged independentlyin the single particle sectors, they are coupled via thephoton field.We show in Fig. 1 the steady state diagram of cou-pled atom cavity system using two different methods: (i)a mean field decoupling of the atomic and the photonicsector considering the fluctuations as a perturbation to-gether with the assumption that within each symmetrysector the atoms thermalize [41, 42] and (ii) a matrixproduct state (MPS) method developed [43] for the nu-merically exact simulation of the time-evolution of thedissipative master equation, Eqs. (1)-(2) (the numericalparameters are given in [41]). For the time-evolution wehave chosen the empty cavity and the ground state in theatomic sector as the initial state.Both methods give in each considered symmetry sec-tor of the strong symmetry a transition from the emptycavity state to the self-organized state with finite cavityoccupation. Importantly, the transitions in the differentsymmetry sectors take place at distinct critical values ofthe coupling strength Ω c . Physically, this can be under-stood in the following way: the self-organization transi-tion arises due to a competition between the ordering ofthe atoms in a density wave induced by interaction withthe photon mode and the kinetic energy of the atomswhich depends on the momentum of the atoms. As ineach symmetry sector the atoms have different momenta,this gives rise to different critical values for the transi-tion. The first approach (i) is a perturbative approachusing fluctuations around the mean-field solution [42]. Inthe presence of the strong symmetry we adapt this ap-proach by using the thermalization of the atoms withineach ( k , k + π (mod 2 π ))-sector with sector-dependenttemperatures. This corresponds to a generalization ofthe ’generalized’ Gibbs ensemble to dissipative systems[41].Thus, at a fixed coupling strength (cf. vertical line inFig. 1) multiple steady states arise depending on the pro-jection of the initial state to the symmetry sectors. Forthe considered strong symmetry, these steady states caneven have very different nature. It can occur that onesector is still in the disordered phase with an empty cav-ity, whereas another sector already is deep in the self-organized phase. This is to be contrasted to the meta-stable states arising for weak symmetries, which are typ-ically connected to a unique steady state.In the limit of large dissipation strength, we can alsodetermine the multiple steady states using the kineticterm of the atoms as a perturbation [43]. In the presenceof interaction, we had shown [30, 43] that the steady statecorresponds to a state with a totally mixed atomic sec-tor, ρ mix = N (cid:80) { n i } | α (∆) , ∆ (cid:105) (cid:104) α (∆) , ∆ | , with the cav-ity field α (∆) = Ω δ − i Γ / ∆, and the odd-even imbalance∆ = (cid:80) L/ j =1 (cid:104)O k j (cid:105) = (cid:80) Lj =1 ( − j (cid:104) n j (cid:105) . Here we have gen-eratlized this approach taking the different conservationlaws into account [41]. This leads to a totally mixed statewithin each symmetry sector which can be very differ-ent in nature compared to the steady state ρ mix in thepresence of a small perturbation breaking the symmetry.Thus, also in the strong dissipation limit we find that theweak breaking of the strong symmetry can lead to drasticeffects in the nature of the steady state. Dissipative freezing . Recently, the effect of the exis-tence of a strong symmetry in Liouvillian on the timeevolution of the quantum trajectories was analyzed andthe phenomenon of ’dissipative freezing’ was shown forsystems in which the Lindblad operator is proportionalto the Hamiltonian ( H ∝ L ) [31]. Dissipative freezingis the phenomenon that single realizations of trajecto-ries, obtained by the stochastic unravelling of the masterequation, can break the strong symmetry. A trajectorywhich is purely in one symmetry sector, will remain inthat sector for the rest of the time evolution and thus,obey the symmetry of the system. However, startingwith an initial state which is a superposition with contri-butions from multiple symmetry sectors, each individualtrajectory will randomly select one of the sectors and FIG. 2. Time evolution of O k for the single quan-tum trajectories sampled in the Monte Carlo averagefor different interaction strengths U , with k = k and k = k . The initial state consists in an equal su-perposition between states for the sectors ( m k = 5) and( m k = 1 , m k = 1 , m k = 1 , m k = 1 , m k = 1). In eachpanel there are 1000 trajectories plotted, the black representthe Monte Carlo average, either for the full set of trajectories,or averaged separately depending on the final value, we shadethe interval of one standard deviation away from the aver-age, with light blue for the full average and light gray for theseparate averages. The parameters used are L = 10, N = 5, (cid:126) δ/J = 2, (cid:126) Ω √ N/J = 4 .
47, and (cid:126) Γ /J = 15. remain there for the rest of the evolution. Further, nointermittency occurs between these trajectories in differ-ent sectors. Thus, even though the Monte Carlo averageexpectation value of the generator of the symmetry isa conserved quantity, this is no longer true at the levelof single trajectories. The single trajectories can break’spontaneously’ the strong symmetry of the model. Thisis an intriguing effect which might have relevance to thesingle realizations of experiments. Further, the interpre-tation of quantities measured in the quantum trajectorymethod which stabilize at different values, should notimply different steady states, but reflect the initial su-perposition.Here we numerically show [Fig. 2(a)-(b)] that even ina system that goes beyond the special case ( H ∝ L ) ofRef. [31], dissipative freezing can occur. In Fig. 2(a)-(b)we plot the expectation value of two of the generators ofthe strong symmetry, (cid:104)O k (cid:105) and (cid:104)O k (cid:105) , in time for 1000single trajectories [41]. The initial state is an equal su-perposition of a state from the sector ( m k = 5) and thesector ( m k = 1 , m k = 1 , m k = 1 , m k = 1 , m k = 1).We can observe that at long times, tJ (cid:38) (cid:126) , all tra-jectories evolved to one of the two symmetry sectors, as (cid:104)O k (cid:105) ( t ) equals the expected occupation in those sectors.The Monte Carlo average of the trajectories stays con-stant throughout the following time-evolution, up to anumerical error. Breaking of the strong symmetry . For any finite inter-action,
U >
0, the operators O k , no longer commute withthe Hamiltonian, thus the strong symmetry of the Liou-villian is broken. We analyze how the system passes overfrom having multiple steady states to a unique steadystate as the on-site interaction is slowly turned on. We fo-cus on the limit of large dissipation, for which the many-body adiabatic elimination predicts the steady state tran-sition between the multiple steady states at U = 0, ρ K, st [41] to a single steady state which is the totally mixedstate, ρ mix [43].In Fig. 3(a)-(b) the behavior of the expectation valueof the photon number and the conserved quantities ofthe symmetry at fixed time tJ = 49 . (cid:126) is plotted as afunction of the interaction strength. We consider initialstates in different symmetry sectors. For these, we ob-serve that at U = 0 multiple steady states are obtainedsignaled by distinct expectation values. However, as theinteraction strength is increased the values of the pho-ton number and (cid:104)O k (cid:105) for the different initial states startto be more and more similar until they agree with eachother and with the values expected for ρ mix , for large val-ues of interaction U . The deviations from the predictedunique steady state ρ mix for small interaction strengthare due to the fact that the system has not yet reachedits steady state at the shown time. This can be observedin time-evolution plots given in Fig. 3(c)-(d) for the pho-ton number and (cid:104)O k (cid:105) − (cid:104)O k (cid:105) ρ mix for different interactionvalues. The expected steady state value for (cid:104) a † a (cid:105) /N isrepresented with a gray line in Fig. 3(c)and (cid:104)O k (cid:105) ρ mix = 1.In order to quantify this we fitted the time dependence of (cid:104) a † a (cid:105) − (cid:104) a † a (cid:105) ρ mix and (cid:104)O k (cid:105) − (cid:104)O k (cid:105) ρ mix with an exponen-tial function, ∝ e − t/τ , and extracted the timescales forreaching the steady state. The fits describe very nicelythe numerical data, which gives strong support that at in-finite time the steady state is given by ρ mix . Additionally,the dependence of the timescale on U is represented inFig. 3(e)-(f) in log-log plots. The timescales exhibits analgebraic dependence on 1 /U . We compare this behav-ior with results from the exact diagonalization of eitherthe full Liouvillian, Eq. (1), (ED) or the many-body adi-abatic elimination equations of motion, with the kineticenergy as a perturbation [43], (AE), for a small system of L = 4. In these case we compute the timescale as the in-verse of the real part of the first excited eigenvalue and weobtain an algebraic dependence ∝ U − α with α ≈ (cid:104)O k (cid:105) are consistently larger than the timescales forthe photon number, which signals that in our simulationthe photon state is reaching the steady state before theatomic one. We attribute this to the spatial extend ofthe atomic system. FIG. 3. The dependence on the interaction strength U of(a) the scaled photon number, (cid:104) a † a (cid:105) /N , and (b) the expec-tation value of O k using tMPS at time tJ = 49 . (cid:126) andmany-body adiabatic elimination (AE). The time evolutionof (c) the scaled photon number, (cid:104) a † a (cid:105) /N , and (d) the ex-pectation value of O k for different values of U . For finite U we fit the time evolution with an exponential decay (blackdashed lines) the difference between the tMPS data and theexpected steady state value, obtained from many-body adia-batic elimination, using the kinetic energy as a perturbation.The timescales obtained from the exponential fits as a func-tion of U for (e) (cid:104) a † a (cid:105) /N and (f) O k j . In (f) we comparethe timescales of O k j with the longest timescale of a smallsystem of N = 2 particles in L = 4 system computed withexact diagonalization (ED) and many-body adiabatic elimina-tion (AE). We fit the timescale dependence on the interactionwith an algebraic decay ∝ U − α and obtain the following ex-ponents: (e) ( m k = 5), α = 1 . ± .
07; ( m k = 1 , m k = 4), α = 1 . ± .
08; (f) j = 2, α = 1 . ± . j = 3, α = 1 . ± . j = 4, α = 1 . ± .
03. The red linesare a guide to the eye of an algebraic decay ∝ U − . The pa-rameters are chosen to be L = 10, N = 5, (cid:126) Ω √ N/J = 4 . (cid:126) δ/J = 2, and Γ /J = 15. Thus, we see that the time-evolution at short timesremembers well the strong symmetry and the mixing ofthe different symmetry sectors only occurs on timescales ∝ /U associated with the scattering of the atoms.A question which arises is how the breaking of thesymmetry affects the phenomenon of dissipative freez-ing discussed before. We study in the following how thisphenomenon is affected by the presence of a weak inter-action, the symmetry breaking term.We observe in Figs. 2(c)-(d) that at U/J = 0 .
01 thetime evolution found for the single trajectories resemblesat early time the one at
U/J = 0. This means that thesingle trajectories break the approximate strong symme-try and approach the two different symmetry sectors. Atintermediate time, 20 (cid:46) tJ/ (cid:126) , many of the quantum tra-jectories spend a long time near the two values expectedfrom the U = 0 symmetry sectors. Only few of the trajec-tories directly show deviations from these values or inter-mittency between the values such that the phenomenonof ’dissipative freezing’ also occurs here to an approxi-mate extend. One has to be careful not to misinterpretthis absence of intermittency as the existence of multiplesteady state. Only starting from U/J (cid:38) .
25 the effect ofthe strong symmetry washes out in the considered timeinterval.To summarize, we analyzed the effects of a strong sym-metry and weak breaking of this symmetry on the dy-namics of a many-body open system consisting of bosonicatoms coupled to an optical cavity. The strong symme-try implies the existence of multiple steady states and weshowed that the dissipative phase transition to the self-organized state can occur at different thresholds in dif-ferent symmetry sectors described by generalized Gibbsensembles in the atomic part. We analyzed how the na-ture of the steady state changes drastically when a smallterm that breaks the strong symmetry is introduced. Thetimescales towards the new unique steady state where found to be proportional to 1 /U associated with thescattering between the atoms. We have shown that evenfor a many-body system with a strong symmetry the phe-nomenon of dissipative freezing can occur when one con-siders the behavior of individual quantum trajectories.It appears such that, at intermediate time, one can stillidentify the effect of dissipative freezing even if the strongsymmetry has been slightly broken. An open questionremains, whether a spontaneous symmetry breaking canalso be observed in single trajectories of an experimen-tal measurement. This would question the interpretationof the absence of intermittency in experimental measure-ments. Acknowledgments
We thank J.-S. Bernier, A. V.Bezvershenko, M. Fleischhauer, A. Rosch and S. Wolfffor stimulating discussions. We acknowledge fundingfrom the Deutsche Forschungsgemeinschaft (DFG, Ger-man Research Foundation) in particular under projectnumber 277625399 - TRR 185 (B3) and project num-ber 277146847 - CRC 1238 (C05) and under Germany’sExcellence Strategy – Cluster of Excellence Matter andLight for Quantum Computing (ML4Q) EXC 2004/1 –390534769 and the European Research Council (ERC)under the Horizon 2020 research and innovation pro-gramme, grant agreement No. 648166 (Phonton). [1] R. Balian,
From Microphysics to Macrophysics (Springer-Verlag Berlin Heidelberg, 2007).[2] V. Y. M. Rigol, V. Dunjko and M. Olshanii,Phys. Rev. Lett. , 050405 (2007).[3] B. Buˇca and T. Prosen, New Journal of Physics ,073007 (2012).[4] V. V. Albert and L. Jiang, Phys. Rev. A , 022118(2014).[5] S. Lieu, R. Belyansky, J. T. Young, R. Lundgren, V. V.Albert, and A. V. Gorshkov, Phys. Rev. Lett. ,240405 (2020).[6] E. M. Kessler, G. Giedke, A. Imamoglu, S. F. Yelin, M. D.Lukin, and J. I. Cirac, Phys. Rev. A , 012116 (2012).[7] F. Minganti, A. Biella, N. Bartolo, and C. Ciuti, Phys.Rev. A , 042118 (2018).[8] H. J. Carmichael, Phys. Rev. X , 031028 (2015).[9] H. Weimer, Phys. Rev. Lett. , 040402 (2015).[10] M. Benito, C. S´anchez Mu˜noz, and C. Navarrete-Benlloch, Phys. Rev. A , 023846 (2016).[11] L. M. Sieberer, S. D. Huber, E. Altman, and S. Diehl,Phys. Rev. Lett. , 195301 (2013).[12] C. S´anchez Mu˜noz, A. Lara, J. Puebla, and F. Nori, Phys.Rev. Lett. , 123604 (2018).[13] M. Biondi, G. Blatter, H. E. T¨ureci, and S. Schmidt,Phys. Rev. A , 043809 (2017).[14] M.-J. Hwang, P. Rabl, and M. B. Plenio, Phys. Rev. A , 013825 (2018).[15] J. J. Mendoza-Arenas, S. R. Clark, S. Felicetti,G. Romero, E. Solano, D. G. Angelakis, and D. Jaksch,Phys. Rev. A , 023821 (2016).[16] H. Wilming, M. J. Kastoryano, A. H. Werner, and J. Eis-ert, Journal of Mathematical Physics , 033302 (2017). [17] J. Hannukainen and J. Larson, Phys. Rev. A , 042113(2018).[18] J. a. S. Ferreira and P. Ribeiro, Phys. Rev. B , 184422(2019).[19] K. Macieszczak, M. Gut¸˘a, I. Lesanovsky, and J. P. Gar-rahan, Phys. Rev. Lett. , 240404 (2016).[20] T. Fink, A. Schade, S. H¨ofling, C. Schneider, andA. Imamoglu, Nature Physics , 365 (2018).[21] C. Carr, R. Ritter, C. G. Wade, C. S. Adams, and K. J.Weatherill, Phys. Rev. Lett. , 113901 (2013).[22] N. R. de Melo, C. G. Wade, N. ˇSibali´c, J. M. Kondo, C. S.Adams, and K. J. Weatherill, Phys. Rev. A , 063863(2016).[23] S. R. K. Rodriguez, W. Casteels, F. Storme, N. Car-lon Zambon, I. Sagnes, L. Le Gratiet, E. Galopin,A. Lemaˆıtre, A. Amo, C. Ciuti, and J. Bloch, Phys. Rev.Lett. , 247402 (2017).[24] F. Letscher, O. Thomas, T. Niederpr¨um, M. Fleis-chhauer, and H. Ott, Phys. Rev. X , 021020 (2017).[25] P. R. Muppalla, O. Gargiulo, S. I. Mirzaei, B. P.Venkatesh, M. L. Juan, L. Gr¨unhaupt, I. M. Pop, andG. Kirchmair, Phys. Rev. B , 024518 (2018).[26] M. J. A. Schuetz, E. M. Kessler, L. M. K. Vandersypen,J. I. Cirac, and G. Giedke, Phys. Rev. Lett. , 246802(2013).[27] M. J. A. Schuetz, E. M. Kessler, L. M. K. Vandersypen,J. I. Cirac, and G. Giedke, Phys. Rev. B , 195310(2014).[28] L. Hruby, N. Dogra, M. Landini, T. Donner, andT. Esslinger, Proceedings of the National Academy ofSciences , 3279 (2018).[29] M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch, and A. A. Houck, Phys. Rev. X , 011016 (2017).[30] C.-M. Halati, A. Sheikhan, H. Ritsch, and C. Kollath,Phys. Rev. Lett. , 093604 (2020).[31] C. S´anchez Mu˜noz, B. Buˇca, J. Tindall, A. Gonz´alez-Tudela, D. Jaksch, and D. Porras, Phys. Rev. A ,042113 (2019).[32] J.-S. Bernier, P. Barmettler, D. Poletti, and C. Kollath,Phys. Rev. A , 063608 (2013).[33] D. Manzano and P. I. Hurtado, Phys. Rev. B , 125138(2014).[34] F. Lange, Z. Lenarˇciˇc, and A. Rosch, Nature Communi-cations , 15767 (2017).[35] Z. Lenarˇciˇc, F. Lange, and A. Rosch, Phys. Rev. B ,024302 (2018) .[36] H. Carmichael, An open systems approach to quantumoptics (Springer Verlag, Berlin Heidelberg, 1991).[37] H. P. Breuer and F. Petruccione,
The theory of openquantum systems (Oxford University Press, Oxford,2002).[38] H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger,Rev. Mod. Phys. , 553 (2013).[39] C. Maschler, I. B. Mekhov, and H. Ritsch, The EuropeanPhysical Journal D , 545 (2008).[40] G. S. D. Nagy and P. Domokos, Eur. Phys. J. D , 127(2008).[41] See Supplemental material.[42] A. V. Bezvershenko, C.-M. Halati, A. Sheikhan, C. Kol-lath, and A. Rosch, (2020), arXiv:2012.11823.[43] C.-M. Halati, A. Sheikhan, and C. Kollath, Phys. Rev.Research , 043255 (2020).[44] J. J. Garc´ıa-Ripoll, S. D¨urr, N. Syassen, D. M. Bauer,M. Lettner, G. Rempe, and J. I. Cirac, New Journal ofPhysics , 013053 (2009).[45] F. Reiter and A. S. Sørensen, Phys. Rev. A , 032111(2012).[46] D. Poletti, P. Barmettler, A. Georges, and C. Kollath,Phys. Rev. Lett. , 195301 (2013).[47] B. Sciolla, D. Poletti, and C. Kollath, Phys. Rev. Lett. , 170401 (2015).[48] A. C. Y. Li, F. Petruccione, and J. Koch, Scientific Re-ports , 4887 (2014).[49] C.-M. Halati, in preparation , Ph.D. thesis, University ofBonn (2021).[50] J. Dalibard, Y. Castin, and K. Mølmer, Phys. Rev. Lett. , 580 (1992).[51] C. W. Gardiner, A. S. Parkins, and P. Zoller, Phys. Rev.A , 4363 (1992).[52] A. J. Daley, Advances in Physics , 77 (2014).[53] S. R. White and A. E. Feiguin, Phys. Rev. Lett. ,076401 (2004).[54] A. J. Daley, C. Kollath, U. Schollw¨ock, and G. Vidal,J. Stat. Mech.: Theor. Exp. P04005 (2004).[55] U. Schollw¨ock, Annals of Physics , 96 (2011),.[56] E. M. Stoudenmire and S. R. White, New Journal ofPhysics , 055026 (2010).[57] M. L. Wall, A. Safavi-Naini, and A. M. Rey, Phys. Rev.A , 053637 (2016). SUPPLEMENTAL MATERIALA. OPEN BOUNDARY CONDITIONS
In the main text the Liouvillian describing the system,Eqs.(1)-(2), was given in momentum space, in real spacethe Hamiltonian of the system reads H = H c + H kin + H ac , (A.1) H c = (cid:126) δa † a,H kin = − J L − (cid:88) j =1 ( b † j b j +1 + b † j +1 b j ) ,H ac = − (cid:126) Ω( a + a † ) L (cid:88) j =1 ( − j n j ,H = H + H int ,H int = U L (cid:88) j =1 n j ( n j − . In our numerical simulations for finite size systems weused open boundary conditions. In this case one uses theFourier sine transform, defined by b k = (cid:114) L + 1 L (cid:88) j =1 b j sin( kj ) , (A.2) b † k = (cid:114) L + 1 L (cid:88) j =1 b † j sin( kj ) , and the unitless momenta are given by k = πmL +1 and m =1 , . . . , L . In this case still L/ | k j (cid:105) and | π − k j (cid:105) , j = 1 , . . . , L/ k j are in the interval [0 , π ], the values π − k j will always be inside the first Brillouin zone. Thesymmetry generator are now given by O k j = b † k j b k j + b † π − k j b π − k j . (A.3)We note that in the calculations presented in Sec. B weemploy open boundary conditions. B. MANY BODY ADIABATIC ELIMINATIONFORMALISM
We employ a variant of the many-body adiabatic elim-ination method [44–47], which provides analytical insightinto the long-time behavior of our system. This approachis a perturbative approach which considers that the ef-fect of one of the terms in the Hamiltonian, H ν , or ofthe fluctuations around the mean-field solution in thedynamics of the system is weak. We used this procedure in Ref. [30, 42, 43] to determine the phase diagram ofsteady states in the presence of interactions.We decompose the Liouvillian as L = L − i (cid:126) [ H ν , · ] intoan unperturbed Linbladian L and a perturbative contri-bution caused by H ν . This approach captures the effec-tive dynamics of the density matrix in the decoherencefree subspace of L , i.e. the space formed by all densitymatrices ρ which are eigenstates of the superoperator L with a vanishing real part of the eigenvalues. Theother subspaces corresponding to the eigenvalues with anon-zero real part are considered via virtual transitionswithin the perturbation theory. The resulting effectivedynamics in the decoherence free subspace is given by[43, 44, 46, 47] ddt ρ = L ρ + 1 (cid:126) P (cid:2) H ν , L − P (cid:2) H ν , ρ (cid:3)(cid:3) , (B.1)where ρ lies in the decoherence free subspace of L and P and P are the projectors onto the decoherence freesubspace and the first excited subspace, respectively. Inthe following we detail how to apply this approach. InSec. B I the perturbation is the fluctuations around themean-field theory and in Sec. B II the perturbation is thekinetic term of the atoms.We note that due to the applied perturbative expan-sion the condition of positive definiteness might not befulfilled for the obtained density matrix [48]. I. Mean field decoupling with thermal fluctuations
Following the approach introduced in Ref. [42] we per-form a mean field decoupling of the term coupling thecavity and the atoms, H ac , and consider the fluctuationsin the coupling as the perturbation in the many body adi-abatic elimination derivation of the effective equations ofmotion. In this situation we have L = − i (cid:126) [ H c + H kin + H MFac , · ] + D ( · ) , (B.2) H ν ≡ δH ac where D ( ρ ) = Γ2 (cid:0) aρa † − a † aρ − ρa † a (cid:1) , and H MFac = − (cid:126) Ω( α + α ∗ ) (cid:88) k b † k b π − k (B.3) − (cid:126) Ω( a + a † )∆ ,δH ac = − (cid:126) Ω (cid:0) a + a † − α − α ∗ (cid:1) (cid:34)(cid:88) k b † k b π − k − ∆ (cid:35) , . Here α is the mean filed value of (cid:104) a (cid:105) which depends onthe mean field value of the imbalance ∆, α (∆) = Ω δ − i Γ / , (B.4)∆ = L/ (cid:88) j =1 (cid:104)O k j (cid:105) = (cid:88) j ( − j (cid:104) n j (cid:105) . Here we see that the atomic and photonic contributionsdecouple beside the self-consistent determination of theparameters. Thus, a general state in the decoherence freesubspace of L is given by ρ = | α (∆) (cid:105)(cid:104) α (∆) | · ρ b , with ρ b = (cid:88) n,m c n,m | n ( α ) (cid:105) (cid:104) m ( α ) | , (B.5)where | n ( α ) (cid:105) are eigenstates of the atomic mean fieldHamiltonian H a = H kin + H MFac and the cavity part isin a coherent state. If we plug in ρ given by Eq. (B.5)in Eq. (B.1) we obtain the equation of motion for theentries c n,m of the density matrix [49]. A substantialsimplification can occur if we consider a generalization ofthe thermal Ansatz introduced in Ref. [42]. As we onlyhave cavity mediated interactions we cannot assume thatthe atomic sector thermalizes as a whole, but we canmake the assumption that the particles in each singleparticle symmetry sector can thermalize and be describedby an effective inverse temperature β j associated withthis sector, i.e. ρ b ∼ L/ (cid:89) j =1 exp [ − β j H a ( k j , π − k j )] , where (B.6) H a ( k j , π − k j ) = − (cid:126) Ω( α + α ∗ ) (cid:16) b † k j b π − k j + b † π − k j b k j (cid:17) − J cos( k j ) (cid:16) b † k j b k j − b † π − k j b π − k j (cid:17) . For example, in Fig. 1, the state from the symmetry sec-tor ( m k = 5) is described by a single temperature asall particles are in the single particle sector with mo-mentum k , but the state from the symmetry sector( m k = 3 , m k = 2) is described by two temperatures aswe have 3 particles in the single particle sector with mo-mentum k and 2 particles in the single particle sectorwith momentum k . This procedure is analogously tothe consideration of different conservation laws in closedsystems, described by generalized Gibbs ensembles. Letus note that a generalized Gibbs ensemble was also in-troduced in Ref. [35] for weakly driven systems in thepresence of approximate conservation laws.As the density matrix is now determined by a smallnumber of parameters, the inverse temperatures β j , it isenough to consider the equations of motion for a reducednumber of observables. Thus, we describe the steadystate of the system with the temperatures for which theAnsatz given by Eqs. (B.5)-(B.6) satisfies the equations (cid:10) ∂∂t H a ( k j , π − k j ) (cid:11) = 0 for all momenta k j and the mean-field self-consistency condition Eq. (B.4). II. Perturbation in kinetic energy
In the following, we consider the perturbation to begiven by the kinetic energy, H ν ≡ H kin , valid in the regime (cid:126) Γ (cid:29) (cid:126) Ω , (cid:126) δ (cid:29) J , thus L is given by L = − i (cid:126) [ H c + H int + H ac , · ] + D ( · ) , (B.7)with D ( ρ ) = Γ2 (cid:0) aρa † − a † aρ − ρa † a (cid:1) the dissipator.We start by constructing the dissipation free subspaceand the excited subspaces of L which will enter our cal-culations, for an arbitrary interaction strength U ≥ ρ = | α (∆); ∆ , u (cid:105) (cid:104) α (∆ (cid:48) ); ∆ (cid:48) , u (cid:48) | (B.8)are right eigenstates of the superoperator L . We notethat we do not assure that these states are physical den-sity matrices. The atomic part is characterized by theeven-odd imbalance ∆ = (cid:80) j ( − j (cid:104) n j (cid:105) and by its totalinteraction energy, u = U (cid:80) j (cid:104) n j ( n j − (cid:105) . Photons arein a coherent state depending on the atomic imbalance α (∆) = Ω δ − i Γ / . (B.9)The eigenvalues corresponding to the eigenvectors fromEq. (B.8) are given by λ = −
12 Ω Γ δ + Γ / − ∆ (cid:48) ) (B.10)+ i (cid:20) Ω δδ + Γ / − ∆ (cid:48) ) − ( u − u (cid:48) ) (cid:21) . We can observe that for ∆ = ∆ (cid:48) the real part of theeigenvalues is zero. Thus, the states in Eq. (B.8) with∆ = ∆ (cid:48) lie in the decoherence free subspace of L .We can include the contributions from the excited sub-spaces that are coupled to the dissipation free subspacevia the perturbation H kin . With a hopping event we cancouple to the subspace spanned by the states in which∆ = ∆ (cid:48) ± N particles in L sites and can show that the mixed stategiven by [43] ρ mix = 1 N (cid:88) { n j } | α (∆); n , . . . , n L (cid:105) (cid:104) α (∆); n , . . . , n L | (B.11)is a steady state of the system. Here N is the numberof ways one can arrange N identical particles in L sites, N = (cid:18) L + N − N (cid:19) .In the rest of this section we deal with the case of non-interacting atoms, L = − i (cid:126) [ H c + H ac , · ] + D ( · ) . (B.12)As we are in the non-interacting case with the momentumlabeling the different symmetry sectors, we first computethe steady state of one particle in a certain symmetrysector and, afterwards, generalize this result to the caseof N particles.
1. The single particle case
For a single particle in the symmetry sector K =( m k = 1), a general state in the dissipation free subspacerestricted to this symmetry sector has the form ρ k ( b ) = (cid:88) i,j odd (cid:18)
12 + b (cid:19) sin( ki ) sin( kj ) |− α ; i (cid:105) (cid:104)− α ; j | (B.13)+ (cid:88) i,j even (cid:18) − b (cid:19) sin( ki ) sin( kj ) | α ; i (cid:105) (cid:104) α ; j | , with b a real parameter, i and j the positions of theparticle, and α = Ω δ − i Γ / the cavity field.The equation of motion for a state ρ = |± α ; i (cid:105) (cid:104)± α ; j | ,with i and j both even or both odd, from the dissipationfree subspace, obtained with the many-body adiabaticelimination is given by [49] ddt |± α ; i (cid:105) (cid:104)± α ; j | = P (cid:2) H kin , L − P [ H kin , |± α ; i (cid:105) (cid:104)± α ; j | ] (cid:3) (B.14)= J λ e − | α | (cid:16) |± α ; i (cid:105) (cid:104)± α ; j | − |∓ α ; i + 1 (cid:105) (cid:104)∓ α ; j + 1 | − |∓ α ; i + 1 (cid:105) (cid:104)∓ α ; j − |− |∓ α ; i − (cid:105) (cid:104)∓ α ; j + 1 | − |∓ α ; i − (cid:105) (cid:104)∓ α ; j − | + |± α ; i + 2 (cid:105) (cid:104)± α ; j | + |± α ; i − (cid:105) (cid:104)± α ; j | + |± α ; i (cid:105) (cid:104)± α ; j + 2 | + |± α ; i (cid:105) (cid:104)± α ; j − | (cid:17) , with λ = − Γ δ +Γ / . From this we can write the equationof motion for the state ρ k ( b ), and for b = 0 we obtain thesteady state for the one particle case ρ k, st = (cid:88) i,j odd sin( ki ) sin( kj ) |− α ; i (cid:105) (cid:104)− α ; j | (B.15)+ (cid:88) i,j even sin( ki ) sin( kj ) | α ; i (cid:105) (cid:104) α ; j | . This state has a fully mixed atomic sector in the momen- tum basis.
2. The two particle case
We consider two particles in the sector K = ( m k =1 , m k = 1). One can determine the basis in the dissi-pation free subspace in the considered symmetry sectorand compute the equations of motion for these states [49].From this one obtains the steady state solutions [49] ρ k ,k , st = (cid:88) i ,i ,j ,j odd sin( k i ) sin( k j ) sin( k i ) sin( k j ) √ n i n j |− α ; i , i (cid:105) (cid:104)− α ; j , j | (B.16)+ (cid:88) i ,i ,j ,j even sin( k i ) sin( k j ) sin( k i ) sin( k j ) √ n i n j | α ; i , i (cid:105) (cid:104) α ; j , j | + (cid:88) i ,j odd i ,j even sin( k i ) sin( k j ) sin( k i ) sin( k j ) √ n i n j | i , i (cid:105) (cid:104) j , j | + (cid:88) i ,j odd i ,j even sin( k i ) sin( k j ) sin( k i ) sin( k j ) √ n i n j | i , i (cid:105) (cid:104) j , j | , i , i , j and j the positions of the two particlesin the ket or bra, and n i the number of particles at site i , and α = δ − i Γ / If we trace out the photon states werecover a fully mixed atomic sectortr photons ρ k ,k , st = 14 ( | k , k (cid:105) (cid:104) k , k | (B.17)+ | k , π − k (cid:105) (cid:104) k , π − k | + | π − k , k (cid:105) (cid:104) π − k , k | + | π − k , π − k (cid:105) (cid:104) π − k , π − k | ) . Thus, similar as in the case of the interacting system, inthe limit of large dissipation and small kinetic energy, theinfinite temperature state of the corresponding symmetry block is reached. N particle case Generalizing our previous finding to the N particlecase, the steady state for N particles will also be thefully mixed state in the different symmetry sectors.For the N particle case the symmetry sectors can beconstructed from different combinations in which one canarrange the particles in the single particle sectors. Thus,the totally mixed state for N particles distributed in the L/ K = (cid:0) m k , ..., m k i , ..., m k L/ (cid:1) ,is given by ρ K,st = 1 N m k (cid:88) i =0 · · · m kL/ (cid:88) i K =0 (cid:12)(cid:12) n k = i , n π − k = m k − i ; · · · ; n k L/ = i L/ , n π − k L/ = m k L/ − i L/ (cid:11) (B.18) (cid:10) n k = i , n π − k = m k − i ; · · · ; n k L/ = i L/ , n π − k L/ = m k L/ − i L/ (cid:12)(cid:12) , with (cid:12)(cid:12) n k = i , n π − k = m k − i ; · · · ; n k L/ = i L/ , n π − k L/ = m k L/ − i L/ (cid:11) ≡≡ M L (cid:88) j ,...,j mk =0 · · · L (cid:88) j L/ ,...,j L/ mkL/ =0 (cid:16) sin (cid:0) k j (cid:1) .. sin (cid:16) k j m k (cid:17) ( − j i + ... + j mk +( m k − i ) (cid:17) × · · ·· · · × (cid:16) sin (cid:16) k L/ j L/ (cid:17) .. sin (cid:16) k L/ j L/ m kL/ (cid:17) ( − j iL/ + ... + j mkL/ +( m kL/ − i L/ ) (cid:17) ×× (cid:112) ( n )! ... ( n L )! (cid:12)(cid:12)(cid:12) α (∆); j , ..., j m k , ..., j L/ , ..., j L/ m kL/ (cid:69) , where n k i is the number of particles with momentum k i , N = (cid:81) L/ i =1 (cid:18) m k i + 1 m k i (cid:19) , j , ..., j m k , ..., j L/ , ..., j L/ m kL/ are the positions of the N particles and n i is the oc-cupation number of each site in real space and theeven-odd imbalance is given by ∆ = (cid:80) m k p =1 ( − j p + ... + (cid:80) m kL/ p =1 ( − j L/ p , and the normalization constant is M = (cid:0) L +12 (cid:1) N/ (cid:113) ( n k )!( n π − k )! ... ( n k L/ )!( n π − k L/ )!.
4. Comparison with numerical exact tMPS results
In Figs. (B1-B4) we present additional data comple-menting Fig. 3 comparing the many body adiabatic elimi-nation results taking the kinetic term as the perturbationand the numerical exact tMPS results at large dissipationstrengths.In Fig. B1 and Fig. B2 the time evolution of the scaledphoton number, (cid:104) a † a (cid:105) /N , and of the conserved quanti-ties, (cid:104)O k j (cid:105) , are shown for different values of the inter- action strength for two additional symmetry sectors areplotted. The dissipation strength has been chosen large,Γ /J = 15, such that we can compare to the results of themany body adiabatic elimination with the kinetic term asthe perturbation. As for the symmetry sector presentedin Fig. 3, we see that, at finite interaction, the late timebehavior is nicely described by an exponential decay to-wards the many body adiabatic elimination state, ρ mix [see Eq. (B.11)]. We capture this by performing an ex-ponential fit for (cid:104)O k j (cid:105) − (cid:104)O k j (cid:105) ρ mix , ∝ e − t/τ . We see thatthe fit work very well in most cases which supports thedecay towards the steady state ρ mix . The decay time τ gives the timescale for reaching the steady state. Thedeviations seen in the curves for the strongest interac-tion are of the order of the statistical uncertainty of theMonte-Carlo sampling of the different trajectories. Thetimescales corresponding to Fig. B1 are shown in Fig. 3,for the photon number, and in Fig. B3, for the conservedquantities, and support the decay of τ ∝ /U . We notethat for the symmetry sector considered in Fig. B1 we donot have enough data to extract the exponent.In Fig. B4 we look at the finite time value at time1 FIG. B1. The time evolution of (a) the scaled photon number, (cid:104) a † a (cid:105) /N , and [(b)-(f)] the expectation value of O k j for differentvalues of U . For finite U we fit the time evolution with an exponential decay (black dashed lines) the difference between thetMPS data and the expected steady state value, obtained from many-body adiabatic elimination. The parameters are chosento be L = 10, N = 5, (cid:126) Ω √ N/J = 4 . (cid:126) δ/J = 2, Γ /J = 15, and the symmetry sector ( m k = 5). tJ = 49 . (cid:126) of the scaled photon number, (cid:104) a † a (cid:105) /N , andall conserved quantities, (cid:104)O k j (cid:105) , as a function of the inter-action strength for different symmetry sectors. Part ofthe data presented is overlapping with the data in Fig. 3(a) and (b), but here additional symmetry sectors areshown and all conserved quantities (cid:104)O k j (cid:105) .In the presence of the strong symmetry, at U = 0(marked by a dashed vertical line in Fig. B4), we comparethe tMPS results with the expectations value computedwith the state ρ K,st , Eq. (B.18), and we obtain very good agreement in all sectors. We note that for the state ρ K,st the expectation of the photon number depends only onthe distribution of the particles in the single particle sec-tors and not the particular momentum values the par-ticles have, i.e. the sectors ( m k = 5) and ( m k = 5)have the same photon number. At finite interaction weexpect an agreement between the tMPS results and thestate ρ mix , Eq. (B.11), as we can observe for U/J (cid:38) . U we attribute to the fact thatthe numerical results are taken at finite time and not in2 FIG. B2. The time evolution of (a) the scaled photon number, (cid:104) a † a (cid:105) /N , and [(b)-(f)] the expectation value of O k j for differentvalues of U . For finite U we fit the time evolution with an exponential decay (black dashed lines) the difference between thetMPS data and the expected steady state value, obtained from many-body adiabatic elimination. The parameters are chosento be L = 10, N = 5, (cid:126) Ω √ N/J = 4 . (cid:126) δ/J = 2, Γ /J = 15, and the symmetry sector ( m k = 3 , m k = 2). the steady state (see Figs. B1, B2), since the exponentialfit in Fig. B1 and Fig. B2 approaches the correct steadystate value. C. DETAILS OF THE tMPS METHOD FORTHE COUPLED PHOTON-ATOM SYSTEM
The numerically exact results are obtained with a ma-trix product state (MPS) method developed for the sim-ulation of the time evolution of the dissipative master equation, Eqs. (1)-(2) in the main text, for the cavity-atoms coupled systems [30, 43]. The details regardingthe implementation and benchmarking of the methodare presented in Ref. [43]. The method is based on thestochastic unravelling of the master equation with quan-tum trajectories [50–52] and a variant of the quasi-exacttime-dependent variational matrix product state (tMPS)employing the Trotter-Suzuki decomposition of the timeevolution propagator [53–55] and the dynamical defor-mation of the MPS structure using swap gates [55–57].The convergence of our results is sufficient for at3
FIG. B3. The timescales obtained from the exponential fitsof O k j as a function of U , for the data presented in Fig. B1.The lines represent a fit of the timescale dependence on theinteraction with an algebraic decay ∝ U − α , we obtain thefollowing exponents: j = 1, α = 1 . ± . j = 2, α =1 . ± . j = 3, α = 1 . ± . j = 4, α = 1 . ± . j = 1, α = 1 . ± .
04. The parameters are chosen to be L = 10, N = 5, (cid:126) Ω √ N/J = 4 . (cid:126) δ/J = 2, and Γ /J = 15. least 500 quantum trajectories in the Monte Carlo sam-pling, the truncation error goal of 10 − , the time-stepof dtJ = 0 . (cid:126) for the parameters used in Fig. 1and dtJ = 0 . (cid:126) for the other parameters consid-ered in this work, and an adaptive cutoff of the localHilbert space of the photon mode between N pho = 20and N pho = 10. D. DISSIPATIVE FREEZING
In this appendix we describe how the phenomenon ofdissipative freezing occurs for the special case that thesystem satisfies (cid:2)
H, L † L (cid:3) = 0.The time-evolution of a single quantum trajectory isgiven by the following time-evolution operator U ( t, t ) = 1 M e − i ˜ H ( t − t N ) / (cid:126) (cid:89) j = N L e − i ˜ H ( t j − t j − ) / (cid:126) , (D.1)where M is the normalization constant, L the jump op-erator, and { t , ..., t N } are the stochastically sampledtimes when a quantum jump occurs. The effective non-Hermitian Hamiltonian is˜ H = H − i (cid:126) Γ L † L. (D.2)In order to analyze the phenomenon of dissipativefreezing we can look at the evolution of one of the sym-metry generators O k in a single quantum trajectory, (cid:104)O k ( t ) (cid:105) traj = (cid:104) ψ | U † ( t, t ) O k U ( t, t ) | ψ (cid:105) (D.3)= (cid:104) ψ | U † ( t, t ) U ( t, t ) O k | ψ (cid:105) , with | ψ (cid:105) the initial state. If | ψ (cid:105) is within one symmetrysector and is an eigenstate of O k , the expectation valuewithin the single trajectory (cid:104)O k ( t ) (cid:105) traj will not evolvein time, as neither the jump operator, or the Hamilto-nian can change the symmetry sector. In contrast, if theinitial state is taken as a superposition with contribu-tions from different symmetry sectors, then, in principle,both the jump operator or the evolution with the effec-tive Hamiltonian can change the weights of these contri-butions. This implies that (cid:104)O k ( t ) (cid:105) traj will evolve in timeusing a single quantum trajectory and only the MonteCarlo average will be constant.In the case the systems satisfies (cid:2) H, L † L (cid:3) = 0 one canget a better insight as U † ( t, t ) U ( t, t ) = e − Γ L † L ( t − t ) ( L † L ) n , (D.4)with n the number of quantum jumps that occur up totime t . Here the evolution of O k will only depend onthe number of quantum jumps that occur up to time t . This includes the particular case of L † L = I when U † ( t, t ) U ( t, t ) = I and (cid:104)O ( t ) (cid:105) traj is trivially constant intime. The system considered in Ref. [31] is also includedin this situation, as the authors prove that dissipativefreezing always occurs if H ∝ L ∝ O .For the coupled atom-cavity system that we considerin this work, Eqs. (1)-(2) in the main text, the condition (cid:2) H, L † L (cid:3) (cid:54) = 0 is not satisfied and the arguments givenabove are not directly applicable. Nevertheless, we shownumerically that the dissipative freezing occurs even forthis more involved case.In Fig. D1 we extend the data presented in Fig. 2,by plotting the expectation value of all the generators ofthe strong symmetry, (cid:104)O k j (cid:105) , j = 1 ..
5, in time for 1000single trajectories. The initial state is an equal super-position of a state from the sector ( m k = 5) and thesector ( m k = 1 , m k = 1 , m k = 1 , m k = 1 , m k = 1).We can observe in the first column of Fig. D1 that thephenomenon of dissipative freezing can be identified inthe evolution of each of the symmetry generators, as fortimes tJ (cid:38) (cid:126) , all trajectories evolved to one of thetwo symmetry sectors and the Monte Carlo average ofthe trajectories stays constant throughout the followingtime-evolution. If we slightly turn on the interaction andbreak the strong symmetry (see second column of Fig. D1for U/J = 0 .
01) we see that at short and intermediatetime scales the behavior of the quantum trajectories isvery similar to dissipative freezing. Thus we can inferthat the approximate strong symmetry still affects theshort-time dynamics. If we increase the interaction evenfurther,
U/J ≥ .
05, the mixing of the trajectories startsearlier and the dissipative freezing effects are washed out.4
FIG. B4. The dependence on the interaction strength U of (a) the scaled photon number, (cid:104) a † a (cid:105) /N , and [(b)-(f)] the expectationvalue of O k j , j = 1 .. tJ = 49 . (cid:126) and many-body adiabatic elimination (AE), for different symmetrysectors. The symbols identifying each symmetry sector are consistent in all panels. The parameters are chosen to be L = 10, N = 5, (cid:126) Ω √ N/J = 4 . (cid:126) δ/J = 2, and Γ /J = 15. FIG. D1. Time evolution of O k for the single quantum trajectories sampled in the Monte Carlo average for different in-teraction strengths U . The initial state consists in an equal superposition between states for the sectors ( m k = 5) and( m k = 1 , m k = 1 , m k = 1 , m k = 1 , m k = 1). In each panel there are 1000 trajectories plotted, the black represent theMonte Carlo average, either for the full set of trajectories, or averaged separately depending on the final value, we shade theinterval of one standard deviation away from the average, with light blue for the full average and light gray for the separateaverages. The parameters used are L = 10, N = 5, (cid:126) δ/J = 2, (cid:126) Ω √ N/J = 4 .
47, and (cid:126) Γ /J = 15. The standard deviationis defined as σ ( O k ( t )) = (cid:113) R (cid:80) Rr =1 ( (cid:104) ψ r ( t ) | O k | ψ r ( t ) (cid:105) − (cid:104)(cid:104)O k (cid:105)(cid:105) ) , where R is the total number of trajectories, | ψ r ( t ) (cid:105) thetime-evolved wave function of the trajectory labeled by r and (cid:104)(cid:104)O k (cid:105)(cid:105)(cid:105)(cid:105)