Dissipation-induced anomalous multicritical phenomena
Matteo Soriente Tobias Donner R. Chitra Oded Zilberberg
cond-mat.quant-gas cond-mat.str-el quant-ph




Tr
tr
Re
Im





†

α_
α_
α_0^re
α_0^im
β_
β_








Institute for Theoretical Physics, ETH Zurich, 8093 Zürich, Switzerland

Institute for Quantum Electronics, ETH Zurich, 8093 Zürich, Switzerland

Institute for Theoretical Physics, ETH Zurich, 8093 Zürich, Switzerland

Institute for Theoretical Physics, ETH Zurich, 8093 Zürich, Switzerland



We explore the influence of dissipation on a paradigmatic driven-dissipative model where a collection of two level atoms interact with both quadratures of a quantum cavity mode. The closed system exhibits multiple phase transitions involving discrete and continuous symmetries breaking and all phases culminate in a multicritical point. In the open system, we show that infinitesimal dissipation erases the phase with broken continuous symmetry and radically alters the model's phase diagram. The multicritical point now becomes brittle and splits into two tricritical points where first- and second-order symmetry-breaking transitions meet. A quantum fluctuations analysis shows that, surprisingly, the tricritical points exhibit anomalous finite fluctuations, as opposed to standard tricritical points arising in ^3He-^4He mixtures. Our work has direct implications for a variety of fields, including cold atoms and ions in optical cavities, circuit-quantum electrodynamics as well as optomechanical systems.


Dissipation-induced anomalous multicritical phenomena
    O. Zilberberg
    December 30, 2023
=====================================================

Dissipation can fundamentally influence quantum many-body systems and their phase transitions in often counter-intuitive ways. Prime examples of open quantum many-body systems are interacting light-matter systems where state of the art experiments are able to engineer dissipation channels <cit.>. They combine in a unique manner the many-body physics of condensed matter systems with quantum optical tools, including driving and well-controlled dissipation <cit.>. In addition to fostering deeper understanding of cooperative phenomena, these systems have potential applications in the realms of quantum computation <cit.>, laser and maser technologies <cit.>, and can potentially generate new states of matter such as light-induced superconductivity <cit.>. New universality classes emerge in driven-dissipative systems <cit.>, and dissipation can generate topological effects <cit.>. Concurrently, due to the inapplicability of the framework of equilibrium statistical physics, our understanding of driven-dissipative models remains limited, motivating further studies in this field.




 


Paradigmatic models of driven-dissipative light-matter systems involve multiple spin-like degrees of freedom that are driven and strongly coupled to bosonic cavity modes (see Fig. <ref>(a)) <cit.>. Such models commonly exhibit quantum phase transitions (QPT) from a normal (NP) to superradiant phases (SP) depending on the coupling between the spins and cavities. Controlled realizations of such models can be achieved in cold atomic quantum gases in high finesse optical cavities <cit.>. These engineered systems can be used to study quantum phase transitions both in and out of equilibrium in an extremely controlled manner, e.g.,  (i) a Z_2 QPT in the so-called driven-dissipative Dicke-model with a single cavity <cit.>, (ii) U(1) supersolid symmetry breaking when coupling two cavity modes to the atoms <cit.>, and (iii) coupling of high-spin atoms to a cavity mode <cit.>. The latter two cases manifest rich phase diagrams where multiple broken symmetry phases meet at a multicritical point. Similar features are also shared by other fundamental models such as the Lipkin-Meshkov-Glick (LMG) model <cit.>, or models with QPTs of tunable symmetries <cit.>. Interestingly, already for the Dicke model, the inevitable coupling to dissipation channels 
was shown to alter the closed system physics <cit.>. In particular, cavity dissipation was shown to lead to shifts of the critical points and to modifications of critical exponents <cit.>.





Here, we study the connection between the closed and open phase diagram of a paradigmatic driven-dissipative model <cit.>, hosting a multicritical point and phase transitions breaking discrete and continuous symmetries. Infinitesimal dissipation  dramatically impacts the  model's phase diagram, resulting in rich phenomena, including a splitting of a multicritical point into two tricritical points, coexistence of phases <cit.>, and relics of the continuous symmetry in rotated order parameters.  In particular,  we analyze the model's quantum fluctuations  and show that these tricritical points  exhibit anomalous finite fluctuations, as opposed to standard tricritical points  <cit.>. 

We consider a bosonic cavity mode coupled to N two-level systems described by the Hamiltonian <cit.>

    H     = ħω_c a^† a + ħω_a  S_z +
       + 2ħλ_x/√(N) S_x (a + a^†)+ 2ħλ_y/√(N) i S_y (a - a^†) ,

where a^† and a are the bosonic creation and annihilation operators of the cavity field, respectively, and the cavity's resonance frequency is ω_c. The collective spin operators S_α = ∑_j=1^N σ_α^j  with α=x,y,z are constituted from the individual Pauli spin operators σ_α^j  describing the identical two-level systems with level spacing ħω_a. The two quadratures of the cavity field couple to different projections of the collective spin operators with arbitrary real couplings λ_x and λ_y. Hence, this model interpolates between two ubiquitous light-matter models, the Dicke model <cit.> limit when either λ_x=0 or λ_y=0 and the Tavis-Cummings model <cit.> for λ_x=λ_y, and will be dubbed here as the interpolating Dicke-Tavis-Cummings (IDTC) model. A schematic illustration of the model is depicted in Fig. <ref>(a). 




The IDTC model has a Z_2×Z_2 symmetry, except for along the diagonal λ_x=λ_y where it has an enlarged U(1) symmetry, as schematically shown in Fig. <ref>(b) <cit.>. For λ_x,λ_y≤λ_c≡√(ω_cω_a)/2, the system has a trivial ground state, dubbed normal phase (NP), that is comprised of an empty cavity and all two-level system in their respective ground states, oriented along the z-axis. Fixing one of the couplings below criticality λ_i≤λ_c and taking the other above it λ_j> λ_c with i,j∈{x,y} and i≠ j, the system undergoes a Z_2-breaking Dicke-like transition from the NP to a superradiant phase (SP), where the cavity features a finite mean population, a≠ 0, and the two-level systems are on average oriented away from the z-axis. Beyond the SP threshold and along the diagonal, λ_x=λ_y, the U(1) symmetry is spontaneously broken marking the Tavis-Cummings QPT <cit.>. The hallmark of the IDTC model is the appearance of a multicritical point  at λ_x=λ_y = λ_c, where all phases meet and the symmetry of the Hamiltonian changes from a discrete to a continuous symmetry <cit.>.

