Superglass formation in an atomic BEC with competing long-range interactions

The complex dynamical phases of quantum systems are dictated by atomic interactions that usually evoke an emergent periodic order. Here, we study a quantum many-body system with two competing and substantially different long-range interaction potentials where the dynamical instability towards density order can give way to a superglass phase, i. e., a superfluid disordered amorphous solid, which exhibits local density modulations but no long-range periodic order. We consider a two-dimensional BEC in the Rydberg-dressing regime coupled to an optical standing wave resonator. The dynamic pattern formation in this system is governed by the competition between the two involved interaction potentials: repulsive soft-core interactions arising due to Rydberg dressing and infinite-range sign changing interactions induced by the cavity photons. The superglass phase is found when the two interaction potentials introduce incommensurate length scales. The dynamic formation of this peculiar phase without any externally added disorder is driven by quantum fluctuations and can be attributed to frustration induced by the two competing interaction energies and length scales.


I. INTRODUCTION
The controllability of individual atomic systems increased tremendously over the past decades [1][2][3][4][5][6]. This enables efficient trapping and cooling of atomic gases or individual atoms, controlling individual photons and tailoring interactions between atoms and photons. One particularly challenging and active modern research direction in this realm are systems with tailored long-range interactions [7]. Prominent examples are, among others, dipolar Bose-Einstein condensates (BECs) [8][9][10][11], ultracold atomic gases in cavities [12][13][14] and individually trapped atoms with Rydberg interactions [15][16][17]. These systems allow for the study of artificial quantum matter in a well-controlled and tunable environment.
Over the last decade the focus of research was on investigating long-range interacting systems with crystalline properties, i. e., periodic systems with long-range order. In this realm first experimental observations of quantum phase transitions like the superfluid to Mott insulator transition [18], the formation of supersolid phases of matter [19][20][21][22][23][24][25][26] or non-trivial spin phases [15,27,28] were realized and understood on a fundamental level. In recent years, the exploration of systems generating more complex patterns and phases became a leading research direction in many-body quantum physics. Particular focus hereby lies on quantum glasses [29][30][31][32], many-body localization [33,34], spin-liquids [17,35] and quasicrystals [36][37][38]. Realizing such complex systems establishes a path towards a deeper understanding of interacting quantum many-body systems, which ultimately fosters potential future applications in modern quantum technologies. Many of the effects mentioned above are related to frustration, i. e., the inability to satisfy all constraints imposed by the system's constituents at the same time. * stefanostermann@g.harvard.edu Figure 1. (a) Considered setup. A Rydberg-dressed BEC is trapped and confined to two dimensions in an optical resonator with lasers impinging from the side. The BEC atoms form an effective V-level structure. The transition |g ↔ |e is coupled to the cavity at a strength Gp. The transition |g ↔ |r is driven by an additional laser in the Rydberg dressing regime, implying that the high-lying Rydberg state |r is only very weakly populated. The Rydberg dressing imposes a long-range soft-core interaction potential between the atoms, see panel (b). In addition, the cavity photons induce nontrivial infinite range interactions between atoms, see panel (c). The top curves in panels (b) and (c) show cuts along the x-axis at y = 0.0λc for the two interaction potentials.
In this work we introduce and study a particularly clean and experimentally well controlled system which can give rise to frustration: a Rydberg-dressed BEC of neutral atoms trapped in an optical standing wave resonator, see Fig. 1(a). We elucidate some of the archetypal effects of the contest among the cavity-induced infinite range interactions and long-range Van-der-Waals interactions among Rydberg atoms and ultimately leading to glassy behavior. Previous work that combined Rydberg atoms and cavities focused on either generating large optical nonlinearities for single photons [39][40][41] or spin models on a lattice where the motional degrees of freedom are frozen [29,41,42]. Here, we focus on the spatial dynamics and resulting density self-ordering of the gas under competing long-and infinite-range interactions. Both constituents considered in this work exhibit a self-ordering phase transition due to roton mode softening [43,44]. These soft roton modes are a direct result of the respective long-range interactions. However, the two types of interactions have substantially different properties. The Rydberg dressing results in spherically symmetric longrange soft-core interactions [see Fig. 1(b)] and the cavity photons induce anisotropic infinite-range periodic interactions between atoms [see Fig. 1(c)]. This distinguishes our work from related studies focusing on a spherically symmetric interaction potential imposing two different length scales [45][46][47]. In fact, while these works use a synthetic theoretically motivated long-range interaction potential, our study is based on an experimentally realistic configuration leading to the non-trivial physics presented below.
As we will show in the following, the particular properties of the light-matter coupled system, i. e., the competition between the two non-trivial interactions at different length scales and the nonlinear interaction between the cavity mode and the BEC lead to the dynamic formation of intriguing phases of quantum matter. In particular, we show how this system can be used to realize a stable glass phase, i. e., an amorphous solid with no long-range density order. The formation of this intriguing phase of quantum matter is attributed to geometrical frustration induced by the two competing interaction length scales and energies. Due to the additional superfluid properties of the BEC the resultant glassy phase is often also referred to as a superglass [30,[48][49][50][51]. The formation of this disordered glassy state in our system is not triggered by an externally imposed disorder or an external lattice geometry which results in frustration. It is the interplay and particular nature of the two considered model constituents (cavity photons in combination with Van-der-Waals interactions) which gives rise to the phases presented below.

II. MODEL
We consider a pancake-shaped BEC of N atoms confined in two dimensions and trapped inside an optical standing-wave resonator, see Fig. 1(a). The BEC atoms' ground state |g is coupled via the cavity to a low-lying excited state |e and, simultaneously, to a highly excited Rydberg s-state |r , realizing the V-configuration illustrated in Fig. 1(a). The BEC atoms are pumped with a laser of frequency ω p impinging from the side, transverse to the cavity axis, which drives the transition |g ↔ |e at a Rabi frequency Ω p . This transition couples to a single cavity mode with a coupling strength G p . We consider the dispersive regime implying that the detuning ∆ p = ω p − ω ge is large compared to the Rabi frequency Ω p (ω ge is the transition frequency between |g and |e ). The transition |g ↔ |r is driven by an additional laser (with frequency ω R ) at a Rabi frequency Ω R in the Rydberg dressing regime implying a large detuning ∆ R Ω R with ∆ R := ω R − ω gr (ω gr is the transition frequency between |g and |r ). This allows the adiabatic elimination of the Rydberg state |r resulting in a effective long-range two-body interaction potential [52][53][54][55][56] U Ryd (r, r ) =C which arises due to the strong Van der Waals (VdW) interactions ∝ C 6 /|r| 6 (r ∈ R 2 ) between Rydberg atoms. The parameters in Eq. (1) are defined asC 6 := . An exemplary plot of the functional dependence of this long-range interaction potential is shown in Fig. 1(b). It was shown in previous work that a mean-field treatment suffices to capture the main features of cavity self-ordering [12,13,57] and Rydberg crystallization induced by the long-range interaction potential given in Eq. (1) [43,54,58]. Therefore, we focus on a mean-field treatment. Here, we present the most important equations used throughout this work. A more detailed discussion of the model can be found in Appendix A. The dynamics is governed by two coupled equations. The equation for the BEC order parameter ψ(r, t) is given as and the dynamics of the mean-field cavity mode amplitude α(t) is governed by In Eq. (2b) we introduced the cavity bunching parameter B[ψ] := V dr|ψ(r, t)| 2 cos 2 (k c x) and the cavity mode order parameter θ[ψ] := V dr|ψ(r, t)| 2 cos(k c x) cos(k c y). The latter is the crucial quantity for understanding the self-ordering phase transition since the cavity mode can only take non-zero values if θ[ψ] = 0, whereas the former only accounts for a cavity resonance frequency shift due to the BEC density. The potential depth of the cavity potential generated by inter-cavity photon scattering is defined as U 0 := G 2 p /∆ p , and the effective pump strength is η := Ω p G p /∆ p in Eq. (2a). In addition, k c = 2π/λ c denotes the cavity wavenumber where λ c is the cavity resonance wavelength. We also introduced the detuning of the pump laser frequency ω p with respect to the cavity resonance frequency ω c as ∆ c := ω p − ω c and the cavity decay rate κ.
A particularly simple and insightful model can be obtained in the low-energy regime by adiabatically eliminating the cavity mode (for details, cf. Appendix A), which results in a new effective equation for the BEC dynamics In Eq. (3) we introduced the cavity two-body interaction potential induced by the cavity photons as is the effective interaction strength. An exemplary plot of this infinite-range interaction potential is shown in Fig. 1(d). It is the competition between the two significantly different types of long-range interaction potentials U Ryd and U cav shown in Figs. 1(b) and (c), which contributes to the non-trivial emergent phases presented below. Note, that the simplified model given in Eq. (3) provides a good intuitive picture about the fundamental interactions induced by the two model constituents. However, this ap- Figure 3. Sketch of the expected phase diagram based on the excitation spectrum. We distinguish four regions: I -no roton softened, i. e., the homogeneous solution is stable, IIthe cavity roton touches zero (indicated by red circle), i. e., acquires an imaginary part, III -the Rydberg roton is unstable at k = k Ryd and IV -both rotons are unstable. The distance ∆k between the two k-values at which the rotons soften can be tuned by tuning Rc in Eq. (1).
proximate model does not cover all features arising from the nonlinear coupled dynamics of the cavity mode and the BEC [see Eq. (2)], as discussed below.

III. ROTON INSTABILITIES
Roton-induced instabilities were first introduced and studied for superfluid Helium-4 [59] but they are a common feature of systems with long-range interactions. A dynamical instability of a certain nonzero k-mode in the collective excitation spectrum usually results in a phase transitions from homogeneous to periodic order. The Bogoliubov excitation spectrum for plane-wave excitations on top of a homogeneous condensate for a system with two-body long-range interaction potentials reads [60] where U Ryd (k) and U cav (k) denote the Fourier Transforms (FTs) of the respective interaction potentials U Ryd (r, r ) and U cav (r, r ). We also introduced the particle density ρ := N/V with V being the total volume of the BEC. The FT of the Rydberg interaction potential cannot be expressed in a comprehensive analytical formula. The FT of the cavity interaction potential (4), however, is given as U cav (k) = I i,j∈{0,1} δ kx,(−1) i kc δ ky,(−1) j kc , where δ k1,k2 denotes the Kronecker-delta. In Fig. 2 an exemplary excitation spectrum is shown. The spherically symmetric Rydberg interaction potential induces roton softening along a circle with a radius k Ryd [indicated by the dashed orange line in Fig. 2(a) and (b)] . The cavity interaction potential, however, results in four δ-shaped peaks/minima which are indicated by the white circles in Fig. 2(a) and (b), and are more prominently visible in panels Fig. 2(c) and (d). Hence, while the two interaction potentials considered in this work both result in dynamical instabilities induced by roton mode softening, the properties of the related excitation spectra differ significantly. The the cavity roton positions are always fixed by the cavity wavelength and are located at a value k cav := √ 2k c while the value of k Ryd can be tuned via the parameter R c in Eq. (1). Since each unstable k-mode can be associated with a characteristic length scale via k Ryd = 2π/l Ryd and k Ryd = 2π/l cav , the role of different emergent length scales on the final steady or ground state can be analyzed by changing R c correspondingly. This can either be achieved by dressing to a different Rydberg state (C 6 scales like n 11 , where n is the principal quantum number) or by changing the detuning ∆ R .
The homogeneous solution gets unstable toward a periodically ordered pattern if the excitation spectrum acquires an imaginary part at a non-zero k-value, i. e., min Re(ε(k))| |k|>0 = 0. This results in two critical val-uesC crit 6 (critical Rydberg interaction strength) and η crit (critical effective cavity pump strength). Since the FT of the cavity interaction potential is solely given as the sum of four Kronecker-delta functions it is possible to find an analytical expression for the threshold for the pure cavity self-ordering case. SettingC 6 to zero and applying the threshold condition for k = k cav in Eq. (5) results in which is the known threshold for self-organization of a BEC in a cavity in 2D [13]. These insights yield a first intuitive phase diagram shown in Fig. 3. We expect the phase diagram to exhibit four different regions. In each region either no, one or two of the respective rotons are unstable (indicated by red circles in Fig. 3). In addition to the two interaction strengths tuned in the phase diagram, the system features another relevant free parameter -the Rydberg interaction range R c . In the following we will restrict ourselves to three cases: A sketch of the unstable roton length scales in k-space for all three cases is shown in Fig. 4. We restrict ourselves to cases where k Ryd < k cav because only in this case the additional Rydberg length scale competes with the length scale set by the cavity potential (see also Appendix B). Note that the first two cases (i) and (ii) are special cases because the two wavenumbers are either the same [case (i)] or k Ryd is commensurate with k c , i. e., the fundamental length scale along the cavity axis [case (ii)]. Case (iii), however, corresponds to one example of the most general incommensurate scenario. We chose a factor two for the Figure 4. The three cases considered in this work. The red dots mark the four k-modes which are unstable due to the cavity interaction potential (see Fig. 2). The blue circle indicates all k-values which are unstable due to the Rydberg dressing potential. case (i): k Ryd = kcav, case (ii): k Ryd = kcav/ √ 2 = kc and case (iii): k Ryd = kcav/2. The cases (i) and (ii) represent special cases, whereas case (iii) corresponds to one example for the most general scenario.
wavenumber ratio in case (iii) but the same qualitative results as presented below were obtained for a variety of different wavenumber ratios.

IV. SELF-CONSISTENT GROUND STATE
One major element of the studied system is the interplay between the dynamic cavity field and the BEC dynamics. While the simplified model given in Eq. (3) provides a good intuitive picture about the features and roles of the two different interactions, it does not fully capture the competing time scales in the system. Hence, we now analyze the full system given in Eq. (2). We calculate the self-consistent ground state for the set of Eqs. (2) via a variational approach. To this end, we employ a selfconsistent iterative algorithm: 1) choose a random initial steady-state value α ss . 2) plug this value into the GPE equation (2a) and perform an imaginary time evolution (t → −iτ ) to find the lowest energy state. 3) calculate a new steady-state of Eq. (2b) by setting ∂ t α(t) = 0 and solving for α. We iterate until we find a state where |α ss | lies within a convergence radius of 10 −6 . Note that this self-consistent algorithm is not necessarily convex. In general, the energy functional which results in the GPE in Eq. (2a), can have multiple local minima for a given value of α. Whether one finds the global minimum or not depends on the initial condition for the imaginary time evolution. Therefore, we perform this algorithm for various different initial conditions and compare the obtained values of |α| and the corresponding ground state energies by evaluating the energy functional in Eq. (7). Only if all initial conditions result in the same energy and mode amplitude the energy minimization problem is convex. For details we refer to Appendix B.

