Resonant light enhances phase coherence in a cavity QED simulator of fermionic superfluidity

Cavity QED experiments are natural hosts for non-equilibrium phases of matter supported by photon-mediated interactions. In this work, we consider a cavity QED simulation of the BCS model of superfluidity, by studying regimes where the cavity photons act as dynamical degrees of freedom instead of mere mediators of the interaction via virtual processes. We find an enhancement of long time coherence following a quench whenever the cavity frequency is tuned into resonance with the atoms. We discuss how this is equivalent to enhancement of non-equilibrium superfluidity and highlight similarities to an analogous phenomena recently studied in solid state quantum optics. We also discuss the conditions for observing this enhanced resonant pairing in experiments by including the effect of photon losses and inhomogeneous coupling in our analysis.

Superconductivity and superfluidity are among the most celebrated predictions of modern condensed matter theory, both for their fundamental importance and for the promise they hold to revolutionize power transmission [1,2]. Recent theory and experimental efforts point at potential non-equilibrium enhancement of superconducting-like phenomena in platforms at the interface of condensed matter and quantum optics, hinting at novel avenues beyond conventional high-temperature superconductors in solid state systems [3][4][5]. These encompass pump and probe experiments in the solid state setting [6][7][8][9][10], as well as proposals to enhance superconducting order using driven photonic cavities coupled to quantum materials [11][12][13][14][15]. The complexity in modelling the physical principles behind these platforms result from the necessity to combine materials science together with an understanding of the role of driven photonic and/or phononic degrees of freedom in manyparticle physics [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33]. It would be therefore desirable to provide an emulator of superconductivity which, although it may simplify the degrees of freedom involved, could shed light on complementary mechanisms for nonequilibrium enhancement of superconducting order. This could then be used as a stepping stone towards richer and more intricate scenarios.
Such an emulator has been proposed in AMO physics for quantum simulation of archetypal s-wave superconductors (for charged particles) or s-wave superfluids (for neutral particles) [34,35] . In these works the dynamics of the superfluid phase coherence, directly related to the Meisner and Anderson-Higgs mechanisms in superconductors [2], can be studied by monitoring the dynamics of the atomic phase coherence. In the QED simulators considered so far, the cavity must be far detuned from atomic frequencies so that photonic degrees of freedom can be integrated out  and so an effective matteronly s-wave model of superconductivity is sufficient to FIG. 1. Top panel: A simple schematic of the distribution of atomic levels (see text) and the spectral response of the cavity with detuning ∆c and linewidth κ. Bottom panel: Time average of the phase coherence S + = 1 N i σ + i as a function of the cavity detuning ∆c/(χN ), for large disorder W/(χN ) = 8. In the adiabatic limit, ∆c/(χN ) → ∞, the simulator shows vanishing coherence, S + (t) → 0, while coherence is maximum when the photon is at resonance with the mean atomic transverse field ∆c ≈ 0 = 6χN (marked by a black dashed line). To make use of an integrability analysis we assume an ideal cavity with κ = 0 for most of the paper, and then confirm that a realistic cavity linewidth, κ/(g √ N ) = 1.4 × 10 −2 , does not significantly modify the resonance phenomenon. describe the dynamics. In such a limit the cavity only contains virtual photons, and their primary purpose is to mediate pairing interactions.
In this Letter, we investigate the effect of real photons on the phase coherence when the cavity detuning to the atomic transition is reduced. In this limit, the single channel s-wave BCS Hamiltonian is no longer an accurate description, and instead the atoms and cavity field simulates the two channel model of the BCS-BEC crossover [2,58]. In this model, the effect of reducing photon detuning on the dynamics are non-trivial because, on one hand, reducing the detuning yields a stronger mediated interaction strength, while on the other hand reducing the detuning leads to retarded photon dynamics where an instantaneous interaction is no longer valid. Here, we find that even when the change in interaction strength is accounted for, the retarded photon dynamics can maintain phase coherence better than the instantaneous interaction. This is demonstrated in Fig. 1, where we show that upon reducing the photon detuning phase coherence increases until resonance, below which the diabatic (small detuning) limit takes over and phase coherence is lost. While these results are mostly obtained by a classical integrability analysis [58][59][60][61][62][63][64][65], we also find via numerical simulation that the phenomenon is robust to the non-integrable effects caused by inhomogeneous couplings and photon loss which are typically present in realistic cavity QED settings.
Simulation of Superfluid Phase Coherence. We consider the simulation of the two-channel model for the BCS-BEC crossover observed in ultracold fermion experiments [2,58]. The model involves fermions (with creation operatorf † k,s with momentum vector k and spin s) that can form Cooper pairs on the BCS side of the crossover or, bind into diatomic bosonic molecules at zero center of mass momentum (with creation operatord † ) on the BEC side of the crossover. Neglecting finite momentum molecular bosons, the dynamics are characterized by the Hamiltonian: whered is the mean molecular field, ∆ c is the molecular binding energy, and g is the coupling strength between fermions and molecules. When the fermions condense into a superfluid on the BCS side of the crossover, they mostly form Cooper pairs [2] quantified by the complex pair amplitudes ρ k = f † k,↑f † −k,↓ . In this Letter, we focus on the dynamics of the superfluid s-wave phase coherence S + = 1 N k ρ k , which quantifies the phase coherence between Cooper pairs with different pairing wave vector k.
Similar to Ref. [35], the Cooper pairs can be simulated by a collection of two level atoms (described by Pauli operatorsσ + i andσ z i ) via the Anderson pseudospin mapping [58,61,62]: where each atom i simulates a pair of fermion momentum modes i → ±k i . The above Hamiltonian can then be simulated by a cavity QED system similar to the experiments described in references [39,40,42,66], in which the internal levels of 2N atoms are encoded in long lived electronic states, e.g the 1 S 0 -3 P 1 states of 88 Sr atoms. The atoms are trapped in an optical lattice and are allowed to interact with a single cavity mode (described by a photon annihilation operatorâ simulating the molecular field,â →d). Such a system is modeled by the Hamiltonian [35,39,40]: where ∆ c is the detuning of the cavity from the mean atomic frequency, 2g i is the single-photon Rabi frequency, and i is an inhomogeneous effective transverse field.
Simulation of H f by the cavity QED system occurs for homogenous light-matter coupling g i = g and for a probability distribution, p( i ), of the inhomogeneous field, i , that is designed to match the density of states for the fermion model. We choose the density of states as is a box distribution with mean x 0 and width α (see Fig. 1). Similar to Ref. [35], such a bimodal distribution is chosen to ensure the possibility of persistent oscillations of the phase coherence (see below) in the W = 0 limit. A possible band structure reproducing this density of states and the superfluidity that would occur in the traditional thermal equilibrium setting is discussed in Ref. [67]. At large detuning, ∆ c g √ N and ∆ c 0 +W/2, the cavity field mediates spin-exchange interactions and an effective spin model can be derived which maps into a one channel BCS model as discussed in Ref. [35]. In this limit, an adiabatic approximation [35,39,40] assumes the state of the light field is in instantaneous equilibrium such that â † eq (t) = − gN ∆c S + , where S + = 1 N k ρ k is both the atomic phase coherence and the simulated superfluid phase coherence. Thus, in the large detuning limit, the photon directly measures the phase coherence S + . Inserting â † eq (t) back into Eq. 2 and taking homogenous couplings, one finds a mediated interaction −χ ii σ + iσ − i with interaction strength χ = g 2 /∆ c and sign which favors effective Cooper pair formation at low temperatures and positive detuning, ∆ c . In this work we will study the dynamics when the photon detuning, ∆ c , is decreased and the adiabatic approximation is no longer valid. One complication to this limit is that when the photon detuning is decreased, the interaction strength χ increases. To isolate this effect we imagine that the experiment simultaneously increases the strength of the atomic energies as the photon detuning is decreased such that 0 /χN and W/χN are held constant.
The dynamical phases are then classified by the required number of collective DOF and the dynamics of the phase coherence S + . First, we consider the resulting collective modes for a quench starting from an initial state with all spins polarized in thex direction, σ x i = 1, and the cavity in the vacuum, a = 0. In the spinonly model, three phases are found [35,58] with at most M = 2. In contrast, we identify a fourth phase with M = 3 upon introducing the photon away from adiabatic elimination. The three phases in the adiabatic limit, ∆ c → ∞, are (for fixed χN and 0 > χN ): • Phase I (M = 0): At large disorder, all phase coherence is lost, and the simulated superfluid enters a normal state: S + (t) → 0.
• Phase II (M = 1): Transition to this phase occurs as disorder is reduced, and involves only one effective degree of freedom (M = 1). In this phase, the magnitude of the phase coherence, |S + (t)|, is constant at late time, and the collective mode corresponds to precession of the phase of S + : S + (t) → |S + | e iµt .
• Phase III (M = 2): This phase occurs at even smaller inhomogeneous atomic broadening, and has M = 2 DOF. The collective mode shows persistent oscillations in |S + (t)| as shown in the lower panel of Fig. 2.
In this adiabatic limit, the critical disorder strengths between phase I and II, and II and III, depends non-trivially on 0 , but are on the order of the interaction strength χN . At finite detuning, ∆ c , the photon becomes another DOF in the collective oscillations of these three phases, and to distinguish the phases of the full model we will write them with a "+1" superscript. The phases I +1 and II +1 show the same qualitative dynamics of S + (t) as the phases II and III respectively, while a new phase III +1 is defined by aperiodic oscillations of |S + (t)| and requires M = 3 collective DOF (two macroscopically coherent spins and a photon). At large but finite ∆ c , the new phase III +1 involves the photon performing fast oscillations around a eq (t), the slowly evolving equilibrium value given by adiabatic elimination (see Fig. 2 for an example). In this limit, the aperiodic contribution to the oscillations of S + becomes small smoothly as function of ∆ c , and thus in the large detuning limit phase III +1 approximates phase III. This limiting behavior is the same for phases I +1 and II +1 which, for large detuning, approximate phase I and II respectively. Upon reducing the detuning, a rich dynamical phase diagram emerges as shown in the upper panel of Fig 2. In the W = 0 limit, there is only phase III +1 , while, at finite W , the cavity field has a broad impact on the dynamical phase diagram. In the diabatic limit, ∆ c /χN 1, the dynamics are much more sensitive to the inhomogeneities due to an inability of the cavity to mediate an effective interaction, and the transition to phase I +1 occurs at much smaller disorder in comparison to the large detuning limit. We also find a region at large disorder, W > 4χN , where phase II +1 occurs when ∆ c ∼ 0 which suggests phase coherence can be enhanced by setting the detuning on resonance with the atoms that have atomic energies close to 0 .
Mechanism of resonant phase coherence. The enhancement of phase coherence is confirmed in Fig. 1, and we explain the formation of this resonance by first considering finite but large detuning, such that 1/∆ c is still the fastest timescale. In this limit the dynamics are in Phase I +1 and the enhancement of phase coherence is very weak at long times, but the following simple picture holds. First, on a timescale of 1/∆ c , the initial polarization of the spins drive the photon into an excited state oscillating around a non-zero a eq (t = 0) = −N g/∆ c . Then, on a timescale of 1/W > ∼ 1/∆ c , the spins mostly dephase and a eq (t → ∞) ≈ 0. Once the spins mostly depolarize to their steady state, the photon remains oscillating around a small equilibrium, a(t) ≈ Ae iµt . The Lax analysis (see [58] and [67]) yields expressions for the frequency and amplitude of these small oscillations which have a simple analytical form when ∆ c > W > 0 : µ = ∆ c and A = χN/g = √ N χN/∆ c . From the perspective of the matter, the photon is effectively an external drive that pumps a small fraction of the spins into a coherent steady state. In the frame of reference of the photon (the effective external drive), the dynamics of each spin is fully described by a constant magnetic field, h = ( i − µ)ẑ + gAx, and we can solve for the steady state as: Since µ ∼ ∆ c , this expression correctly predicts the loss of coherence, S + , in the adiabatic limit shown in Fig. 1. Further away from the adiabatic limit, the separation of time scales, 1/W > ∼ 1/∆ c , that yields the simple picture above is no longer valid. Regardless, the Lax analysis still produces the same expression, Eq. 3, for the phase coherence in phase I +1 but now with a different A and µ that must be numerically determined by solving for the roots of a Lax vector (see [58] and SM). Since µ gives the precession frequency of the photon, it is expected to be close to the detuning µ ≈ ∆ c and this is what we find numerically. Eq. 3 therefore predicts the atom at site i will be in resonance when ∆ c = i . The coherence is then maximally enhanced when most spins are driven close to resonance and occurs when the drive, ∆ c , is at the center of the band of atomic frequencies ∆ c ≈ 0 . This approximation is confirmed by the peak in coherence shown in Fig 1. Although Eq. 3 provides an intuitive picture, simliar to a single particle resonance, when the system is in phase I +1 , the relevant enhancement of coherence at the resonance happens in phase II +1 where the cavity field and atomic coherence must both be treated as dynamical variables. As shown in Fig. 2, their dynamics in this regime show coupled nonlinear oscillations [58].
Experimental Realization. In the experiments of Refs. [39,40] an optical lattice is used to trap Sr atoms,  [39,40]. The bottom panel is computed assuming inhomogeneous couplings and an initial state prepared by a coherent drive through the cavity. Even though inhomogeneities and cavity losses reduce the coherence, we can observe a signature of the resonance as a minimum [75] in the photon density. Note that in contrast to Fig. 2 and Fig. 1, the initial state depends on the initial number of photons driven into the cavity which scales as 1/∆c. Both figures were obtained by numerical simulation of the Lindblad equations of motion at mean-field [76]. In the top panel, the decay rate of S + (t) is constant with ∆c/χN and proportional to 1/κ, but appears to increase with ∆c/χN in the figure because the unit for time, 1/χN , decreases with ∆c when g is fixed.
featuring a long-lived electronic clock transition with atomic decay rate of γ. The optical lattice is placed inside a standing wave optical cavity with linewidth κ. While both γ and κ destroy phase coherence at long times, we find that the effect of resonant phase coherence is still observable on times O(1/κ) provided we operate at large collective cooperativity (N g 2 κγ). Given that for long-lived Sr atoms, κ γ we neglect atomic decay. Fig. 1 shows the dependence of |S + (t)| on ∆ c , and demonstrates that the resonant enhancement can be maintained even with cavity loss. Furthermore, Fig. 3 depicts how the dynamics in Fig. 2 simply features a slow decay for moderate κ.
The experiments in [39,40] also have inhomogeneous couplings g i = g cos(k 0 a 0 i) with k 0 a 0 = 3.7 due to an incommensurability between the optical lattice spacing, a 0 , and the cavity wavelength, 2π/k 0 . The inhomogeneous couplings will disrupt the effect discussed in this work if we start in a homogenous state, since the couplings will no longer excite the photon. However, as long as the initial state is generated by coherently driving the optical cavity, inhomogeneities do not play a detrimen-tal effect. In this case, the initial state involves all spins aligned with the inhomogeneities sgn(σ x i ) = sgn(g i ) such that cavity will be coherently pumped by the atoms. The resulting simulations show signature of resonant phase coherence as a minimum [75] of the time averaged photon density shown in Fig. 3. Note that both dissipation and inhomogeneities break Lax integrability.
Conclusion. Our work demonstrates that dynamical fluctuations of a mediating field can produce enhancement of phase coherence in cavity QED simulators of superconductivity and superfluidity. The generality of our result based on a resonance argument, would suggest a natural extension to a broad variety of platforms such as trapped ions or quantum optics in waveguides, both of which serve as tunable simulators of non-equilibrium quantum many body physics, employing mediating photons or phonons [77,78]. It also suggests a promising direction for cavity enhanced superconductivity in real materials. Such a possibility requires extending the phenomenon to charged superfluids in which the lightmatter couplings are structurally different from the atommolecule couplings of Eq. 1. Our results offer the possibility of studying novel regimes of enhanced cooperative light-matter, and hint that quantum many-body optics with active light and matter degrees of freedom has the potential to become a blossoming area of quantum simulation in the near future.