The relevant complex order parameter which captures these transitions is ⟨ a ⟩ = √(N)α, where α= + i. At a Dicke-like phase transition either ≠ 0 and = 0, or = 0 and ≠ 0. Along the Tavis-Cummings line, the U(1)-symmetry is broken and both , ≠ 0. The Ginzburg-Landau  energy potential of the order parameter is schematically plotted in Fig. <ref>(b) and was calculated in Ref. <cit.>. It shows (i) a single minimum in the normal phase, (ii) two minima along either  the real or the imaginary axis marking the Z_2 Dicke-like symmetry breaking, and (iii) a “sombrero-hat” potential with an enlarged U(1) symmetry on the diagonal.

The closed system phase diagram [Fig. <ref>(b)], however, fundamentally changes if one includes dissipation in the model, which will be relevant for any experimental implementation of such a system. In the presence of both cavity and global atomic dissipation, the driven and dissipative nature of the system is described by a Liouvillian equation for the density matrix ρ_ sys of the system <cit.>

    d ρ_ sys/dt   =    - i/ħ [H(t), ρ_ sys]  + κ[ 2 a ρ_ sys a^† - { a^† a , ρ_ sys}]  
          + 
    γ/N [ 2 S_- ρ_ sys S_+ -{S_+ S_-, ρ_ sys}] ,

where S_±=S_x± iS_y are ladder operators. The first term on the r.h.s. describes the standard Hamiltonian evolution while the other two terms represent the Markovian dissipation for both the cavity and the collective spin in Lindblad form with rates κ and γ, respectively. It should be noted that the master equation (<ref>) is valid in the rotating frame of driven systems with weak cavity-spin coupling, cf. Refs. <cit.>. The ultrastrong coupling regime  should  generally be investigated using  dressed operators <cit.>. Here, we use Eq. (<ref>) and discuss the validity of our results in the strong coupling regime in <cit.>. In the following, we set γ=0 and focus mainly on cavity dissipation. 

 
The mean-field equations governing the different ordered phases can be derived from Eq. (<ref>),

    ω_c - κ -2λ_y Y    = 0 , 
    ω_c  +κ +  2 λ_x X     = 0 ,
    ω_a Y + 4 λ_y   Z     = 0 , 
    ω_a X - 4 λ_x  Z     = 0 ,

where we defined ⟨ S_x⟩ = NX, ⟨ S_y⟩ = NY and ⟨ S_z⟩ = NZ and have taken the steady-state limit. Additionally, for the case of the global dissipation considered here, we have the spin-conservation law X^2 + Y^2 + Z^2 = 1/4 that is used to solve the mean-field equations analytically <cit.>.

Solving Eqs. (<ref>)-(<ref>), we find that dissipation stabilizes a richer phenomenology inaccessible in the closed system paradigm, with multiple bifurcations  and coexisting many-body phases, see Fig. <ref>(c) <cit.>. Specifically, the mean-field equations can be manipulated to obtain an equation for Z: κ^2 ω_a^2+ [(8 Z λ_x^2 + ω_c ω_a) (8 Z λ_y^2 + ω_c ω_a)] =0 <cit.>. All valid solutions are then subject to a constraint 4 κ^2 λ_x^2 λ_y^2 ≤ (λ_x^2 -λ_y^2)^2ω_c^2 <cit.>. However, only a subset of these solutions are compatible with the aforementioned spin-conservation law. Note that, the spin conservation is trivially violated in the presence of single spin dephasing leading to the destruction of Dicke superradiant phases, <cit.>, which is beyond the scope of this work. 


The implications of this constraint are numerous, for example, (i) it is violated at the Tavis-Cummings line, λ_x=λ_y, implying that the U(1) QPT is destroyed by dissipation, as can be easily seen using adiabatic elimination <cit.> [cf. Ref. <cit.>]; (ii) it is additionally violated in a κ-dependent sliver around the Tavis-Cummings line, see Fig. <ref>(c). Consequently, the NP percolates through this sliver beyond the multicritical region of the closed system. Away from the sliver, solving the mean-field equations leads to multiple stable and unstable solutions <cit.>. In particular, the independent Dicke-like phases boast two stable solutions corresponding to the Z_2-broken SP states, as well as an unstable NP solution, a feature  seen  also in  <cit.>. Interestingly, in the regime where both couplings are above λ_c, coexistence of both stable SP and NP states appears, supported by two unstable solutions, see Fig. <ref>(c) and <cit.>. 

Dissipation leads to additional important features in the IDTC as evinced by Fig. <ref> where the order parameters,  and , for one of the stable  symmetry broken solutions are  plotted. Clearly,   and  leak across the NP diagonal sliver, unlike the nondissipative case [cf. Fig. <ref>]. Consequently, the cavity field α is complex within the Z_2 SP regions, which can be attributed to  the remnant memory of the underlying U(1) symmetry in the problem. Similarly,  changes sign as a function of λ_x and λ_y, thus showing that the mean field solutions are rotating within the complex plane. Importantly, the order parameter components,  and  evolve continuously from zero across the Dicke-like phase boundaries, but show a discontinuous behavior along the dissipation-induced NP sliver. Therefore, along the sliver edges, there are two second-order transition lines, which morph into two first-order transition lines in the vicinity of the nondissipative multicritical point, λ_x=λ_y=λ_c, i.e., these lines meet at new dissipation-induced tricritical points,  which separate the continuous and discontinuous symmetry breaking transitions in the system, see Figs. <ref>(c) and <ref>. For ω_a=ω_c=ω, we find that the separation between the two out-of-equilibrium tricritical points scales linearly with κ for κ<<ω and is given by δ = ω√(1 + ( κ/ω) ^2 -√(1 + ( κ/ω) ^2))≈κ/√(2) <cit.>.






To summarize, the dissipation renders the multicriticality of the IDTC model brittle. Nevertheless, signatures of the Tavis-Cummings U(1) symmetry manifest in the splitting of the multicritical point into two new tricritical points with coexisting phases, and the two closed-system SP phases are separated by an emerging NP sliver. Importantly, the multicriticality is sensitive to even infinitesimally small cavity dissipation. This is radically different from standard driven-dissipative systems such as the Dicke model, where cavity dissipation or a global atomic dissipation engender a small modification of phase boundaries <cit.>.  


