Emerging dissipative phases in a superradiant quantum gas with tunable decay

Exposing a many-body system to external drives and losses can transform the nature of its phases and opens perspectives for engineering new properties of matter. How such characteristics are related to the underlying microscopic processes of the driven and dissipative system is a fundamental question. Here we address this point in a quantum gas that is strongly coupled to a lossy optical cavity mode using two independent Raman drives, which act on the spin and motional degrees of freedom of the atoms. This setting allows us to control the competition between coherent dynamics and dissipation by adjusting the imbalance between the drives. For strong enough coupling, the transition to a superradiant phase occurs, as is the case for a closed system. Yet, by imbalancing the drives we can enter a dissipation-stabilized normal phase and a region of multistability. Measuring the properties of excitations on top of the out-of-equilibrium phases reveals the microscopic elementary processes in the open system. Our findings provide prospects for studying squeezing in non-Hermitian systems, quantum jumps in superradiance, and dynamical spin-orbit coupling in a dissipative setting.


I. INTRODUCTION
Open many-body systems can annul fundamental laws that typically govern the physics of systems in thermal equilibrium. In the idealized case of an ensemble of interacting particles, isolated from the environment and at zero temperature, the ground state is set by energy minimization and phase transitions arise from competing energy contributions [1,2]. In open systems however, the interplay between coherent dynamics within the system and its interaction with the environment gives rise to a much richer set of phenomena [3][4][5][6][7]. Such interaction is not only unavoidable, but can be exploited via the engineering of external drives and coupling to specific baths [8][9][10][11]. The experimental access to manybody ground states provided by ultra-cold atoms [12][13][14] laid the foundation for a recent revival of interest in many-body systems interacting with their environment. Experimental observations that are specifically due to the system's openness include the emergence of bistability [15][16][17][18], the stabilization of insulating phases [19,20], the access to absorbing-states phase transitions [21], the appearance of dissipation-induced instabilities [22] and time crystals [23], and the change in correlation properties [24,25] that can signal non-Hermitian phase transitions [26].
Besides their fundamental interest, non-equilibrium phenomena bear the prospect of becoming a powerful tool for engineering new materials ranging from exciton condensates to light-induced superconductors [27][28][29][30][31]. The properties of these phases of matter emerge from tuning the elementary excitations by hybridization with the light field [32,33], which provides a natural coupling to the external environment in presence of optical drives and losses [34,35]. To gain further insight into this phenomenology, it is desirable to achieve good experimental control over coherent and dissipative channels and at the same time to gain access to the microscopic properties lying at the origin of the macroscopic phases [36].
In this work, we engineer a driven-dissipative manybody system with global-range interactions, that is subject to tunable coherent and dissipative channels. Our realization employs a quantum gas strongly coupled to an optical cavity [37,38]. Building on schemes that have been extensively exploited with thermal atoms where the atomic spin is coupled to light fields [39][40][41][42][43][44], our implementation also involves the density degree of freedom of the gas [45][46][47]. We employ two Raman laser drives to control the strength of the co-and counter-rotating terms of the resulting light-matter coupling independently. In combination with photon loss from the cavity, this allows us to explore different regimes of competing coherent coupling and dissipation. Schematically, one can identify the following processes ( Fig. 1): the combined action of the two drives coherently mixes two many-body atomic states (|0 a , |1 a ) with the cavity field into polariton modes. As the strength of the drives increases, the excited polariton mode |1 softens to the lowest-energy one |0 , and a second-order phase transition occurs from a normal to a superradiant phase that is phase-locked to the drives. In addition to this coherent process, each individual drive can induce transitions from one atomic state to the other (|0 a |1 a ), where photons are scattered into the cavity mode and successively lost. As a result, adjusting the relative strength of the two lasers yields to tuning the effective polariton dissipation. This leads to qualitative changes in the phase diagram of the system, with the appearance of a dissipation-stabilized phase and a discontinuous phase transition in a multistable region. We relate these phenomena to the changing properties of the polaritonic excitations, which we characterize experimentally and theoretically.
FIG. 1. Competing coherent coupling and dissipation at a superradiant phase transition. (a) A quantum gas interacts coherently with a cavity mode via two tunable drives with mean coupling strengthη and imbalance ∆η, giving rise to two low-energy polariton modes |0 and |1 , corresponding to decoupled and coupled light-matter modes, respectively. Increasingη softens the energy of |1 (black line); cavity dissipation is responsible for the effective damping (γ ↓ ) and amplification (γ ↑ ) of the soft-mode polariton |1 . For small enough ∆η/η, the rates γ ↓ , γ ↑ are balanced, and the mode softening to zero energy atη = ηc is accompanied by a phase transition from the normal phase populating only |0 (gray shade) to the superradiant phase (green shade), where |1 is occupied. By increasing ∆η/η, the dominating damping γ ↓ of the soft mode leads first to bistability (green-gray hashed region), and finally to the suppression of the superradiant transition. (b) Sketch of the experimental setup and (c) corresponding level scheme. A BEC inside a high-finesse cavity (with resonant frequency ωc and field decay rate κ) is illuminated transversally by two Raman lasers with coupling strengths η b(r) and frequencies ω b(r) . The BEC (|0a ) couples to a spatially-modulated state |1a of the neighbouring Zeeman sublevel, separated by an energy ωz = (ωz − 2ωrec), with Zeeman splitting ωz and recoil energy ωrec. In the superradiant phase, a coherent field at frequencyω (green) builds up in the cavity. The twophoton transitions driven by each pump and the cavity field atω are detuned from the bare atomic states |0a , |1a by ∓ω0, as indicated by the lower dashed lines. The dissipative channels between modes |0 and |1 shown in (a) are due to Raman scattering of photons from each single drive into the cavity (blue and red wiggly arrows), and subsequent photon loss. The Lorentzian density of states of the cavity is sketched in orange.