A. Case (i) -equal length scales
We first consider the special case where k Ryd = k cav (see Fig. 4), which implies equal length scales l Ryd = l cav (R c = 0.51λ c ). In this case, the regions II and IV in the phenomenological phase diagram in Fig. 3 are expected to merge. This is verified by the numerical phase diagram shown in Fig. 5(a). It exhibits three different phases a homogeneous phase (HOM), a triangular phase (TRI) [see Fig. 5(b)] and a checkerboard lattice phase (CB) [see Fig. 5(c)]. The TRI state only forms if the Rydberg induced roton is unstable, i. e., beyond the critical Rydberg interaction strength indicated by the red dashed line in Fig. 5(a). For increasing cavity pump, the system organizes in a checkerboard pattern, which is the configuration realized in cavity self-organization without Rydberg dressing [12,13]. Note, however, that the softening of the Rydberg roton affects the threshold for cav-ity self-ordering (blue line in Fig. 5(a)). This implies that the effect of the Rydberg roton can be observed already far below the critical Rydberg interaction strengthC crit 6 . This is a remarkable result as the observation of density modulation due to Rydberg dressing has proven elusive in experiments [56] due to experimental complications such as unfavorable scaling with the density and long time scales required for experiments. Our results suggest a way to observe effects of Rydberg dressing indirectly for much smaller Rydberg fractions (i. e., smallerC 6 ). In Fig. 5(d) we show the obtained ground state energy obtained from Eq. (7) and the cavity mode amplitude for all seven employed initial conditions. We shifted all energies with respect to the homogeneous state energy E hom which is obtained by evaluating Eq. (7) for α = 0.0 and a homogeneous BEC density. Energies are given in units of the recoil energy E rec := 2 k 2 c /(2m) = ω rec throughout this work. We find that all curves for all initial conditions lie on top of each other, implying that this particular iterative energy minimization for case (i) is convex in ψ and α.