A more thorough characterization of the phase diagram is provided by an analysis of the steady state fluctuations and their scaling beyond mean field.  In all regimes of the parameter space, the fluctuation Hamiltonian reduces to a problem of two coupled linear oscillators <cit.>. In the thermodynamic limit (large-N), it takes the generic form <cit.>

    H_ fl =   ħω_c c^† c +ħΩ_a d^† d 
       +(Γ_1 c^† d + Γ_2 c d +Γ_3  d^†^2 + h.c.) ,

where we have defined a=α√(N)+c with c the bosonic cavity fluctuation operator, and have used the Holstein-Primakoff representation for the spins, S_+=  b^†√(N - b^† b) and S_z= -N/2 + b^† b. Here, b=β√(N)+d is a bosonic annihilation operator describing deviations away from the z-axis, with β=(X-iY)/√(1/2 - Z) its mean value and d the atomic fluctuation operator. The parameters Ω_a, Γ_1, Γ_2 and Γ_3 are determined by the mean-field solutions <cit.>.

A normal mode transformation on Eq. (<ref>) yields the excitation spectra of the problem as was studied for the closed system case in Ref. <cit.>. To analyze the nature of the phase transitions in the open system dynamics, we calculate, using the Liouvillian (<ref>), the time evolution of the equal-time two-operator correlation functions. The resulting equations of motion form a closed set of ten coupled differential equations for the correlators, ⟨ c^† c⟩,  ⟨ c^† c^†⟩,  ⟨ c c⟩,  ⟨ d^† d⟩, ⟨ d^† d^†⟩, ⟨ dd⟩, ⟨ c^† d⟩,  ⟨ c^† d^†⟩,  ⟨ c d⟩, and ⟨ d^† c⟩ <cit.>. In the steady state, this set conforms to a set of linear inhomogeneous equations, which can then be solved for λ_x λ_y to obtain the various correlators <cit.>. All correlators display similar features and in the following we shall focus on the photon number correlator ⟨ c^†  c⟩, see Fig. <ref>.

It is instructive to study the expectation value for the photon number in the normal phase. For ω_a=ω_c=ω, it is given by 
 
    ⟨ c^†  c⟩=   (λ_x^2 - λ_y^2)^2/2[κ^2 λ_xλ_y + (λ_x+λ_y)^2ω^2]
       ×ω^2 ( ω^2 + κ^2 + 4 λ_x λ_y)/[16 λ_x^2λ_y^2 +ω^2(ω^2 +κ^2 - 4(λ_x^2+λ_y^2))] .

In the Dicke-limit, where λ_i →λ_c while λ_j=0 with i≠ j, the photon number diverges with an exponent of 1, as expected <cit.>. Note, however, that away from the Dicke-limit, the critical regime of the Dicke-like transitions shrinks. This is consistent with the fact that the multicritical point is expected to have zero fluctuations <cit.>. For the case λ_x=λ_y, we find that the set of ten equations is no longer invertible. However, a reduced solvable set exists for which the photon fluctuations vanish. This is a remarkable result which shows that even weak cavity dissipation which technically preserves the length of the spin destroys an ordered phase.

The specific features of the QPTs in the IDTC can be inferred by considering photon number fluctuations along three representative cuts in the λ_x- λ_y plane shown  in Fig. <ref>(a). Along I, where at least one of the couplings is below λ_c (top panel of Fig. <ref>(b)), as expected for a continuous  Dicke-like QPT, cavity fluctuations diverge continuously across both sides of the transition with critical exponent 1. Across cut II (bottom panel of Fig. <ref>(b)),  where  either λ_x or λ_y is greater than λ_c, the loss of the broken U(1) phase results in discontinuous first-order transitions between the NP and the two SP phases concomitant with coexistence regions. Cavity fluctuations, though enhanced, remain finite and exhibit a jump across the phase boundaries. The widths of the coexistence regimes effectively indicate the size of hysteresis loops that will appear under scans of the couplings. In contrast to standard tricritical points that arise in systems as diverse as ^3He-^4He mixtures <cit.> and high-Tc superconducting vortex lattices <cit.>, the fluctuations at the two out-of-equilibrium tricritical points of the IDTC remain finite. The scaling of these fluctuations with dissipation is illustrated in Fig. <ref>(c) for cut III, which shows that fluctuations diminish as the two tricritical points approach the original multicritical point.


We have shown that weak dissipation can dramatically alter the paradigm of standard continuous symmetry breaking phase transitions in a model system exhibiting multicriticality. Additionally, the dissipation induced tricritical points are characterized by anomalous quantum fluctuations. We expect our results to be qualitatively valid for other dissipation channels provided the spin is conserved. Extending our work to variants of the IDTC with (i) higher-spin systems, (ii) negatively detuned frequencies where interesting oscillatory behavior are expected [cf. Ref. <cit.>], (iii) additional cavity fields [cf. Ref. <cit.>], (iv) dephasing and non spin-conserving dissipation, and (v) quenched dynamics will further reinforce our predictions for existing experiments <cit.>. Our work also motivates a study of potential brittle multicritical phenomena in quantum engineered systems, and in out-of-equilibrium matter systems as well as the influence of non-Markovian noise on such phase diagrams.



  






We would like to thank I. Carusotto,  A. Imamoglu, J. Keeling, and M. Landini for useful discussions. We acknowledge financial support from the Swiss National Science Foundation (SNSF), Division 2 and through the SNSF DACH-project “Quantum Crystals of Matter and Light”.









42
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Syassen et al.(2008)Syassen, Bauer, Lettner, Volz, Dietze, García-Ripoll,
  Cirac, Rempe, and Dürr]Syassen1329
  author author N. Syassen, author D. M. Bauer,
  author M. Lettner, author T. Volz, author
  D. Dietze, author J. J. García-Ripoll, author J. I. Cirac, author G. Rempe,  and author S. Dürr, 10.1126/science.1155309
  journal journal Science volume 320, pages 1329 (year
  2008)NoStop
[Baumann et al.(2010)Baumann, Guerlin, Brennecke, and  Esslinger]Baumann2010
  author author K. Baumann, author C. Guerlin,
  author F. Brennecke,  and  author T. Esslinger, http://dx.doi.org/10.1038/nature09009 journal journal Nature volume 464, pages
  1301 (year 2010)NoStop
[Barreiro et al.(2011)Barreiro, Müller, Schindler,
  Nigg, Monz, Chwalla,
  Hennrich, Roos, Zoller, and  Blatt]Barreiro2011
  author author J. T. Barreiro, author M. Müller, author P. Schindler, author D. Nigg,
  author T. Monz, author
  M. Chwalla, author M. Hennrich, author C. F.  Roos, author P. Zoller,  and author R. Blatt, http://dx.doi.org/10.1038/nature09801 journal journal Nature volume 470,  pages 486 (year 2011)NoStop