II. SETUP AND TUNABLE DECAY
In our experiments, we trap a Bose-Einstein condensate (BEC) inside a high-finesse optical cavity and pump the atoms with a two-frequency optical standing wave, perpendicular to the cavity axis [ Fig. 1(b)]. The BEC is formed by N = 10 5 atoms of 87 Rb, prepared in the m F = −1 sublevel of the F = 1 ground state hyperfine manifold. A magnetic field along the z-direction is applied to generate a Zeeman splitting ω z = 2π · 48 MHz between the initial state and the m F = 0 sublevel. The driving fields are far red-detuned from the atomic resonance with frequencies ω b , ω r that lie on opposite sides of the dispersively shifted cavity resonance ω c , i.e., ω r < ω c < ω b , and ω b − ω r ∼ 2ω z . The standing-wave modulations of the two drives overlap at the position of the atoms, forming a one-dimensional lattice potential with spacing λ/2 = 784.7/2 nm. Each pump beam realizes a cavity-assisted Raman coupling between the m F = −1 and m F = 0 levels, as sketched in Fig. 1(c). The resulting system is effectively described using two atomic modes: the initial ground state of the BEC |0 a , and the excited-momentum state of the neighbouring Zeeman sublevel |1 a ∝ cos(k recx ) cos(k recẑ )F + |0 a , witĥ F + being the raising spin-operator in the F = 1 manifold, and k rec = 2π /λ the recoil momentum. When increasing the driving strength, the ground state |0 a evolves from a harmonically confined BEC to a stack of pancake-shaped BECs trapped in the maxima of the standing-wave drives.
In a rotating frame defined by the driving frequencies, the coherent dynamics of our system is described by the many-body Hamiltonian as detailed in [48]. Here,Ĵ i=x,y,z are the components of the pseudo-spin operator associated with the two manybody states |0 a ( Ĵ z = −N/2) and |1 a ( Ĵ z = N/2), andâ (â † ) is the annihilation (creation) operator of the relevant cavity mode. The effective atomic frequency is ω 0 = (ω b − ω r )/2 − ω z + 2ω rec , with the recoil energy ω rec , and the cavity detuning ∆ c =ω − ω c , wherē ω = (ω b + ω r )/2 is the mean of the pump frequencies.
The light-matter coupling strengths are parametrized byη = (η b + η r )/2 and ∆η = (η b − η r )/2, with the cavity-assisted Raman coupling η b(r) arising from the ω b(r) pump and the cavity mode. These two-photon Raman couplings implement the independently tunable co-and counter-rotating terms of the light-matter interaction. The dynamics of the open system due to photon losses is well described by a master equatioṅ ρ = −i/ [Ĥ,ρ] +L[ρ], where the Lindblad superoper-atorL[ρ] = κ 2âρâ † − {â †â ,ρ} accounts for the cavity field decay at rate κ = 2π · 1.25 MHz. The model introduced here is a generalized Dicke model that is predicted to exhibit rich phenomenology, see Ref. [49] and the more recent Refs. [50][51][52]. Correspondingly, first ex-periments exploring effective Dicke models with tunable co-and counter-rotating terms and different beam geometries have been carried out using thermal atoms [41,53], whose motional state is not well defined.
When the co-and counter-rotating coupling terms are balanced (∆η = 0), Eq. (1) reduces to the Dicke Hamiltonian [54,55]. In this limit, at low couplingη, the system is in the normal phase with Ĵ x = 0, â = 0, and the lowest-energy polariton mode |0 is mostly occupied. By increasing the couplingη, the energy of the first excited polariton |1 , admixing both atomic states |0 a , |1 a and the cavity photons, softens. As the energy of the polariton |1 reaches zero, the system undergoes a second-order phase transition to the self-organized superradiant phase with Ĵ x = 0, â = 0 [ Fig. 1(a)] [38,45,56]. The transition occurs at a collectively-enhanced critical coupling , which is only slightly shifted from the closed-system critical point for our parameters [48]. On the other hand, if an imbalance ∆η between the coupling of co-and counter-rotating terms is introduced, the effect of dissipation on the system becomes qualitatively different. Specifically, due to cavity loss, each Raman drive η b(r) is responsible for an effective decay γ ↓(↑) of the polariton mode |1 (|0 ) towards mode |0 (|1 ). In the parameter regime of our experiment ω 0 κ, the effective decay rates take the form which we derive using an effective Keldysh action for the polariton modes [48]. We identify that the microscopic process corresponding to the effective decay γ ↓(↑) is a collectively-enhanced Raman scattering driven by the ω b(r) pump beam, into the dissipatively-broadened density of states of the cavity, as sketched in Fig. 1(c). This mechanism is analogous to the Raman decay lying at the heart of superradiant Raman lasers [57][58][59][60]. Note that the effective decays (2) are independent of the phase of the cavity field. This is different from the coherent Raman couplings leading to the superradiant phase, where the intra-cavity field is always in-or outof-phase with the effective driving fieldω [61]. The two processes γ ↓(↑) act against each other, either damping or amplifying the population of mode |1 , such that global dissipative effects on the system are enhanced when the pumps are not balanced. In particular, as we experimentally demonstrate, the effective damping generated by these processes leads to a dramatic modification of the superradiant phase transition, as well as to new regions of multistability and hysteresis.

III. PHASE DIAGRAM
We restrict the experiments to the parameter space 0 ≤ ∆η ≤η; the properties of the system in the re-  (2), a coherent cavity field builds up in the cavity (superradiant phase) above a critical coupling (gray vertical line). (b) Green: intra-cavity photon number n ph as a function of the couplingsη, ∆η, obtained by implementing the protocol shown in (a) for 51 different values of ∆η/η, and averaging over 5 repetitions. Labels NP, SP and DSNP stand for normal, superradiant and dissipation-stabilized normal phase, respectively. The dots indicate the critical point extracted from the measured photon traces, and averaged within each subset of constant ∆η/η. The slope (∆η/η)DSNP of the dashed line corresponds to the smallest value of ∆η/η at which the SP is not present (the uncertainty on the slope is marked as a shaded gray region around the line). The black line is a guide to the eye through the critical points (obtained as a sliding average over 7 points) and along the measured boundary of the DSNP. We mark the boundary of the SP obtained from a mean-field stability analysis (solid, red), and from a numerical simulation of the measurement protocol (dashed, red) [48,62]. For comparison, the analytical mean-field result for the closed system is also plotted (blue, dotted). The arrow indicates the measurement path followed in (a). Inset: slope (∆η/η)DSNP extracted from phase diagrams measured at different cavity detunings ∆c, and plotted vs. gion 0 ≤η ≤ ∆η are mirrored, cf. Eq. (1). We map out the phase diagram of the system by ramping up the power of the pump beams while keeping the ratio ∆η/η constant, and monitoring the cavity output with a heterodyne detection, see Fig. 2(a). The onset of a superradiant phase is signalled by the build-up of a coher-ent cavity field with frequencyω above a critical coupling strength. We show in Fig. 2(b) the measured mean intra-cavity photon number n ph in the (η, ∆η) parameter space. At small imbalances ∆η η, the phenomenology of the Dicke phase transition is observed, and the value of the couplingη ≈ η c at which the transition occurs depends only weakly on ∆η. In contrast, at larger ratios ∆η/η > 0.71 (2), the superradiant phase transition is suppressed and the system remains in the normal phase at values ofη largely above η c . We compare the measured critical couplings with the phase boundaries obtained from a mean-field treatment of our drivendissipative model, and observe consistency between the experiment and our theoretical description. The phase boundaries are calculated from both the steady-state solutions and numerical simulations including time-varying couplings, with no free parameters (see [48] for details).
The existence of a dissipation-stabilized normal phase near the line ∆η/η = 1 is a consequence of the open character of our system, as pointed out in previous theoretical works [50,62]. To further characterize this phase, we measured the full phase diagram of the system for different cavity detunings ∆ c [48]. We observe that the width of the dissipation-stabilized normal phase increases for smaller detunings [inset of Fig. 2(c)]. This agrees with the predictions of our theoretical model, by which we find that the slope (∆η/η) DSNP of the boundary between the dissipation-stabilized normal phase and the superradiant phase at large couplingsη η c is given by (∆η/η)