B. Cases (ii) and (iii) -different length scales
The main property of the considered system is that two interactions with different functional dependence and different intrinsic length scale are combined. Here, we show that this feature results in highly non-trivial physics which goes beyond what is achievable with only a single model constituent. Fig. 6 exhibits the two obtained ground state phase diagrams together with the additional non-trivial density distributions realized in these regimes. Since now two fundamental length scales become unstable the phase diagram resembles the intuitive picture from Fig. 3 and contains four different phase regions. The critical values obtained via Eq. (5) in section III, however, do not coincide with the numerical values. This shows that the simplified model discussed in section III for the adiabatically eliminated cavity dynamics [see Eq. (3)] provides a good first intuition about the contributing interactions and the expected phase diagram, but fails for different competing length scales due to a break down of the adiabatic solution. This is due to the influence of the additional Rydberg length scale on the mode dynamics which is taken into account in the full numerics resulting in the phase diagrams shown in Fig. 6. The instability in this system is actually induced by the combined collective instability of the BEC density and the cavity mode amplitude. The fluctuations of the BEC couple to the fluctuations of the cavity mode and vice versa. This results in actual thresholds for the transitions TRI ↔ SQ and CB ↔ SQ in Fig. 6(a) and TRI ↔ RC and CB ↔ RC in Fig. 6(c) that are below the values anticipated by the simplified excitation spectrum in Eq. (5). This again could facilitate the experimental observation of effects induced due to Rydberg induced roton mode softening at smaller Rydberg interaction strengthsC 6 . For case (ii) we find that, in addition to the TRI and CB phases, a square (SQ) lattice phase emerges [cf. Fig. 6(e)]. This phase is a direct result of the non-trivial interplay between the two different interactions at these particular length scales. All seven initial conditions again result in the same ground state energies and mode amplitudes [see Fig. 6(b)]. This feature changes in the incommensurate case (iii) -k Ryd = k cav /2. We identify an additional rotated chain (RC) phase [see Fig. 6(c) and (f)], where the iterative energy minimization in the selfconsistent algorithm is no longer convex. There are two possible realizations of the RC phase. In Fig. 6(f), we show one possible realization. However, the given density pattern rotated by π/2 is another ground state with the same energy. This can be seen by the two degenerate lowest energy curves in Fig. 6(d) which correspond to the two possible realizations of the RC phase. The two degenerate states maximize the population of the cavity mode [see Fig. 6(d)] while fulfilling the length scale restrictions imposed by the Rydberg interactions, cf. Appendix B for details. Such degeneracy in combination with non-convex features is a prime indicator for geometrical frustration ultimately leading to glassiness as we show in the next section. Note that we restricted our analysis on three cases for a concise presentation of the results. It should be mentioned that cases (i) and (ii) are indeed two special cases, which lead to a convex problem. Any other choice (case (iii) is one possibility) of length scale ratios results in a non-convex energy landscape. The non-convex na-ture of the energy landscape in case (iii) is the crucial feature which results in the superglass formation presented below.