[Barontini et al.(2013)Barontini, Labouvie, Stubenrauch,
  Vogler, Guarrera, and Ott]PhysRevLett.110.035302
  author author G. Barontini, author R. Labouvie,
  author F. Stubenrauch, author A. Vogler, author
  V. Guarrera,  and author
  H. Ott, 10.1103/PhysRevLett.110.035302 journal journal
  Phys. Rev. Lett. volume 110, pages
  035302 (year 2013)NoStop
[Carusotto and Ciuti(2013)]Carusotto2013
  author author I. Carusotto and author C. Ciuti, 10.1103/RevModPhys.85.299 journal
  journal Rev. Mod. Phys. volume 85,  pages 299 (year 2013)NoStop
[Hartmann et al.(2008)Hartmann, Brandão, and Plenio]Hartmann2008
  author author M. Hartmann, author F. Brandão,  and author M. Plenio, 10.1002/lpor.200810046 journal
  journal Laser & Photonics Review volume 2, pages 527 (year 2008)NoStop
[Duan et al.(2001)Duan,
  Lukin, Cirac, and Zoller]Zoller2001
  author author L. M. Duan, author M. D. Lukin,
  author J. I. Cirac,  and  author P. Zoller, http://dx.doi.org/10.1038/35106500 journal journal Nature volume 414, pages
  413 (year 2001)NoStop
[Bohnet et al.(2012)Bohnet,
  Chen, Weiner, Meiser,
  Holland, and Thompson]Bohnet2012
  author author J. G. Bohnet, author Z. Chen,
  author J. M. Weiner, author D. Meiser, author
  M. J. Holland,  and author
  J. K. Thompson, http://dx.doi.org/10.1038/nature10920 journal journal Nature volume 484, pages
  78 (year 2012)NoStop
[Fausti et al.(2011)Fausti,
  Tobey, Dean, Kaiser,
  Dienst, Hoffmann, Pyon,
  Takayama, Takagi, and Cavalleri]cavalleri
  author author D. Fausti, author R. I. Tobey,
  author N. Dean, author
  S. Kaiser, author A. Dienst, author M. C.  Hoffmann, author S. Pyon, author T. Takayama,
  author H. Takagi,  and author A. Cavalleri, 10.1126/science.1197294 journal journal
  Science volume 331, pages 189
  (year 2011)NoStop
[Tomadin et al.(2011)Tomadin, Diehl, and Zoller]Tomadin2011
  author author A. Tomadin, author S. Diehl,   and author P. Zoller, 10.1103/PhysRevA.83.013611 journal journal Phys. Rev. A volume 83, pages 013611 (year 2011)NoStop
[Diehl et al.(2010)Diehl,
  Tomadin, Micheli, Fazio,  and Zoller]Diehl2010
  author author S. Diehl, author A. Tomadin,
  author A. Micheli, author R. Fazio,  and author
  P. Zoller, 10.1103/PhysRevLett.105.015702 journal journal
  Phys. Rev. Lett. volume 105, pages
  015702 (year 2010)NoStop
[Nagy et al.(2011)Nagy,
  Szirmai, and Domokos]Nagy2011
  author author D. Nagy, author G. Szirmai,   and author P. Domokos, 10.1103/PhysRevA.84.043637 journal journal Phys. Rev. A volume 84, pages 043637 (year 2011)NoStop
[Diehl et al.(2011)Diehl,
  Rico, Baranov, and Zoller]Diehl2011
  author author S. Diehl, author E. Rico,
  author M. A. Baranov,  and  author P. Zoller, http://dx.doi.org/10.1038/nphys2106 journal journal Nature Physics volume 7, pages 971 (year 2011)NoStop
[Zhiqiang et al.(2017)Zhiqiang, Lee, Kumar, Arnold, Masson, Parkins, and  Barrett]Barrett2017
  author author Z. Zhiqiang, author C. H. Lee,
  author R. Kumar, author K. J. Arnold, author
  S. J. Masson, author
  A. S. Parkins,  and author
  M. D. Barrett, 10.1364/OPTICA.4.000424 journal journal Optica  volume 4, pages 424 (year
  2017)NoStop
[Baksic and Ciuti(2014)]Baksic2014
  author author A. Baksic and author C. Ciuti,  10.1103/PhysRevLett.112.173601 journal
  journal Phys. Rev. Lett. volume 112,  pages 173601 (year 2014)NoStop
[sup()]supmat
  @noop  note See Supplemental Material for additional
  details, which includes Ref. <cit.>NoStop
[Dimer et al.(2007)Dimer,
  Estienne, Parkins, and Carmichael]Dimer
  author author F. Dimer, author B. Estienne,
  author A. S. Parkins,  and  author H. J. Carmichael, 10.1103/PhysRevA.75.013804 journal journal Phys. Rev. A volume 75, pages 013804 (year 2007)NoStop
[Nagy et al.(2010)Nagy,
  Kónya, Szirmai, and Domokos]Domokos
  author author D. Nagy, author G. Kónya,
  author G. Szirmai,  and  author P. Domokos, 10.1103/PhysRevLett.104.130401 journal journal Phys. Rev. Lett. volume 104, pages 130401 (year 2010)NoStop
[Chirolli et al.(2012)Chirolli, Polini, Giovannetti, and  MacDonald]MacDonald
  author author L. Chirolli, author M. Polini,
  author V. Giovannetti,  and  author A. H. MacDonald, 10.1103/PhysRevLett.109.267404 journal journal Phys. Rev. Lett. volume 109, pages 267404 (year 2012)NoStop
[Klinder et al.(2015)Klinder, Keßler, Wolke, Mathey, and Hemmerich]Klinder2015
  author author J. Klinder, author H. Keßler,
  author M. Wolke, author L. Mathey,  and author
  A. Hemmerich, 10.1073/pnas.1417132112 journal journal
  Proceedings of the National Academy of Sciences volume 112, pages 3290 (year
  2015)NoStop
[Léonard et al.(2017)Léonard, Morales, Zupancic,
  Esslinger, and Donner]Leonard2017
  author author J. Léonard, author A. Morales, author P. Zupancic,
  author T. Esslinger,  and  author T. Donner, http://dx.doi.org/10.1038/nature21067 journal journal Nature volume 543, pages
  87 (year 2017)NoStop