Lax analysis for quench dynamics
In the main text, we studied the dynamics of an initial state with σ x i = 1 and with the photon in an empty cavity state. Here, we consider a more general state used in [35] where the spins 1 . . . N point in thex cos(∆φ 0 /2) + y sin(∆φ 0 /2) direction, while the spins N + 1 . . . 2N point in thex cos(∆φ 0 /2) −ŷ sin(∆φ 0 /2) direction, and the photon is either in the vacuum, a(t = 0) = 0 or at equilibrium with matter, a(t = 0) = a eq (0) = − 1 ∆c i g i σ − i (t = 0). We now compute the lax vector for such a state: Since the initial state is uniform for spins 1 . . . N and for N + 1 . . . 2N , the sums over i splits into two sums which can be approximated by disorder average: where p( ) is the disorder distribution for transverse fields: i = /2 + x i for i ∈ (1 . . . N ) and i = − /2 + x for i ∈ (N + 1 . . . 2N ), and x i are drawn from a uniform distribution with zero mean and width W/2. The integral results in a logarithm whose branch cut must be chosen to match the continuum of poles that develop when taking the N → ∞ limit of the summation [63][64][65]. This results in: where the branch cut for ln(x) is Re(x) ∈ [0, −∞), and the branch cut for ArchTanh(x) is Re(x) ∈ [−1, −∞) and Re(x) ∈ [1, ∞). The square lax vector is therefore computed as: where s = 0 if the cavity field is a(t) = 0 at t = 0, or s = 1 if the cavity field is in equilibrium with matter a(t = 0) = a eq (t = 0) at t = 0.
In the main text, we used the numerical search described below to study the dynamical phases as a function of ∆ c /χN and W/χN , for fixed 0 /χN = 6, ∆φ 0 = 0 and photon starting in the vacuum. Such a phase diagram does not qualitatively change by increasing ∆φ, but if the cavity field is in equilibrium with matter a(t = 0) = a eq (t = 0) at t = 0, and ∆φ 0 = π we see larger region for phase II +1 (see Fig 4). FIG. 4. Phase diagram as in main text but for an initial state with ∆φ0 = π and a(0) = aeq(0). Interestingly, for large W , the region of Phase II +1 shifts to higher frequencies. While this region is in phase II +1 where the amplitude S + oscillates, the oscillations and over all phase coherence are relatively small at large W . Similarly, the Lax analysis of phase I +1 still predicts a resonance in phase coherence S + at ∆c = 0 (marked by a black dashed line); we have confirmed this numerically.