IV. PROBING THE EXCITATION SPECTRUM
The effective damping induced by the relative coupling imbalance ∆η/η leads to a profound change in the superradiant phase transition and even suppresses it. This observation is closely linked to the system's spectrum of collective excitations. We implement an experimental protocol that allows for non-destructive, real-time monitoring of the free evolution of the excited polariton mode |1 , both in amplitude and in frequency. We promote a small population of mode |0 to mode |1 by means of a Bragg scattering process involving the transverse pump beams and a 1 ms-long laser pulse injected into the cavity at frequencyω + ω 0 + δ probe . This small occupancy of the excited mode produces a scattering of a weak field from the pumps into the cavity also after the pulse ended, which we monitor with the frequency-resolving heterodyne detector [63]. We shine the excitation pulse at a fixed time while ramping up both couplings at constant ratios ∆η/η. At the falling edge of the pulse, the mean pump intensities correspond toη/η c ≈ 0.6 [ Fig. 3(a)].
In Figs. 3(b-e), we present exemplary spectrograms of the cavity field showing the excitation pulse and the subsequent evolution of the system for increasingly larger values of ∆η/η. In the Dicke limit [∆η/η = 0, Fig. 3 Time (ms) the main components of the spectrum evolve towardsω as the couplingη is swept to larger values, reflecting the softening of the excited mode energy. At the critical point η = η c , the energy gap between mode |0 and the soft mode |1 vanishes, and the superradiant phase transition occurs, signalled by the build-up of a strong coherent field at frequencyω. As the relative imbalance ∆η/η is increased, the mode softening is accompanied by a faster decay of the excitation amplitude during the experiment. The 1/e-lifetime τ of the free excitation extracted from the integrated spectrograms decreases rapidly as the relative imbalance ∆η/η is increased [ Fig. 3(f)]. At the same time, the superradiant phase transition occurs at a couplingη that depends only weakly on the coupling imbalance ∆η until, for large enough ∆η/η, the transition is fully suppressed. The connection between the damping of the excitations and the suppression of the superradiant phase transition can be understood by analyzing the excitation spectrum of the open system. We linearize the mean-field equations of motion around the normal phase and study the low-energy eigenfrequencies ω ± of the corresponding dynamical matrix, as a function of the couplingη/η c , and for different values of ∆η/η [ Fig. 3(g,h)]. The real part of the spectrum captures the energy of the excited polariton, while a negative (positive) imaginary part signals damping (amplification). At first order in ω 0 /κ 1, the eigenfrequencies are given by as illustrated in [48]. As the couplingη increases towards the critical point η c , the phase transition occurs whenever any of Im[ω ± ] becomes positive [55,64,65]. At large enough ratios ∆η/η, the damping rate γ ↓ of the soft mode is dominant; this counteracts the coherent build-up of superradiance and stabilizes the normal phase. We report in Fig. 3(f) (blue shaded region) an estimate of the quasi-particle lifetime τ = − (2Im[ω ± ]) −1 obtained from the spectrum of the eigenvalues. For this estimation we assume that the couplingsη, ∆η are kept fixed at the end of the excitation pulse; this simplification provides a valid approximation for the measured lifetime where τ varies only slightly during the decay, i.e., as long as ∆η/η 0.05. For a closer comparison with the experimental data, we perform a numerical simulation of our experimental protocol including the excitation pulse [red shade in Fig. 3 To account for collisional interactions and spin dephasing, the theoretical estimations include a phenomenological atomic damping (see [48]), which we assume constant throughout the dynamics.

V. BISTABILITY AND HYSTERESIS
We focus now on the boundary between the superradiant phase and the dissipation-stabilized normal phase observed at large imbalances ∆η/η. Previous theoretical works predicted an intermediate region where both phases are stable (bistability) [50,62]. Accordingly, the occurrence of a first-order phase transition is expected. The bistability can be understood in terms of a competition between the coherent and dissipative processes described above. By increasing the relative imbalance ∆η/η, the damping γ ↓ fosters the population of mode |0 and acts against the coherent coupling responsible for superradiance. In the limit of dominant dissipation, such damping makes the dissipation-stabilized normal phase the unique stable steady state. Conversely, in the regime of comparable coherent coupling and dissipation, the phase to which the system converges depends on its initial preparation. If the system is prepared in the normal phase, it remains stable because of the aforementioned damping. On the other hand, if the system is initially in the superradiant phase, the dissipative damping towards mode |0 is counteracted by the presence of a coherent intra-cavity field that contributes to keep mode |1 significantly populated. In other words, preparing the system in the organized, symmetry-broken superradiant phase makes it more rigid against transitions assisted by cavity dissipation.
To explore the boundary between the dissipationstabilized normal phase and the superradiant phase, we initialize the system in the former phase by preparing the BEC and ramping up the couplingη above η c at fixed imbalance ∆η/η = 0.78. From this initial state, the transition to the superradiant phase is crossed by reducing ∆η at constantη. Then, within the same experimental realization, the direction of the ∆η sweep is reversed and the dissipation-stabilized normal phase is retrieved. The comparison between the forward and backward paths shows a hysteretic behavior at the phase boundary, in agreement with the expected discontinuous character of the transition, see Fig. 4. By performing the hysteresis measurement at different couplingη, an experimentally accessible region where the superradiant and the dissipation-stabilized normal phase are both stable is mapped out [ Fig. 4(c)]. We verified that implementing the parameter loop in the opposite direction also gives rise to hysteresis [48].

VI. CONCLUSION AND OUTLOOK
We showed that varying the imbalance of co-and counter-rotating coupling terms between a quantum gas and a lossy optical cavity engenders a tunable competition between coherent and dissipative processes across the superradiant phase transition, leading to the emergence of a dissipation-stabilized phase and hysteresis. Combining the control over dissipative and coherent couplings with real-time access to the dynamics allowed us to identify the underlying microscopic processes determining the observed phase diagram. We note that, if the cavity dynamics is adiabatically traced out, our system maps to a driven-dissipative version of the anisotropic Lipkin-Meshkov-Glick (LMG) model [66][67][68][69][70]-an important reference model for quantum magnetism that describes a many-body spin system with all-to-all interactions. Furthermore, interesting phenomena are expected near the boundary of the normal phase [52]; as visible in Figs. 3(g,h), here the real parts of the eigenfrequencies of the low-lying polariton merge, and their imaginary part bifurcates such that one squeezes, while the other broadens. In the Dicke limit, this phenomenon is limited to the coupling interval between the bifurcation and the transition point to the superradiant phase, which is very narrow for typical experimental parameters [55]. Squeezing of fluctuations can be obtained on a much wider range of parameters at the boundary between the normal phase and the dissipation-stabilized normal phase, where dissipation suppresses the change of the system's steady state. Characterizing the fluctuations of the normal modes in this regime will shine light on the generation of squeezing in the vicinity of exceptional points in non-Hermitian systems [71][72][73]. Moreover, performing experiments with small atom or photon numbers would unveil effects beyond mean-field, such as quantum jumps in the bistability region, as predicted recently in Refs. [74,75]. Furthermore, combining our findings on prominent dissipative effects with the generation of cavity-mediated spin-orbit interaction [76,77] opens a way to study spin-orbit coupling in a dissipative setting.

ACKNOWLEDGMENTS
We are grateful to M. Landini  We prepare a Bose-Einstein condensate (BEC) of 87 Rb atoms in the m F = −1 Zeeman sublevel in the F = 1 hyperfine manifold of the 5 2 S 1/2 electronic level, using radio-frequency assisted evaporation in a magnetic quadrupole trap. The atom cloud is optically transported and confined at the center of the cavity mode by an optical crossed dipole potential V ext , with frequencies [ω hx , ω hy , ω hz ] = 2π · [220(3), 24.6(8), 170.1(3)] Hz.
We apply a magnetic field B = B z e z , with B z < 0. To measure the Zeeman shift ω z between the sublevels m F = −1 (high-energy level) and m F = 0 (low-energy level), we employ cavity-assisted Raman transitions. We illuminate the BEC with the transverse pump with frequency ω r < ω c . Close to the two-photon resonance ω r ≈ ω c − ω z , a large fraction of the BEC is transferred to m F = 0 via two-photon processes involving absorption of photons from the red pump and emission into the cavity mode, with an increase of the kinetic energy by 2 ω rec . In this Raman process, light is scattered into the cavity at frequencyω = ω r + (ω z − 2ω rec ), fulfilling energy conservation. We infer the Zeeman splitting ω z by measuring the frequencyω of the photons leaking out of the cavity using a heterodyne detection scheme.

B. Transverse pump characterization
The two transverse pumps (drives) are derived from the same laser source. Their frequencies ω r and ω b are independently adjusted by employing double-pass acoustic optical modulators (AOMs) in different optical paths and recombining the beams afterwards. A small fraction of the beam is split close to the vacuum chamber and overlapped with an optical local oscillator at frequency ω LO , also derived from the same laser, on an AC photodiode. The beat notes at frequencies ω r − ω LO and ω b − ω LO are electronically separated and employed for intensity stabilization of the individual pumps. The distance between the retro-reflecting mirror and the atomic cloud is carefully adjusted such that the two standing waves overlap maximally at the position of the atomic cloud. The lattice depth associated to each pump is calibrated by means of Kapitza-Dirac diffraction [1]. For all the measurements presented in the main text, we increase the transverse pump powers via ramps of the form is the final lattice depth of the ω r,b pump and t ramp is the duration of the ramp.

C. Heterodyne detection
We monitor the photon field leaking out of the cavity by separating on a polarizing beam-splitter (PBS) the y− and z−polarization modes, and detecting each of them with heterodyne setups. The detection branch for the z−polarization is used to produce the data discussed in this work. The auxiliary detection setup for the y−polarized mode is used to probe the cavity resonance at the end of each experimental cycle.
In the heterodyne setup for the relevant z−polarized mode, the light field from the cavity is interfered with a local oscillator laser at frequency ω LO , and the high detection bandwidth of 250 MS/s allows for an all-digital demodulation of the beat note at ω − ω LO over a broad frequency range of [0, 125] MHz. In order to calibrate the photon number, we inject an on-resonance laser field into the empty cavity. We find the conversion factor between the demodulated heterodyne signal and the mean intra-cavity photon number by measuring the power after the PBS with a powermeter, and using the knowledge on the losses of the cavity mirrors.
The complex intra-cavity field α(t) = X(t) − iY (t) is obtained from the quadratures X(t) and Y (t) after digital demodulation at frequency δω D =ω − ω LO . The power spectral density (PSD) is calculated as PSD(f )=|FFT(α)| 2 (f ) using a fast Fourier-transform of the form FFT(α) [2], where t i is the time of the i th step and N is the total number of steps in the integration window. To construct the spectrograms, the traces are divided in time intervals of T = 150 µs with an overlap of 50% between subsequent intervals. Finally, we calculate the photon number spectrograms asñ ph = PSD(f )/T .

SII. DERIVATION OF THE HAMILTONIAN
In this Section, we derive the Hamiltonian in Eq. (1) from the closed-system dynamics of our spinor BEC coupled to the cavity.

A. Single-atom Hamiltonian
The Hamiltonian of a single atom coupled to the cavity mode readŝ (S1) In the dispersive regime of atom-light interaction [3,4], the optically excited states of the atom can be adiabatically eliminated, and the Hamiltonian of the bare atomĤ at can be written in terms of the ground states levels |F, m F only:Ĥ wherep and m are respectively the momentum and the mass of the atom, V ext (x) is the trapping potential, which is kept fixed, and ω F,m F is the energy of the |F, m F atomic level, with F denoting the hyperfine manifold, and m F the Zeeman sublevel. In our experiment, the 87 Rb atoms are initialized in |F = 1, m F = −1 , and near-resonant two-photon Raman transitions couple them to |F = 1, m F = 0 . Transitions to the F = 2 manifold are off resonance by the hyperfine splitting ω HF = 2π · 6.834 GHz, and can be neglected. The internal dynamics of each atom can then be described in terms of the spin operatorF = (F x ,F y ,F z ) T , with F = 1. The energy difference between the Zeeman sublevels is determined by first-and second-order Zeeman shifts ω z < 0 and ω z > 0, such that Eq. (S2) can be written asĤ The Hamiltonian of the bare cavity mode readsĤ where the operatorâ † creates z-polarized photons in the TEM 00 cavity mode with resonance frequency ω c . In the dispersive regime at which we operate, the interaction between the light fields and the atom takes the form whereÊ (−) (Ê (+) ) is the negative (positive) part of the total electric field at the position of the atom, withÊ (+) = (Ê (−) ) † , and α s , α v are respectively the scalar and vectorial polarizability at the frequency of the driving lasers [3][4][5], with α s < 0 and α v > 0. In Eq. (S5) we are neglecting an additional rank-2 tensor contribution to the polarizability, which is justified for 87 Rb at the wavelength λ = 784.7 nm at which we operate.
We consider classical y-polarized transverse pump fields propagating along the z-direction at frequencies ω r,b , and a quantized cavity field. The negative partÊ (−) of the total electric field is given bŷ with unit vectors e j (j ∈ {x, y, z}) and spatial mode profiles f r (x), f b (x), g(x). The two laser drives with amplitude E r , E b originate from the same optical fiber, and their standing-wave modulations overlap in phase at the position of the trapping potential. Given the small frequency difference ω b − ω r = 2π · 96 MHz, we can consider the same wavevector k =ω/c for the two drives interacting with the atoms, withω = (ω b +ω r )/2, and restrict to a single spatial profile f r (x) = f b (x) = f (x) = exp[−2x 2 /w 2 x − 2y 2 /w 2 y ] cos(kz). We also take g(x) = exp[−2(y 2 + z 2 )/w 2 c ] cos(kx) for the cavity mode profile. The waist sizes of the modes are [w x , w y , w c ] = [24,27,25] µm. The amplitude of the cavity field per photon is defined by the frequency and volume of the mode, resulting in E 0 = 403 V/m. We introduce the auxiliary HamiltonianĤ rot = ωâ †â − ω zFz , and perform the unitary transformationÛ = exp[ i Ĥ rot t], with ω z = (ω b − ω r )/2. By making use of the rotating wave approximation, we obtain the timeindependent HamiltonianĤ with cavity detuning ∆ c =ω − ω c and effective linear shift δ z = ω z + ω z . The interaction part has a scalar and a vectorial contributionĤ s ,Ĥ v given byĤ respectively. Note that we applied a global rotation of the cavity field of the formâ →âe iπ/2 . The first term in the scalar interactionĤ s [cf. Eq. (S10)] describes the attractive potential created by the transverse drives, giving rise to a one-dimensional lattice along the z-direction, with on-axis depth V TP = −α s (E 2 b + E 2 r )/4, and to an additional confinement along the x− and y−direction. The second term inĤ s is responsible for the dispersive shift of the cavity resonance and for a weak one-dimensional lattice potential along x-direction. We define the maximal dispersive frequency shift per atom as U 0 = α s E 2 0 / . The vectorial interactionĤ v in Eq. (S11) produces spin-changing Raman transition between the Zeeman sublevels of the F = 1 manifold. The spin-changing termsF x ,F y are mediated by orthogonal quadratures of the cavity field and can be tuned by the sum or difference of the two pump fields, respectively.

B. Many-body Hamiltonian
We derive the Hamiltonian for the many-body system of N atoms in a degenerate Bose gas using the secondquantization formalism. Using the single-atom results from the previous section, the many-body Hamiltonian can be written asĤ where theΨ(x) is the spinor atomic field operatorΨ(x) = Ψ +1 (x),Ψ 0 (x),Ψ −1 (x) T , fulfilling the bosonic commutation relations Ψ i (x),Ψ † j (x ) = δ ij δ(x − x ) and Ψ i (x),Ψ j (x ) = 0, with i, j = +1, 0, −1. At this level, we neglect collisional interactions assuming low densities.
We set the half-frequency difference ω z between the drives close to the energy separation between levels |m F = −1 and |m F = 0 , i.e., ω z ≈ ω z , with ω z = −ω (1) z + ω (2) z = 2π · 48 MHz. Thus, spin-changing Raman transitions to |m F = +1 are off-resonance by ∆ +1 ≈ 2ω (2) z = 2π · 0.7 MHz. The large detuning ∆ +1 determines the fastest timescale of the atomic evolution. This allows to adiabatically eliminate the atomic operatorΨ +1 , and restrict the dynamics to the two Zeeman sublevels with m F = 0, −1. The effect of the eliminated state leads to a sub-kHz energy shift of state |m F = 0 , which we neglect. We map our system to an effective generalized Dicke model by further restricting the Hilbert space to two spinmomentum modes, in the same spirit as previous works [6,7]. In the normal phase, the BEC, prepared in m F = −1, occupies the ground state of the total trapping potential, resulting from the combination of the external trap V ext and the attractive lattice potential V TP of the laser drives [cf. Eq (S10)]. We label this ground state as |0 a , with corresponding wave function Φ 0a (x). The cavity-mediated spin-changing interaction couples |0 a to a densitymodulated state |1 a in m F = 0, with wave function Φ 1a (x) = N Φ 0a (x) cos(kx) cos(kz), with N being a normalization factor. Within this two-mode description, the spinor field operator takes the formΨ( whereĉ 0a ,ĉ 1a are bosonic annihilation operators for the respective atomic modes. The corresponding expression for the many-body Hamiltonian iŝ We indicate with ω 0 (V TP ) the energy difference between the bare atomic modes. The quantities I(V TP ) and M(V TP ) are overlap integrals defined by I(V TP ) = 0 a | g(x) 2 |0 a /N and M(V TP ) = 0 a | f (x)g(x) |1 a /N . In writing Eq. (S13), we neglected the dependence of the dispersive cavity shift onĴ z , which is a valid approximation whenever the system is not deep into the superradiant phase.
By considering ω 0 (V TP ) as a constant, and taking its value in the limit of small V TP , i.e., ω 0 = (ω z − ω z + 2ω rec ), we can write the Hamiltonian in Eq. (S13) in terms of the parameters defined in the main text: withĴ ± =Ĵ x ± iĴ y . From Eq. (S15), it is apparent that the couplings η b , η r tune the strength of the co-and counterrotating term of the light-matter interaction, respectively. In the limit η r = 0, the Hamiltonian (S15) reduces to the Tavis-Cummings model [8].

C. Mapping between Hamiltonian couplings and experimental parameters
We describe the mapping between the measured power of the transverse pump beams and the Raman couplings η b , η r introduced in the main text. The BEC is trapped in the combined potential of the harmonic confinement V ext (x) and of the attractive potential created by the transverse pumps V TP (x) = −V TP f (x) 2 , which has contributions from the two drives, i.e. V TP = V b + V r . We monitor the power of each drive and extract the corresponding value of V b(r) in real time as described in the previous section.
To calculate the wave function Φ 0a , we consider spin-independent s-wave scattering and employ a Thomas-Fermi approximation for the interacting BEC in the total trapping potential [9]. Spin-changing collisions can be neglected due to the large second-order Zeeman shift ω (2) z = 2π · 0.35 MHz at which we operate [10]. In addition, we treat the one-dimensional lattice created by the transverse pump in the limit of large depth, and approximate the lattice as a succession of independent harmonic traps. This is justified by the fact that, in our experiments, the superradiant phase transition occurs at large V TP 25 ω rec . We then calculate the three-dimensional overlap integrals I(V TP ), M(V TP ) by using the expressions of the mode functions f (x) = exp[−2x 2 /w 2 x − 2y 2 /w 2 y ] cos(kz) and g(x) = exp[−2(y 2 + z 2 )/w 2 c ] cos(kx), with waist sizes [w x , w y , w c ] = [24,27,25] µm. The divergence of each Gaussian mode over the extension of the BEC is negligible. The Raman couplings η b , η r are then found where α v /α s = −0.928 and sgn[α s ] = −1 at the wavelength of the laser drives.
The overlap integral M(V TP ) converges to M max = 0.68 at large lattice depths. In the regime V TP 25 ω rec at which the superradiant phase transition occurs, M(V TP ) deviates from M max by less than 2%. For simplicity, we then assume M(V TP ) = M max when applying the conversion in Eq. (S16) throughout the paper.

A. Analytical calculation of the steady state
Starting from the Hamiltonian in Eq. (1) of the main text, we consider dissipation due to photons leaking out from the cavity at rate κ in the form of a Lindblad operator (S17) We disregard the spin decay rate due to the negligible spontaneous emission between Zeeman sublevels. At this level, we also neglect the role of spin dephasing, which would lead to a negligible shift of the phase boundaries for our experimental parameters. Using the master equation we obtain mean-field equations of motion (EOMs) of the form where the mean-field order parameters are â = √ N α, Ĵ x = N X, Ĵ y = N Y and Ĵ z = N Z. Introducing the renormalized couplingsη N = √ Nη , ∆η N = √ N ∆η, and imposing the spin constraint X 2 + Y 2 + Z 2 = 1 4 , we can solve analytically for α, X, Y, Z in the steady state. The normal phase corresponds to the trivial steady-state solution α Re = α Im = X = Y = 0, Z = −1/2, with α Re = Re[α] and α Im = Im[α]. The non trivial solutions of Eq. (S19) [11] read (S21)

B. Analytic expressions for the phase boundaries
We find an analytic expression for the slope of the dissipation-stabilized normal phase starting from the expression of Z from Eq. (S20): Requiring Z to be real results in the condition The equality provides the slope of the boundary between the superradiant phase and the dissipation-stabilized normal phase, cf. Fig. 2 in the main text: for ∆ c < 0. Moreover, from Eq. (S22) we also obtain the stability boundary of the normal phase. Specifically, in the normal phase we set Z = −1/2 on the left hand side of Eq. (S22) and square both sides, we then solve for ∆η N obtaining the following expression which allows us to find the boundaries of the bistability region shown in Figs. 4(b,c) in the main text.

C. Numerical simulations of the mean-field dynamics
In order to simulate the time evolution of the system, we numerically solve the semi-classical EOMs (S19). For this purpose, we use the MATLAB built-in 'ode45' solver which is based on a Runge-Kutta (4,5) method [12]. It employs variable time step sizes and the error tolerance in each step is constrained to 10 −8 . To sample the fluctuations on top of the mean-field observables and allow for a phase transition to take place, we assume an initially small photon field of the form α(t = 0) = [randn(0, 0.5) + i · randn(0, 0.5)]/ √ N with pseudo-random numbers randn(0, 0.5) sampled from a normal distribution with (µ, σ)=(0,0.5). This assumption is compatible with an initial coherent vacuum state for the cavity field since We plot in Fig. S1 the phase diagrams obtained from the numerical simulations and analytic steady-state calculations for the experimental parameters of Fig. 2 in the main text. The red lines in each plot indicate the boundaries of the superradiant phase assuming a threshold photon number of n ph,th = 5. We attribute the small shift of the phase boundaries of the numerical simulations from the analytical results to residual non-adiabatic effects.
D. Measurement of the phase diagram

Experimental protocol and data processing
To measure the phase diagram, we prepare a BEC in m F = −1 and ramp up the power of the driving lasers at constant cavity detuning ∆ c . The calibrated heterodyne signal is used to construct photon number spectrogramsñ ph (f, t) as described in the next section. We integrate them in a narrow frequency range of P = [ω/2π−2.5 kHz,ω/2π+2.5 kHz] to obtain the photon traces n ph (t) = f ∈Pñ ph (f, t) of the cavity field at the frequencyω characteristic of the superradiant phase.
The phase diagram in Fig. 2(b) of the main text is obtained by combining measurements for 51 different ratios ∆η/η, with 5 realizations each. From each realization, we extract the coupling rampsη(t), ∆η(t) by monitoring in real time the power of the driving lasers, as well as the time trace of the mean cavity photons n ph (t) [cf. Fig. 2(a) in the main text]. After parametrizing n ph as a function ofη and ∆η, the colorplot in Fig. 2(b) of the main text is obtained by averaging the different experimental realizations in the (η, ∆η) parameter space within squared bins with size 2π · 24 Hz.
From each experimental realization, we extract the time t th at which the superradiant phase transition occurs by fitting n ph (t) with a piecewise linear and power law function. The critical couplings are obtained as (η th , ∆η th ) = (η(t th ), ∆η(t th )). The dots in Fig. 2(b) are the average of the critical couplings (η th , ∆η th ) from measurements taken with the same ratio ∆η/η. The errorbars are the corresponding standard error of the mean. The experimental boundary between superradiant phase and the dissipation-stabilized normal phase [dashed line in Fig. 2(b) in the main text] corresponds to the smallest value of ∆η/η at which the phase transition is not observed in at least one of the experimental realizations, which we define as (∆η/η) DSNP . The upper (lower) boundary of the gray shaded region around this line marks the smallest (largest) ratio ∆η/η at which the phase transition is absent (present) in all realizations. These bounds provide an uncertainty to the slope (∆η/η) DSNP .

Phase diagrams for different cavity detunings
We record experimental phase diagrams for three different cavity detunings ∆ c and display the results in Fig. S2. From each phase diagram, we extract the slope of the phase boundary (∆η/η) DSNP and its uncertainty as discussed above. The results of (∆η/η) DSNP vs. −κ/∆ c are shown in the inset of Fig. 2(b) in the main text.

Absence of additional self-organization processes
The drive at frequency ω r , red-detuned from cavity resonance, may induce spurious atomic self-organization that does not involve a change of the m F state [6]. If present, this process would be accompanied by the build-up of a coherent cavity field with the same polarization (along y) and frequency (ω r ) as the driving laser. During our measurements, we continuously monitor such cavity field with the auxiliary heterodyne setup, and never observe any signal above noise. The suppression of spin-preserving self-organization is due to the large detuning |ω r − ω c | and to the presence of the second drive at frequency ω b .

A. Spectrum of the open system
To find the spectrum of the open system, we consider the system's EOMs (S19) and expand the order parameters as α = α 0 + δα, X = X 0 + δX, Y = Y 0 + δY, Z = Z 0 + δZ around the steady states (α 0 , X 0 , Y 0 , Z 0 ), chosen to be either the normal or the superradiant phase. A similar treatment has been used in Refs. [11,13]. The linearized EOMs read d dt with where for completeness we phenomenologically included atomic dephasing at rate Γ φ . This damping term is compatible with a Lindblad term of the form . From the diagonalization of the dynamical matrix in Eq. (S27), we obtain the eigenfrequencies and eigenmodes of the system around the steady states. Since in our experiment |∆ c | ω 0 , a clear separation between the photon-like and atom-like polariton modes exists. The polariton mode |1 discussed in the main text is the atom-like mode obtained by linearizing around the normal phase. The eigenfrequencies associated to this polariton mode are the eigenvalues ω ± presented in Figs. 3(g,h). The ω − branch corresponds to the annihilation of a particle in the unexcited mode |0 and the creation of a particle in the polariton mode |1 . The ω + branch corresponds to the opposite process.

B. Polaritonic decay rates γ ↓(↑)
The decay rates γ ↓(↑) of the low-energy polariton modes can be derived from the diagonalization of the dynamical matrix in Eq. (S27). Here, we derive the simplified analytical expression of γ ↓(↑) given in Eq. (2) in the main text, which results from a perturbative expansion of the eigenvalues of the system in the small parameter ω 0 /κ 1, which is well justified for our experiment. We do not consider spin dissipation and use the Keldysh action formulation [14,15]. First, we bosonize the spin using Holstein-Primakoff transformation,Ŝ z =b †b − N 2 ,Ŝ + = N −b †bb , witĥ b being a bosonic annihilation operator. We then write the Keldysh action in frequency domain and integrate out the cavity degree of freedom, obtaining the spin only action where the 4-component Nambu-spinor is given byb i = (b i (ω)b * i (−ω) ), i = c, q, and the inverse Green's functions and Keldysh component are We perform a first order expansion in ω 0 /κ 1 by approximating ω + iκ ≈ iκ and solve Eq. (S31) for ω. We obtain We note that this result can be obtained also from a linear analysis after adiabatic elimination of the cavity field. The rates γ ↓(↑) describe the dissipative damping (amplification) of the polariton mode |1 , as discussed in the main text. They can be re-written in the form where ρ(ω) = κ/[(ω − ω c ) 2 + κ 2 ] is the density of states of the cavity at the frequencyω of the cavity field, cf. Fig. 1 in the main text. The expression (S34) indicates that the mechanism at the origin of damping (amplification) γ ↓(↑) is the scattering of photons from a single drive with strength η b(r) into the bath of vacuum modes provided by the cavity, accompanied by a transfer of population from mode |1 (|0 ) to mode |0 (|1 ).
The expression in Eq. (S34) is valid in the limit ω 0 κ. A more accurate estimation for γ ↓(↑) can be obtained using Fermi's golden rule [16], and the limit in which only a single drive is present, i.e., η b = 0 or η r = 0. The result is withω ± =ω ±ω 0 . The frequency of the field scattered into the cavity by each drive deviates fromω by ±ω 0 , according to energy conservation [see also Fig. 1(c) in the main text for a schematic visualization].
The correction obtained by using Fermi's golden rule becomes particularly relevant near the Dicke limit η b = η r . Without this correction, the rates γ ↓(↑) compensate each other [cf. Eq. (S34)], and the imaginary part of Eq. (S33) is zero for allη <η c , resulting in a vanishing damping rate. The result obtained with Fermi's golden rule allows to account for higher orders of ω 0 /κ, leading to a nonzero damping. By plugging Eq. (S35) into the expression for the eigenvalues Eq. (S33), we find that, in the Dicke limit, the effective damping rate of the polariton mode |1 is where we again made use of ω 0 κ. This result is in agreement with previous derivations [17] and provides insights into the physical origin of a finite effective polariton damping in the driven-dissipative Dicke model.
We point out that the theoretical results used for the analysis of excitations (Fig. 3 in the main text) have been obtained by exact diagonalization of the dynamical matrix (S27), and therefore they do not suffer from truncation effects.

Experimental protocol
To measure the evolution of the excitation spectra, we prepare a BEC of N = 9.6(4) × 10 4 atoms in m F = −1 and ramp up the coupling strengths η r,b (t) within t ramp = 9.1 ms. The experimental parameters for these measurements are ω 0 = 2π · 48(4) kHz and ∆ c = −2π · 5.8(1) MHz. While ramping up the coupling, we inject an excitation field through the cavity between t ∈ [3.0, 4.0] ms. The amplitude of the excitation field corresponds to 7.2(1) intra-cavity photons; its frequency is chosen to be close to the polariton resonanceω + ω 0 . By this method, we typically transfer < 10% of the atomic population in the excited polariton mode |1 . After the end of the excitation pulse, the polariton mode evolves freely according to the dynamics of the open system. We monitor this free evolution in real time via the spectrum of the associated photon field.

Relation between the polariton dynamics and the cavity spectrogram
Here, we show how the dynamics of the polaritonic excitation can be derived from the associated cavity field, detected with a heterodyne setup, as done in Sec. IV of the main text. Since the population in mode |1 prepared by the excitation pulse is small, we can linearize the mean-field dynamics of the system around the normal phase. By substituting (α Re , α Im , X 0 , Y 0 , Z 0 ) = (0, 0, 0, 0, −1/2) in Eq. (S27) and linearizing the pseudo-spin conservation δZ = −(X 0 δX+Y 0 δY )/Z 0 = 0, we obtain a 4x4 stability matrix M in the basis (δα Re , δα Im , δX, δY ) T . Diagonalization of this matrix leads to where ∆ + = −∆ * − and ω + = −ω * − because M has real coefficients. The pairs ∆ ± , ω ± correspond to the photon-like and atom-like polariton modes, respectively, with |Re[∆ ± ]| ≈ |∆ c | and |Re[ω ± ]| ≤ ω 0 [18]. The eigenfrequencies depend on the couplingsη, ∆η, which are time-dependent in our experimental protocol. The coupling sweeps are however slow enough to allow the system to evolve adiabatically, and the polariton modes evolve independently of each other. By decomposing into polariton modes, the evolution of the cavity field quadratures can be written as where c lj are complex coefficients, with c l+ = c * l− and l = p, a denoting photon and atom, respectively. Due to the large separation |∆ c | ω 0 , the photon-like mode cannot be excited by the external pulse on resonance with the atom-like mode, and can be neglected. We thus find The initial conditions c 0 , φ 0 are determined by the externally induced excitation process. An analogous expression holds for α Im , allowing to directly relate the real and imaginary part of the polariton frequency to the measured cavity output. Specifically, from Eq. (S39) we find that, since Im[ω ± ] < 0, the amplitude of the field n ph decays as n ph ∝ e 2Im[ω±]t , with a corresponding 1/e-decay time τ = −(2Im[ω ± ]t) −1 , which is the result used in Sec. IV of the main text.

Data processing and comparison to theory
The experimental values of the excitations lifetime shown in Fig. 3(f) of the main text are extracted from the heterodyne measurement of the cavity output in the following way. We consider the photon number spectrogram n ph (ω, t) and set t = 0 to the end of the excitation pulse. We first extract the time t c at which the superradiant phase transition occurs by integratingñ ph (ω, t) on a frequency interval −ω lim ≤ ω −ω ≤ ω lim with ω lim = 2π ·10 kHz ≈ 0.2ω 0 and setting a transition threshold of 100 intra-cavity photons. Such threshold is large enough to capture only the coherent field in the superradiant phase, but still small enough to detect reliably the critical point. For large ratios ∆η/η leading to the dissipation-stabilized normal phase, t c is taken as the average of the transition times extracted where the superradiant phase builds up. To study the time evolution of the polaritonic excitations, we integratẽ n ph (ω, t) on a larger frequency range ω lim = 2π · 150 kHz ∼ 3ω 0 to get the photon trace n ph (t), with a time resolution of 10 µs. We extract the lifetime τ from the cumulative signal s(t) = t 0 n ph (t )dt . For n ph ∝ e −t/τ , s(t) takes the form s(t) = s max (1 − e −t/τ ). Since for most of the data points the lifetime is significantly shorter than the transition time t c , we determine τ by the time at which s(t) has reached a fraction (1 − e −1 ) of its maximum below t c . The results of this estimation are in agreement with the ones obtained from a fit of s(t), but more robust especially for weak signals. As an exception, for the single dataset at ∆η/η = 0 the condition τ t c is not fulfilled, and τ is extracted from a fit of s(t) with the exponential model. We neglect the experimental realizations in which the atomic response during the excitation pulse is below the noise level. The data shown in Fig. 3(f) of the main text are averaged values of τ over 10 to 25 realizations, with the error bar representing the maximum between standard error of the mean and the time resolution of the photon trace.
We compare the values of τ extracted experimentally with the theoretical expectations from the excitation eigenfrequencies of the system. According to the description given above, we expect that n ph (t) ∝ e 2Im[ω±]t , which provides an analytical estimation of the lifetime τ an = − (2Im[ω ± ]) −1 . If the couplings vary in time, the decay of the photon number n ph (t) is in general non-exponential. However, a meaningful estimation for the measured lifetime τ at large enough imbalance ratio ∆η/η 0.05 is provided by the value τ an obtained for the instantaneous couplingsη, ∆η just after the excitation pulse. [19] To account for processes leading to dephasing of the individual atomic spins, such as collisions, we introduce a phenomenological atomic dephasing rate Γ φ , as described in SIV A. In Fig. 3(f) of the main text, the blue shaded region shows the lifetime estimated from the eigenvalues, and assuming Γ φ = 0 (upper bound) and Γ φ = 2π · 500 Hz (lower bound, corresponding to the estimated collision rate in the total trapping potential).
For a closer comparison to the experiment, a numerical simulation including time varying coupling is performed, using the method described below in SIV D [red shaded region in Fig. 3(f) of the main text]. The shaded region includes results of simulations performed for different initial phase φ probe ∈ [0, 2π) of the excitation drive, which we do not control in the experiment, and an atomic dephasing rate Γ φ varying in the same interval described in the previous paragraph.

D. Numerical simulations of the probing method
We simulate numerically the experimental protocol that we implemented to probe the excitation spectrum of our system, and described in Sec. IV of the main text. For this purpose, we extend the theoretical model and incorporate an additional intra-cavity probe beam. We consider a classical z-polarized electric field propagating along the cavity axis E probe (t, x) =Ẽ probe n probe (t) cos(krecx)e −iω probe t−iφ probe , with ω probe =ω + ω 0 + δ probe . Here, n probe (t), φ probe and δ probe are the average intra-cavity photon number, relative phase and detuning with respect to the cavity field associated with the polariton branch ω + at low couplings (ω + ω 0 ). For the numerical simulation, we choose (ω0, ω, κ, Γ) = 2π · (48 kHz, 5.8 MHz, 1.25 MHz, 100 Hz) and N = 9.6 × 10 4 atoms. The couplings are increased via an s-shaped ramp within 9.1 ms toη < 2π · 1.16 kHz at fixed ratio ∆η/η. Moreover, a blue detuned probe with δ probe = 2π · 2 kHz, n probe = 7.2 photons and φ probe = 0 illuminates the system between −1 ms < t < 0 ms.
Moreover,Ẽ probe is the electric field per photon in this beam.
Following an analogous approach to the derivation of the HamiltonianĤ in Sec. SII, we obtain a time-dependent many-body Hamiltonian describing the interaction of the light-matter system with the probe field: H exc =Ĥ MB + 4 η n probe (t) sin [(ω 0 + δ probe )t + φ probe ]Ĵ x + 4 ∆η n probe (t) cos [(ω 0 + δ probe )t + φ probe ]Ĵ y . (S41) The probe beam drives the atomic coherencesĴ x,y , similar to cavity-enhanced Bragg spectroscopy techniques [20]. Hence, we expect to coherently transfer non-negligible atomic populations to the excited state if we approach the low-coupling two-photon resonance δ probe → 0.
We derive mean-field equations of motion from the Hamiltonian in Eq. (S41) and numerically solve them to simulate the experimental conditions of our excitation probing method, cf. Fig. 3 in the main text. We plot the resulting spectrograms of the PSD in Fig. S3. The coupling imbalances ∆η/η and the parameters (ω 0 , ∆ c , κ, N ) are chosen in accordance to the experimental observations reported in the main text [Figs. 3(b)-(e)]. Depending on the choice of ∆η/η we transfer between 8% and 12% of the atoms to the excited state at the end of the probe pulse. The results from the numerical simulations are in good agreement with the experimental results at different coupling imbalances: while in the Dicke limit the excitations are long-lived and a complete mode softening towards the superradiant phase is observed, the polaritonic excitations damp faster for increasingly large coupling imbalances, as observed in the experiment.

A. Data Processing
For every single hysteresis measurement, we fix the coupling strengthη and record the photon number n ph extracted from the heterodyne detector as a function of the other coupling strength ∆η. To further process the raw data, we smoothen it by applying a moving average over 51 subsequent points. Subsequently, we define a threshold photon number for the detection of a stable superradiant phase.
In order to extract the threshold of the superradiant region, we set a threshold of 36 mean photons, which is 5 times above the noise level. Around this value, the width of the hysteresis is independent on the choice of the threshold. We determine and compare the critical coupling strength for the forward (backward) path, as shown in purple (orange), in Fig. 4 of the main text and in Fig S4 below. The critical couplings are used to map out the hysteresis region for different coupling strengthsη. For every data point we take on average 15 measurements, with at least 12 and at most 18 repetitions. ). An artificial offset inη has been introduced between the forward and backward paths for better visibility. (c) Boundaries of the normal phase detected during the forward (purple) and backward (orange) path for differentη. The position of the boundaries is determined from the photon traces by setting a threshold of 36 mean intra-cavity photons, as indicated with a gray line in (a). As guide to the reader, in the background of (b,c), the phase diagram from analytical calculations shows the region of stable normal phase (white), stable superradiant phase (dark green) and bistability (light green). The theoretical boundaries have been rescaled to the experimental data, with a single factor applied to both couplings. This scaling factor is chosen to overlay the theoretical phase boundary between the superradiant and bistability regions (dashed red line) and the corresponding experimental datapoint with the largest couplingη. For these measurements, we employ the following experimental parameters N = 1.10(8) × 10 5 , ∆c = −2π · 3.0(5) MHz and ω0 = 2π · 40(5) kHz.

B. Hysteresis loops in opposite directions
Our hysteresis measurement is potentially sensitive to atom loss and heating during the experimental protocol. These processes affect the collective atom-cavity coupling, and consequently shift the stability boundaries of the different phases. To ensure that the measured bistability region is not substantially biased by a variation of the collective coupling due to these effects, we complement the measurement shown in Fig. 4 of the main text with the result of a hysteresis loop performed in the opposite direction, see Fig. S4.
Hysteresis is observed also in this second measurement protocol, confirming that the effect of atom loss and heating is not substantial.