[Morrison and Parkins(2008)]Morrison2008
  author author S. Morrison and author A. S. Parkins, 10.1103/PhysRevLett.100.040403 journal journal Phys. Rev. Lett. volume 100, pages 040403 (year
  2008)NoStop
[Fan et al.(2014)Fan,
  Yang, Zhang, Ma,
  Chen, and Jia]Fan2014
  author author J. Fan, author Z. Yang, author Y. Zhang, author
  J. Ma, author G. Chen,  and author S. Jia, 10.1103/PhysRevA.89.023812
  journal journal Phys. Rev. A volume 89, pages 023812 (year
  2014)NoStop
[Kirton and Keeling(2017)]Keeling2017
  author author P. Kirton and author J. Keeling, 10.1103/PhysRevLett.118.123602 journal journal Phys. Rev. Lett. volume 118, pages 123602 (year
  2017)NoStop
[Larson and Irish(2017)]Larson2017
  author author J. Larson and author E. K. Irish, http://stacks.iop.org/1751-8121/50/i=17/a=174002
  journal journal Journal of Physics A:
  Mathematical and Theoretical volume 50, pages 174002 (year 2017)NoStop
[Nagy and Domokos(2015)]Nagy2015
  author author D. Nagy and author P. Domokos,  10.1103/PhysRevLett.115.043601 journal
  journal Phys. Rev. Lett. volume 115,  pages 043601 (year 2015)NoStop
[Kulkarni et al.(2013)Kulkarni, Öztop, and Türeci]Kulkarni2013
  author author M. Kulkarni, author B. Öztop,   and author H. E. Türeci,  10.1103/PhysRevLett.111.220408 journal
  journal Phys. Rev. Lett. volume 111,  pages 220408 (year 2013)NoStop
[Brennecke et al.(2013)Brennecke, Mottl, Baumann, Landig, Donner, and Esslinger]Brennecke2013
  author author F. Brennecke, author R. Mottl,
  author K. Baumann, author R. Landig, author
  T. Donner,  and author
  T. Esslinger, 10.1073/pnas.1306993110 journal journal
  Proceedings of the National Academy of Sciences volume 110, pages 11763 (year
  2013)NoStop
[Keeling et al.(2010)Keeling, Bhaseen, and Simons]Keeling2010
  author author J. Keeling, author M. J. Bhaseen,  and author B. D. Simons, 10.1103/PhysRevLett.105.043001 journal journal Phys. Rev. Lett. volume 105, pages 043001 (year
  2010)NoStop
[Griffiths(1970)]griffiths
  author author R. B. Griffiths, https://link.aps.org/doi/10.1103/PhysRevLett.24.715
  journal journal Physical Review Letters  volume 24, pages 715 (year
  1970)NoStop
[Dicke(1954)]Dicke
  author author R. H. Dicke, 10.1103/PhysRev.93.99 journal
  journal Phys. Rev. volume 93,  pages 99 (year 1954)NoStop
[Emary and Brandes(2003)]Emary2003
  author author C. Emary and author T. Brandes, 10.1103/PhysRevE.67.066203 journal journal Phys. Rev. E volume
  67, pages 066203 (year 2003)NoStop
[Tavis and Cummings(1968)]TavisCummings
  author author M. Tavis and author F. W. Cummings, 10.1103/PhysRev.170.379 journal
  journal Phys. Rev. volume 170,  pages 379 (year 1968)NoStop
[Vogel and Welsch(2006)]vogel2006
  author author W. Vogel and author D.-G. Welsch, @noop  title Quantum optics  (publisher John Wiley & Sons, year
  2006)NoStop
[Shirai et al.(2014)Shirai,
  Mori, and Miyashita]Miyashita2014
  author author T. Shirai, author T. Mori,  and  author S. Miyashita, http://stacks.iop.org/0953-4075/47/i=2/a=025501 journal
  journal Journal of Physics B: Atomic, Molecular and Optical
  Physics volume 47, pages 025501
  (year 2014)NoStop
[Lolli et al.(2015)Lolli,
  Baksic, Nagy, Manucharyan,  and Ciuti]Lolli2015
  author author J. Lolli, author A. Baksic,
  author D. Nagy, author
  V. E. Manucharyan,  and author
  C. Ciuti, 10.1103/PhysRevLett.114.183601 journal journal
  Phys. Rev. Lett. volume 114, pages
  183601 (year 2015)NoStop
[dre()]dressed
  @noop  note The master equation description
  [Eq. (<ref>)] is valid in the Born-Markov limit of weak
  dissipation. In the limit of strong cavity-atom coupling it is advisable to
  use dressed operators for the Lindblad dissipators <cit.>,
  which introduce new dissipation channels, e.g., spin decay and multi-photon
  processes. Nevertheless, the sensitivity of the Tavis-Cummings line to
  infinitesimal cavity-dissipation hints to the fact that dissipation in the
  ultra-strong-coupling regime would only affect our reported phenomena in a
  quantitative way.Stop
[Dalla Torre et al.(2016)Dalla Torre, Shchadilova, Wilner,
  Lukin, and Demler]Torre2016
  author author E. G. Dalla Torre, author Y. Shchadilova, author E. Y. Wilner, author M. D. Lukin,   and author E. Demler, 10.1103/PhysRevA.94.061802 journal journal Phys. Rev. A volume 94, pages 061802 (year 2016)NoStop
[Öztop et al.(2012)Öztop, Bordyuh, Müstecaplıoğlu, and Türeci]Oztop2012
  author author B. Öztop, author M. Bordyuh,
  author Ö. E. Müstecaplıoğlu,  and author H. E.  Türeci, http://stacks.iop.org/1367-2630/14/i=8/a=085011 journal
  journal New Journal of Physics volume
  14, pages 085011 (year 2012)NoStop
[Konczykowski et al.(2006)Konczykowski, van der Beek, Koshelev,
  Mosser, Dodgson, and Kes]vortex
  author author M. Konczykowski, author C. J. van der Beek, author A. E. Koshelev, author V. Mosser,
  author M. Dodgson,  and  author P. H. Kes, 10.1103/PhysRevLett.97.237005 journal journal Phys. Rev. Lett. volume 97, pages 237005 (year 2006)NoStop
[Moodie et al.(2018)Moodie,
  Ballantine, and Keeling]Keeling2018
  author author R. I. Moodie, author K. E. Ballantine,  and author J. Keeling, 10.1103/PhysRevA.97.033802 journal journal Phys. Rev. A volume
  97, pages 033802 (year 2018)NoStop