Numerical Search for the roots
The dynamical phases are characterized by the number of roots of L 2 (u) = 0 as discussed in the main text. To find these roots for a range of W/χN and ∆ c /χN as in the main text, we employ the NLsolve Julia library. The algorithm requires an initial seed for the roots and then uses the steepest decent to find the numerical value. Therefore, to find all the unique roots L 2 (u) = 0 and thereby correctly identify the dynamical phase we perform the following process. First, we plot the complex magnitude of L 2 (u) and identify an approximate location of the 6 roots from the plot. Then we use NLsolve function to find a precise location of the root with a relative and absolute numerical error tolerance of 10 −4 . For the first root, this is done for W = 0, a specific initial state, a large value of ∆ c /χN , and the remaining Hamiltonian parameters fixed.
Since the location of the roots are smoothly connected to each other as ∆ c /χN is smoothly decreased, we next find the location of roots for W/χN = 0 and different values of ∆ c /χN . This is done by starting with roots found for W = 0 and the large value of ∆ c /χN and then finding the roots sequentially decreasing ∆ c /χN . The seeds used for the NLsolve function are the roots obtained from the previous step in the sequence with slightly larger ∆ c /χN . Then, starting with the root found for a specific ∆ c /χN and W/χN = 0, we find the roots for W/χN > 0 performing a similar sequential process, but this time increasing W/χN and keeping ∆ c /χN fixed. After finding all the roots, we identify the unique ones, up to an error tolerance of 10 −4 , to identify which phase we are in.
A problem that can arise in this process is two roots closely approach each other, as W/χN is increased. If this occurs, the seeds used for steepest decent can become identical to each other in the sequential update of seeds based on the previous roots. In order to prevent such issue, we keep a list of fixed seeds that we also use when finding distinct roots at each point in the phase diagram.
Analytic solution for phase I +1 roots As discussed in the main text, the location of the two roots in phase I +1 gives useful information on the dynamics of the photon. These roots can be found if when the cavity is initially empty, and when ∆φ 0 = 0. In this limit, the ArchTanh can be expanded, and the Lax squared vector yields: Defining u = u ∆ c /2, the condition L 2 (u) = 0 becomes: The quadratic equation solves as When ∆ c = ∆ c /χ/N is very large the small angle approximation yields u ≈ 1 ± i 2 ∆ c , or u = ∆c 2 ± iχN . The Lax analysis argues [58] that in Phase I +1 , the real and imaginary parts of this root pair gives the frequency, µ, and amplitude, A, of the photon oscillations respectively: a ≈ Ae iµt . Therefore we find that when, ∆ c > W > 0 the photon frequency is given as µ = ∆ c the photon amplitude as A = χN/g = √ N χN/∆ c . as in the main text. This is confirmed in Fig. 5.