V. DYNAMICS
It is a priori unclear whether the self-consistent ground state can be reached by dynamically evolving Eq. (2) in real time. While this is usually the case for systems with a unique lowest-energy state, ground-state degeneracy in combination with higher lying energy states may block the dynamical realization of the RC phase in case (iii). In fact, we will show below that the steady-state for case (iii) is an amorphous solid with no long-range density order.
In Fig. 7, the time evolution of the cavity mode amplitude and the total energy together with the steady-state density distribution at t = 500/ω rec is shown for case (ii). We find that the square lattice obtained in the previous section is a stable steady state, which can be obtained dynamically when starting from a homogeneous BEC. In contrast, the steady state for case (iii) is stable but not the rotated chain state as the self-consistent phase diagram suggests. The steady state is an inhomogeneous density pattern with no long-range density order. In Fig. 8(a), we show three outcomes for three independent runs of the time evolution for the same parameters.  Every run results in a different density distribution. This feature, in addition to the BEC's superfluidity, qualifies this phase as a so-called superglass [48][49][50]. The glassy behavior can be directly understood from the nonconvex properties and the degeneracy of the energies found in the previous section [see Fig. 6(f)]. Each run results in a different mixture of the degenerate ground states giving rise to the glassy steady state. Remarkably, however, the resulting different density distributions still have exactly the same steady-state cavity mode amplitude. Hence, the different glassy states are not distinguishable via the cavity fields. Note that the formation of this superglass state is a direct consequence of the two particularly chosen model constituents and cannot be straightforwardly realized by combining any two long-range interaction potentials.