[Bhaseen et al.(2012)Bhaseen, Mayoh, Simons, and Keeling]Bhaseen2012
  author author M. J. Bhaseen, author J. Mayoh,
  author B. D. Simons,  and  author J. Keeling, 10.1103/PhysRevA.85.013817 journal journal Phys. Rev. A volume 85, pages 013817 (year 2012)NoStop









Supplemental Material for:



Dissipation-induced anomalous multicritical phenomena



 Matteo Soriente,^1 Tobias Donner,^2 R. Chitra,^1 and Oded Zilberberg^1



^1Institute for Theoretical Physics, ETH Zürich, 8093 Zürich, Switzerland

^2Institute for Quantum Electronics, ETH Zurich, 8093 Zürich, Switzerland























§ I. ADIABATIC ELIMINATION
  
For the case of symmetric coupling (λ = λ_x = λ_y), the Hamiltonian can be written as

    H =    -ħω_c a^† a + ħω_a S_z + λ (a S_+ + a^† S_-) .

From this we derive the equations of motion for the expectation values of the operators S_z and a including cavity decay at rate κ,

    ⟨Ṡ_z⟩=    2i λ( ⟨ a^†⟩⟨ S_-⟩ - ⟨ a⟩⟨ S_+⟩) ,
    ⟨ȧ⟩ =    i ω_c ⟨ a⟩ - i λ⟨ S_-⟩ - κ⟨ a⟩ ,

where the factorization ⟨ a S_i⟩ = ⟨ a⟩⟨ S_i⟩ has been imposed.
Consdering the steady state for the cavity field (⟨̇ ̇ȧ⟩̇=0) gives

    ⟨ a⟩ =    -2λ⟨ S_-⟩/ω_c - i κ
    ⟨ a^†⟩ =    -2λ⟨ S_+⟩/ω_c + i κ .


We then look at the steady state for the spin population, i.e. ⟨Ṡ_z⟩=0, which results in

    ⟨ a⟩⟨ S_+⟩ = ⟨ a^†⟩⟨ S_-⟩ ,

 i.e., in detailed balance for the excitations in the system.
Making use of the expressions (<ref>) and (<ref>) for the light fields, we find

    ⟨ S_-⟩⟨ S_+⟩/ω_c-iκ =   ⟨ S_-⟩⟨ S_+⟩/ω_c + iκ .

Any finite cavity decay κ≠ 0 thus gives rise to the breakdown of detailed balance, leaving only the trivial solution ⟨ a⟩=0, and therefore hindering the Tavis-Cummings phase transition.



§ II. MEAN-FIELD SOLUTION
 


In order to solve the mean-field Eqs. (3)-(6) in the main text, we first obtain expressions for X and Y as functions of ,  and Z, for the case of vanishing atomic dissipation (γ=0). Eqs. (3) and (4) lead to

    X=    -ω_c + κ/2 λ_x ,
     Y=   ω_c - κ/2 λ_y ,

while Eqs. (5) and (6) to

    X=   4 λ_x Z /ω_a 
     Y=    -4 λ_y Z /ω_a .

Substituting (<ref>, <ref>) into (<ref>, <ref>), we obtain a homogeneous system for , parametrized by Z. Requiring that the system has non-trivial solutions, we find the constraint on Z presented in the main text. This constraint yields two solutions 

    Z_±=    -(λ_x^2+λ_y^2)ω_c ±√((λ_x^2-λ_y^2)^2ω_c^2-4κ^2λ_x^2λ_y^2)/16λ_x^2λ_y^2ω_a .

We shall find it useful in the following to define Z_±=𝒜_±ω_a. 

We insert the solution Z_±, as well as, Eqs. (<ref>)-(<ref>) into the spin length normalization constraint, X^2+Y^2+Z^2=1/4, and find two inhomogeneous quadratic equations for  and ,

    c=   a_1 ^2 + a_2 ^2  
     c=   b_1 ^2 + b_2  + b_3 ^2

with coefficients depending on which solution of Z-constraint is taken, a_1=16𝒜_±^2λ_x^2, a_2=16𝒜_±^2λ_y^2, b_1=(κ^2/λ_y^2+ω_c^2/λ_x^2)/4, b_2=κω_c(1/λ_x^2-1/λ_y^2)/2, b_3=(κ^2/λ_x^2+ω_c^2/λ_y^2)/4, and c=1/4 - 𝒜_±^2 ω_a^2.

Equations (<ref>) and (<ref>) are solved at once yielding four possible solutions. We note, however, that our procedure has increased the space of solutions, since we have moved from a first-order system for  and  to a second-order one. This means that not all of the obtained solutions are physical. Imposing that the solutions preserve real spin order parameters and that the total spin is conserved, we find that the validity of the solutions is dependent on the coupling strengths λ_x, λ_y. Specifically, we find that two of the solutions are physical for λ_x > λ_y while the other two are physical for λ_y > λ_x. This is readily combined and the 

















solutions can be written as

    =
       1.!±√(c)√(2 a_2^2 b_1 + a_2 b_2^2 - 2 a_2 b_3 (a_1 + b_1) + 2 a_1 b_3^2 + λ_y - λ_x a_2 |b_2| √(b_2^2 - 4 (a_1 - b_1)(a_2 - b_3) )/2 (a_2^2 b_1^2 + a_1^2 b_3^2 + a_1 a_2 (b_2^2 - 2 b_1 b_3)) ) ,
       =b_2^2 - λ_y - λ_x |b_2|√(4 (b_1 - a_1)(a_2 - b_3) + b_2^2)/2 b_2 (a_2 - b_3) .

Remarkably, in the respective regions of validity, the solutions are the opposite of one another, reflecting the underlying broken Z_2 symmetry.




§ III. STABILITY ANALYSIS


 

In the previous section, we have obtained solutions to the mean-field equations of our model [cf. Eqs. (3)-(6) in the main text]. These solutions, however, can be stable or unstable. To check the linear stability of our solutions, we employ standard stability analysis <cit.>. In particular, we find the stability matrix associated with our system of equations and diagonalize it at a given solution. Once we have its eigenvalues, we study their behavior as a function of the coupling strengths, λ_x, λ_y. If at least one eigenvalue has non-vanishing positive real part then the solution is unstable.