Minimum in photon occupation at resonance
As seen in Fig. 5, and in Fig. 3 of the main text, the photon density 1 T N T 0 dt a † a shows a minimum as a function of the detuning, ∆ c , near resonance ∆ c = 0 . This phenomenon is explained in a similar way to the resonance maximum seen by the phase coherence |S + | and it is related to the conservation of total excitations: S z + a † a. In Fig. 5, the initial state is polarized in thex direction and with no photons in the cavity such that S z + a † a = 0. Therefore the photon density is fixed to the J z polarization of the spins: a † a = −S z . Using the same argument in the main text that, in phase I +1 , the photon simply acts as an external drive with a = Ae iµt , we find that Therefore the i th spin contributes least to the total S z polarization (and photon amplitude) when it is in resonance with the cavity field i = µ ≈ ∆ c . This implies a minimum in a † a when the spins on average are in resonance with the cavity field: 0 = ∆ c . Furthermore, a spin i with frequency i > µ will have oppositeẑ polarization than a spin j with frequency j < µ resulting in overall reduction of their net contribution to the total spin polarization S z . This deconstructive interference between spins with different relative frequency i − µ will shift the minimum in a † a to a smaller detuning than ∆ c = 0 due to the contributions of spins in the negative band of frequencies: The argument is similar for an initial state discussed in Fig. 3 of the main text, which also has S z (t = 0) = 0 but with photons in the cavity, |a(t = 0)| = |a 0 | > 0. The only difference is the conservation condition gives an additional contribution to the photon density: Band Structure We consider a dispersion, k,b , that has two bands (indexed by b = ±1) centered at b 0 /2 with bandwidths max k |k| − b 0 /2 = W/2, and a constant density of states within the two bands (see Fig. 6). Such a constant density of states can occur with Dirac cones and is required to be easily simulated by a uniform distribution for the inhomogeneous broadining on the atoms (see Fig. 6). Furthermore, we work in the limit W < 2 0 such that there are two bands separated by a band gap of 0 − W/2. Note that, while the energy difference between the atomic excited state and its ground state is 2 i , each atom represents two fermion modes with momentum −k and k such that i → k,b is the energy of a fermion at either of those momenta.