VI. CONCLUSIONS AND OUTLOOK
We showed that combining the spherically symmetric long-range interactions induced in a Rydberg-dressed BEC with infinite-range cavity-induced interactions establishes a versatile platform for studying the interplay between different interaction types and landscapes on the phases of quantum matter. Most strikingly, we point out a way towards realizing a superglass phase in a robust and well controlled system. The disorder in this superglass phase is not imprinted from outside via some external potential or lattice but it is solely driven by quantum fluctuations and the interplay of the two interaction potentials. This dynamical, fluctuation induced formation of a superglass phase is closely related to the theoretically predicted phase transition from a liquid to a superglass after a temperature quench in Helium-4 [48], which has not been realized experimentally or found elsewhere yet. Our findings hint towards a realizable system to dynamically study phase transitions from a homogeneous or solid (CB phase) to a superglass by tuning the Rydberg interaction strength. The studied system is highly non-linear and the formation of the intriguing phases presented in this work involves the interplay between the BEC wave function ψ and the cavity mode α. This gives way to ease the observation of Rydberg induced instabilities experimentally because the non-linear coupling of atomic density fluctuations to the cavity mode facilitates the softening of the Rydberg roton. Our findings also have potential applications in modern quantum technologies: The formation of a density ordered pattern in this hybrid system ultimately resembles a quantum optimization problem. Hence, the system could be tailored to particular optimization problems. In addition, glassy states are promising candidates for efficient quantum memories.
Our work lays the ground for a variety of further studies combining long-range interacting systems theoretically as well as experimentally. This research avenue is fostered by the new experimental and theoretical opportunities opening up in recently established experimental setups all over the world. Our proposed setup should be realizable in state-of-the-art cavity QED setups. In particular, the realization of the glassy phase where the parameters are favorable for experiments could readily be realized. However, the setup presented here is not the only system which is expected to exhibit such intriguing properties. Other promising avenues are combining other systems exhibiting light-induced instabilities in free space [32,[61][62][63][64][65][66] with Rydberg dressing or long-range dipole interactions. More control over a wide range of realizable density patterns is expected if one assumes even more complex cavity-induced interaction potentials or Rydberg interactions. This can be achieved by changing to more complex cavity geometries such as multi-mode resonators [67][68][69] or dressing to spherically asymmetric Rydberg p-states [70,71]. From a quantum optical viewpoint the dynamic transverse coupling of Rydberg atoms to a cavity could result in a highly nonlinear mechanism to enhance the achievable nonlinearities with state-of-the-art setups. The platform presented in this work could also be extended to include spin degrees of freedom, and, serve as a viable tool to gain deeper understanding of spin glasses or even spin liquids [17,72,73].
While this work focused on highlighting the remarkable features of this system and their intriguing consequences, several directions remain open for further research. One open and exciting question is the connection of the formation of the glassy phase to Anderson localization or even many-body localization. Also, the role of beyond mean-field effects and temperature needs further investigation. To answer these questions, alternative theoretical techniques, which go beyond the scope of the present manuscript, have to be developed and applied [74][75][76][77]. In any case, the results presented in this work open up exciting avenues in the growing research field of hybrid quantum systems with long-range interactions. In the following we provide a concise derivation of the model presented in section II. For a more detailed discussion of the underlying physics we refer to Refs. [79] and [54]. The full many-body Hamiltonian for the considered system can be written as the sum of five Hamil- Here, Ψ (Ψ † ) are annihilation (creation) operators of bosonic ground state atoms, i. e., Ψ(r), Ψ † (r ) = δ(r − r ) and a (a † ) are annihilation (creation) operators of the photonic cavity mode fulfilling a, a † = 1. U p := Ω 2 p /∆ p is the potential depth of the lattice generated by the two interfering pump beams and U 0 := G 2 p /∆ p is the depth of the optical potential generated by photon scattering from the atomic density distribution inside the cavity. The effective cavity pump strength is given as η := Ω p G p /∆ p and k c = 2π/λ c denotes the wavenumber of the cavity mode. The Hamiltonian in Eq. (A1d) takes into account local two-body interactions between atoms and is omitted throughout this work. This is a reasonable assumption since this interaction strength can be tuned to be small via, e. g., a Feshbach resonance such that cavity and Rydberg induced long-range interactions dominate the dynamics [54,80]. The dynamics of the hybrid atom-cavity system is governed by the Heisenberg-Langevin equations i ∂ t Ψ(r, t) = [Ψ(r, t), H] and i ∂ t a(t) = [a(t), H] − i κa(t) where the decay of the cavity mode at a rate κ is included. Calculating these equations of motion and performing a mean-field approximation via Ψ(r, t) → Ψ(r, t) := ψ(r, t) and a(t) → a(t) := α(t), results in the two coupled cnumber equations for the BEC order parameter ψ(r, t) and the cavity mode amplitude α(t) given in Eq. (2) of the main text.
To adiabatically eliminate the cavity mode we calculate the equation of motion for the field operator a and solve for its steady-state (∂ t a = 0). This results in a ss = η drΨ † (r, t) cos(k c x) cos(k c y)Ψ(r, t) ∆ c − U 0 drΨ † (r, t) cos 2 (k c x)Ψ(r, t) + iκ .
(A2) Apart from the fundamental cavity parameters η, ∆ c and κ the steady-state value for the cavity mode is determined (A4) These two quantities obviously depend on the BEC state which exhibits the non-trivial coupling between the cavity mode and the BEC. While B only acts as an effective shift of the cavity resonance frequency that comes into play as soon as a ss = 0, Θ is the crucial parameter when it comes to understanding the cavity self-ordering phase transition. The cavity mode is nonzero if Θ = 0. Hence this parameter is crucial for the instability described in the main text. To simplify the model we replace B with its value for the homogeneous condensate B = N/2 while keeping the full functional dependence of Θ. To eliminate the cavity field we plug the resultant steady-state solution into the many-body Hamiltonian. Since we are solely interested in a model capturing the dynamic instability which is governed by terms ∝ cos(k c x) cos(k c y) (see argument above) we keep only these terms. This results in the effective interaction Hamiltonian H cav int = I drdr Ψ † (r)Ψ † (r ) cos(k c x) cos(k c x ) cos(k c y) cos(k c y )Ψ(r )Ψ(r).
where we used the symmetry of the cosine function in the second line and introduced the cavity induced interaction strength as The many-body Hamiltonian for the BEC then reduces toH = − drΨ † (r) 2 ∇ 2 2m Ψ(r)+H cav int +H Ryd int . The meanfield equation of motion given in Eq. (3) of the main text is found by calculating the Heisenberg equation of motion for Ψ(r, t) and again performing the mean-field approximation.
Note that in both GPE equations presented in the main text [Eqs. (2a) and (3)] we neglected the potential term ∝ cos 2 (k c y), which arises due to interference of the pump beams [see Eq. (A1a)]. This term is irrelevant for understanding the collective instabilities and the main fundamental results presented in this work also hold if this term is included as we checked by including it in the numerics. However, neglecting this term facilitates the interpretation of the results and provides a method to understand the fundamental physics from an intuitive semi-analytical standpoint as it is outlined in section III of the main text. In fact, the interference between the two pump beams could even be suppressed in real experimental setups. This can either be achieved by choosing a large enough frequency difference between the two pump beams, or by using two counterpropagating beams with orthogonal polarization. In this case, the model in Eq. (2) becomes exact.

Appendix B: Variational algorithm
To find the self-consistent ground state in section IV we employ a variational algorithm. Here we discuss this algorithm in more detail. In particular, we focus on our choice of initial conditions and argue why the chosen initial conditions are good guesses for the final ground states. In Fig. 9, exemplary plots for the total interaction potential U tot (r, r ) = U cav (r, r ) + U Ryd (r, r ), defined as the sum of the two individual interaction potentials given in Eq. (1) and (4) are shown. These potentials provide a good intuition about the anticipated ground states. Hence, we use these potentials as a guideline to define the initial conditions for the imaginary time evolution. We see that for the two special cases (i) (k cav = k Ryd ) and (ii) (k cav = k Ryd / √ 2), the interaction potentials suggest a checkerboard pattern and a square lattice as a ground state which is indeed the self-consistent ground state in these particular cases [cf. Figs. 5(d) and Fig. 6(e)]. In addition, to those cases we choose the hexagonal closed packed (hcp) lattice as an additional initial condition for all cases because this would be the ground state for a Rydberg-dressed BEC without cavity interactions [43]. From Fig. 9 we see that case (iii) is more complex than the previous two cases. In this case, the ground state could be a checkerboard lattice with larger period than in case (i) [see red dots in Fig. 9(c)] or one of the two lattice configurations indicated by the light green and blue circles in Fig. 9(c). Based on the above arguments it also becomes clear that only the cases where k Ryd < k cav results in modified results. If this condition is not fulfilled the length scale imposed by the Rydberg interaction potential is smaller than the cavity interaction potential which only modifies the corresponding interaction potential locally [see Fig. 9(a)].
Based on these intuitive arguments we generate the respective initial states ψ(r) = i φ σ (r − R i ) by superimposing Gaussian wave functions φ σ (r) = C exp[−|r| 2 /(2σ 2 )]. The resultant non-trivial initial conditions (we refrain from showing the homogeneous density) are shown in Fig. 10. We run the variational algorithm discussed in section IV for all these initial conditions for all cases (i)-(iii). We then compare the groundstate energies calculated via Eq. (7) and the resultant cavity mode amplitude to estimate whether the respective optimization problem is convex or not.