More rigorously, we expand the order parameters to linear order in small 
fluctuations, i. e., α = α_0 + δα, X = X_0 + δ X, Y = Y_0 + δ Y, Z = Z_0 + δ Z. Plugging this ansatz into the mean-field equations, we obtain

    δ̇α̇    = - 2 κ(α_0 + δα) - i ω_0(α_0 + δα) +  
          - 2 i λ_x (X_0 + δ X) - 2 λ_y (Y_0 + δ Y) , 
    δ̇ ̇Ẋ    = - ω_z δ Y + 2 i λ_y (α_0 - α_0^*) δ Z +   
           + 2 i λ_y Z_0 (δα - δα^*)  ,
    δ̇ ̇Ẏ    = ω_z δ X - 2 λ_x (α_0 + α_0^*) δ Z 
           - 2 λ_x Z_0 (δα + δα^*) , 
    δ̇ ̇Ż    = - 2 λ_x (α_0 + α_0^*) δ Y - 2 λ_x Y_0 (δα + δα^*) +  
           + 2 i λ_y (α_0 - α_0^*)δ X + 2 i λ_y X_0 (δα - δα^*)  ,

where α_0, X_0, Y_0, Z_0 are the mean-field steady-state order parameters. 

We use the steady-state solutions from the previous section and obtain

    δ̇ =  - κδ + ω_0 δ - 2 λ_y δ Y , 
       δ̇ =  - κδ - ω_0 δ - 2 λ_x δ X , 
       δ̇ ̇Ẋ =  - ω_z δ Y - 4 λ_y δ Z - 4 λ_y Z_0 δ , 
       δ̇ ̇Ẏ =  ω_z δ X - 4 λ_x δ Z - 4 λ_x Z_0 δ , 
       δ̇ ̇Ż =  - 4 λ_x δ Y - 4 λ_x Y_0 δ + 
              - 4 λ_y δ X - 4 λ_y X_0 δ .

To reduce the number of equations of our system, we linearize the spin constraint, X^2 + Y^2 + Z^2 = 1/4, which allows us to get rid of the δ Z dependence:

    δ Z = - X_0 δ X + Y_0 δ Y/Z_0 .

Substituting the linearized condition into the previous equations, we can rewrite them in matrix form as follows










    [    δ̇;    δ̇; δ̇ ̇Ẋ; δ̇ ̇Ẏ ] = 
      
      
      
      
      
      M[   δ;   δ; δ X; δ Y ] ,

where the matrix
    M ≡[                 - κ                 ω_0                   0             - 2 λ_y;               - ω_0                 - κ             - 2 λ_x                   0;                   0         - 4 λ_y Z_0       4 λ_y X_0/Z_0 4 λ_y Y_0/Z_0 - ω_z;         - 4 λ_x Z_0                   0 4 λ_x X_0/Z_0 + ω_z       4 λ_x Y_0/Z_0 ] ,
is the stability matrix. We now diagonalize this matrix and check the positivity condition for the real part of the eigenvalues. The λ_x, λ_y regions for which the real part of an eigenvalue is zero define the stability lines of Fig. 1(c) in the main text.

Performing this analysis, we find that the nontrivial solutions associated with the Z_- constraint are always stable. The ones associated with the Z_+ constraint, even though physical, are never stable. Interestingly, performing the stability analysis for the normal phase, we find that it is stable when both couplings are above threshold.





§ IV. SEPARATION BETWEEN THE TWO OUT-OF-EQUILIBRIUM TRICRITICAL POINTS


 

We first identify the two tricritical points. This can be done in different ways, e.g. considering the constraints in Eq. (<ref>) or the stability conditions, we chose to use the constraint Z_-. We can restrict to the normal phase without loss of generality. The tricritical points are then given by the intersection of two curves, namely
    ( -(λ_x^2+λ_y^2)ω_c + √((λ_x^2-λ_y^2)^2ω_c^2-4κ^2λ_x^2λ_y^2)/16λ_x^2λ_y^2ω_a)^2     = 1/4
    (λ_x^2-λ_y^2)^2ω_c^2-4κ^2λ_x^2λ_y^2     = 0 .
The first one is just the spin constraint in the normal phase, Z^2=1/4, while the second one gives the boundary for which the square root becomes imaginary and therefore non physical. We solve this system for λ_x and λ_y and find that among all the possible solutions the tricritical points are given by

    [ λ_x1; λ_y1 ]    = [ 1/2√(-√(κ ^2 ω_a^2 (κ ^2+ω_c^2))/ω_c+κ ^2 ω_a/ω_c+ω_c ω_a);                   1/2√(κω_a (√(κ ^2+ω_c^2)+κ)/ω_c+ω_c ω_a) ]
    [ λ_x2; λ_y2 ]    = [ 1/2√(√(κ ^2 ω_a^2 (κ ^2+ω_c^2))/ω_c+κ ^2 ω_a/ω_c+ω_c ω_a);                1/2√(ω_a (-κ√(κ ^2+ω_c^2)+κ ^2+ω_c^2)/ω_c) ] .

In order to distinguish which solutions correspond to a tricritical point we required λ_x,λ_y to be positive real numbers for the parameters used in the main text. We can now evaluate the distance between these two points

    δ    = √(|λ_x1-λ_x2|^2 + |λ_y1-λ_y2|^2)
        = √(ω_a (κ ^2/ω_c+ω_c-√(κ ^2+ω_c^2))) ,

that for ω = ω_c = ω_a reduces to

    δ = ω√(1 + κ ^2/ω ^2-√(1 + κ ^2/ω ^2)).


Expanding this expression for small κ values κ≪ω we get

    δ = κ/√(2) + 𝒪(κ ^2) .





§ V. LARGE-N EXPANSION
 


Taking the IDTC model [Eq. (1) in the main text] and applying the Holstein-Primakoff transformation for the spins, S_+=  b^†√(N - b^† b) and S_z= -N/2 + b^† b, we obtain

    H     = ħω_c a^† a + ħω_a  (-N/2 + b^† b) +
       + ħλ_x/√(N)(b^†√(N - b^† b) + √(N - b^† b)b) (a + a^†)
       + ħλ_y/√(N)(b^†√(N - b^† b) - √(N - b^† b)b) (a - a^†) .