Equilibrium Phase Diagram
While in the main text we consider quench dynamics from a state in which every possible Cooper pair is condensed, intuition about superconducting systems is often developed near thermal equilibrium. To make connection to that limit, we consider the superconducting order of the grand canonical state ρ µ,T = exp[−(H − 2µ(S z + a † a))/K b T ] with chemical potential µ and temperature T → 0, and FIG. 6. An example band structure with density of states matching the distribution of the inhomogeneous atomic energies shown in the left panel. In this diagram, 0 > W/2 such that there is a band gap of size 0 − W/2. Note that center of the bands are at ± 0/2 for the fermion band structure, and not ± 0 as for the atomic level distribution. This is because there are two fermions per cooper pair, and it is cooper pairs that map to atomic levels.
This model has a U (1) symmetry generated by the conservation of total excitations i 1 2 σ z i + n. At high T and µ = 0, the system enters a symmetric phase with σ → 0 and with photon occupation scaling monotonically with K b T . At low temperature, a symmetric phase with all spins polarized on theẑ axis, σ z i = sign( i − µ) (a normal state insulator or conductor), competes with a symmetry broken phase can where | S + | > 0 and a > 0. Away from the quantum critical point separating these two phases, a mean field approximation (i.e.
is valid, and we again adopt the notation: O 1 ≡ O 1 . At zero temperature, energy is minimized at the fixed points of the classical dynamics, ∂ t σ = 0 and ∂ t a = 0. The later condition fixes the photon in a way mathematically equivalent to the adiabatic condition: Defining the gap as ∆ = g 2 ∆c−2µ i σ + i and inserting Eq. 14 into Eq. 13, one finds that the spins minimize energy in an effective magnetic field h = ( i − µ)ẑ + Re(∆)x + Im(∆)x. Minimizing the energy of the spins in this field we find: using g 2 = χN ∆ c we get the gap equation: where the integral is over the region defining the constant density of states: In addition to the gap equation, we can also restrict the total number of excitations to match that of the initial state in the paper: In the adiabatic limit, ∆ c /(χN ) 1, the photon amplitude goes to zero, and thus S z must also limit to zero. This is ensured by a chemical potential µ = 0 such that i − µ is positive for half the spins and negative for the other half. The resulting zero temperature phase diagram is calculated numerically and is shown in Fig. 7. At finite temperature, the superconducting coherence |S + | is reduced in the ordered phase and remains 0 in the disordered phase. Weather the disordered phase is a insulator or conductor at low temperature depends on if the chemical potential µ = 0 is located in a band or band gap. Upon consideration of Fig. 6, we see the disordered phase is a two-band insulator when 0 > W/2 and a single band conductor for larger dispersion W . Upon inspection of Fig. 7, we see that for the adiabatic limit studied in the paper, the system would traditionally be considered a two-band insulator.
The phase diagram shown in the left panel of Fig. 7 holds for all values of ∆ c /χN assuming the chemical potential is fixed to µ = 0. If, instead, µ is chosen to fix the total number of excitations, it will depend on the detuning of the photon since, in the ordered phase, the photon occupation increases monotonically with χN/∆ c . Therefore, to maintain a fixed number of total excitations, the chemical potential must increase to ensure a finite and negative FIG. 7. Left Panel : Equilibrium phase diagram at T = 0 and with chemical potential fixed to µ = 0. In the paper we study the line along 0/χN = 6 which has S + = 0 in the ground state and at all finite temperatures. Right Panel : Zero temperature coherence, S + , versus photon detuning, ∆c/χN , with µ chosen to fix the total number of excitations, a † a + 1 2 i σ z i = 0. This panel is for, 0/χN = 6 and W/χN = 8 as in the paper and shows a similar resonance peak.
S z = −a † a. This can lead to instabilities and produce an ordered phase where there was none in the adiabatic limit. Fixing a † a + 1 2 i σ z i = 0, we find a similar resonant order at zero temperature as was seen out-of-equilibrium (see right panel of Fig. 7). Note, such an equilibrium state is unlikely to be observed in the cavity QED system due to the effects of dissipation which lead to a disordered state with no photons in the cavity at long times.

Lindblad Dynamics
In order to investigate the effect of photon loss, we study dynamics of the spin and cavity evolving under the Liouvillian with jump operator L = κ/2a. In the Heisenberg picture operators evolve following where κ is the loss rate. In the mean field limit (i.e. O 1 (t)O 2 (t) = O 1 (t) O 2 (t) ), we can truncate the hierarchy of equations generated by Eq. 17. We then numerically evolve the closed set of equations for the spin and photon variables.