We can rewrite the operators in terms of a mean-field classical part plus quantum fluctuations, a=α√(N)+c and b=β√(N)+d. Inserting these expressions into Eq. (<ref>), we can perform a large-N limit expansion, i.e., collect only terms that do not vanish nor diverge in this limit. We obtain a fluctuation Hamiltonian of the form of Eq. (7) in the main text, where the prefactors are


    Ω_a    =ω_a - λ_x (1+ 3(1-|β|^2)) /(1-|β|^2)^3/2-λ_y (1+ 3(1-|β|^2)) /(1-|β|^2)^3/2
       =ω_a +  2 (-5 + 6 Z) (X λ_x  - Y λ_y )/(1 - 2 Z)^2 ,
    Γ_1   =λ_x(1-|β|^2- β^*)/√(1-|β|^2)+λ_y(1-|β|^2- iβ^*)/√(1-|β|^2)
       =λ_x (- 4 X (X + i Y) + (1 - 2 Z)^2)  + λ_y(- 4 Y (- i X + Y) + (1 - 2 Z)^2)/√(2 - 4 Z)(1 - 2 Z) ,
    Γ_2    =λ_x(1-|β|^2- β^*)/√(1-|β|^2)-λ_y(1-|β|^2- iβ^*)/√(1-|β|^2)
       =λ_x (- 4 X (X + i Y) + (1 - 2 Z)^2)  - λ_y(- 4 Y (- i X + Y) + (1 - 2 Z)^2)/√(2 - 4 Z)(1 - 2 Z) ,
    Γ_3    =  -β[λ_x(2-2|β|^2+ β)/2(1-|β|^2)^3/2+iλ_y(2-2|β|^2- iβ)/2(1-|β|^2)^3/2]
        =- 2 (X - i Y) ((- 2 X Y + i (2 Y^2 + (1 - 2 Z)^2)) λ_y  + (2 X (X - i Y) + (1 - 2 Z)^2) λ_x /(1 - 2 Z)^3 ,


and for the second equality we have used the relation between the classical Holstein-Primakoff spin and the mean-field solutions, β=(X-iY)/√(1/2 - Z).



§ VI. FLUCTUATION EQUATIONS OF MOTION
 


We are interested in the equal-time two-operator correlation functions. Taking, for example, the photon number fluctuation ⟨ c^†(t) c(t)⟩, its time-dependence can be obtained using Heisenberg's equation of motion

    d/dt  c^†(t) c(t) =     -i/ħ [H_ fl, c^† c]  - 2κ c^† c ,

where we have also taken the large-N limit of the Lindblad terms in the Liovillian [Eq. (2) in the main text]. 

Note that the commutation of the photon number with H_ fl couples it to the other two-operator correlation functions. Hence, in order to solve the time-evolution of the photon number, we need to evaluate the Heisenberg's equation of motion for all ten possible two-point correlations functions, and obtain


    d/dt c^† c^ =   - i Γ_1 c^† d^ + i Γ_1^* c^ d^† + i Γ_2 c^ d^ - i c^† d^†Γ_2^* - 2 κ c^† c^ ,
    d/dt  c^ c^ =    - 2 i ω_c c^ c^  -2 i Γ_1 c^ d^ - 2 i Γ_2^* c^ d^†  - 2 κ c^ c^ ,
    d/dt c^† c^† =    2 i ω_c c^† c^† + 2 i Γ_1^* c^† d^† + 2 i Γ_2 c^† d^ - 2 κ c^† c^† ,
    d/dt c^ d^ =    - i (ω_c+ω_a) c^ d^ - i Γ_1 d^ d^ - i Γ_1^* c^ c^ - i Γ_2^*(c^† c^ + d^† d^+1) - 2 i Γ_3 c^ d^† - κ c^ d^ ,
    d/dt c^† d^† =    i (ω_c+ω_a) c^† d^† + i Γ_1^* d^† d^† + i Γ_1 c^† c^† + i Γ_2 (c^† c^ + d^† d^+1) + 2 i Γ_3^* c^† d^ - κ c^† d^† ,
    d/dt c^ d^† =     i c^ d^† (ω_a - ω_c ) + i Γ_1 (c^† c^ - d^† d^) + i Γ_2 c^ c^  - i Γ_2^* d^† d^†  +  2 i Γ_3^* c^ d^  - κ c^ d^† ,
    d/dt c^† d^ =    i c^† d^ (ω_c - ω_a ) - i Γ_1^* (c^† c^ - d^† d^) - i Γ_2^* c^† c^†  + i Γ_2 d^ d^  -  2 i Γ_3 c^† d^†  - κ c^† d^ ,
    d/dt d^ d^ =    - 2 i ω_a d^ d^ - 2 i Γ_1^* c^ d^ - 2 i Γ_2^* c^† d^ - 2 i Γ_3(1 + 2 d^† d^)  ,
    d/dt d^† d^† =    2 i ω_a d^† d^† + 2 i Γ_1 c^† d^† + 2 i Γ_2 c^ d^† + 2 i Γ_3^*(1 + 2 d^† d^) ,
    d/dt d^† d^ =    i Γ_1 c^† d^ - iΓ_1^* c^ d^† + i Γ_2 c^ d^ - i Γ_2^* c^† d^†  - 2 i Γ_3 d^† d^† + 2 i Γ_3^* d^ d^ .


In the steady-state limit, the l.h.s. is set to zero and we obtain a set of ten coupled homogeneous equations. In the IDTC model, for all values λ_x≠ 0 and/or λ_y≠ 0, Γ_1≠ 0 [cf. Eqs. (<ref>)-(<ref>) and (<ref>)]. The set of ten equations in invertible whenever Γ_2≠ 0 and/or Γ_3≠ 0, resulting in an analytical expression for all ten correlators in the steady-state limit. For brevity, we only write the photon number fluctuation


    c^† c=
       1.!2  ( κ^2 + ω_c^2 + |Γ_1|^2 - |Γ_2|^2 ) 
                   [ 2  {Γ_1^2  Γ_2^2  Γ_3^2} 
                     - 2  ω_a  {Γ_1  Γ_2  Γ_3}(|Γ_1|^2 + |Γ_2|^2)
                     + ω_a^2  |Γ_1|^2  |Γ_2|^2 
                     + |Γ_1|^4  |Γ_3|^2
                     + |Γ_2|^4  |Γ_3|^2 ] /[|Γ_1|^2  (κ^2 + (ω_c + ω_a)^2 - 4  |Γ_3|^2)
                      - |Γ_2|^2  (κ^2 + (ω_c - ω_a)^2 - 4  |Γ_3|^2) 
                      - 8  ω_c  {Γ_1  Γ_2  Γ_3}]
                     [( κ^2 + ω_c^2)  (ω_a^2 - 4  |Γ_3|^2)
                       - 2  ω_c  ω_a  (|Γ_1|^2 + |Γ_2|^2)
                       + (|Γ_1|^2  - |Γ_2|^2)^2                   
                       + 8  ω_c  {Γ_1  Γ_2  Γ_3}]   ,


which in the normal phase simplifies to Eq. (8) in the main text.

