Tunable-range, photon-mediated atomic interactions in multimode cavity QED

Optical cavity QED provides a platform with which to explore quantum many-body physics in driven-dissipative systems. Single-mode cavities provide strong, infinite-range photon-mediated interactions among intracavity atoms. However, these global all-to-all couplings are limiting from the perspective of exploring quantum many-body physics beyond the mean-field approximation. The present work demonstrates that local couplings can be created using multimode cavity QED. This is established through measurements of the threshold of a superradiant, self-organization phase transition versus atomic position. Specifically, we experimentally show that the interference of near-degenerate cavity modes leads to both a strong and tunable-range interaction between Bose-Einstein condensates (BECs) trapped within the cavity. We exploit the symmetry of a confocal cavity to measure the interaction between real BECs and their virtual images without unwanted contributions arising from the merger of real BECs. Atom-atom coupling may be tuned from short range to long range. This capability paves the way toward future explorations of exotic, strongly correlated systems such as quantum liquid crystals and driven-dissipative spin glasses.


I. INTRODUCTION
Cavity QED provides strong light-matter coupling [1]. For example, exotic nonlinear optical properties arise in cavity systems with atom-mediated photon-photon interactions [2]. Realizations of topologically nontrivial states of interacting photons are within reach [3]. Adiabatically eliminating the photonic field, rather than the atomic, yields photon-mediated atom-atom interactions. These interactions may be sufficiently strong to create novel quantum phases of matter [4]. Indeed, single and fewmode cavity QED in the optical domain have already provided demonstrations of supersolidity [5,6] and exotic Mott physics [7,8] in addition to supermode-densitywave-polariton condensation [9]. Moreover, the drivendissipative, open-quantum-system nature of cavity QED can change the character of quantum phase transitions, providing a new window into quantum nonequilibrium physics [10,11].
An outstanding challenge has remained to create many-body cavity QED systems whose description requires physics beyond mean-field approximation. Doing so enables, e.g., exploration of spin glass physics beyond the Sherrington-Kirkpatrick model where meanfield, replica-symmetry-breaking solutions may no longer hold [12][13][14], or quantum liquid crystals and intertwined orders such as those found in strongly correlated materials like high-T c superconductors [15][16][17][18]. More generally, strongly fluctuating, inhomogeneous (and frustrated) systems may organize in unexpected ways and the resulting surprises may lead to a deeper understanding of how quantum matter organizes. A crucial limitation to exploring such physics using cavity QED stems from the fact that the single or few-mode cavities employed so far admit photon-mediated interactions that are all- (a) A 87 Rb BEC (red circle) is trapped at the cavity waist z = 0 at a location x1 relative to the cavity center. The transverse pump beam (red beam) propagates alongx; the system undergoes a superradiant, self-organization phase transition above a critical field strength Ωc. Photons scattered off the BEC into the modes of the confocal cavity (green) create a virtual image (not shown) of the BEC at −x1. The distance δL indicates the tunable offset of the mirror from the confocal configuration. Emission of intracavity photons can either be sent to a single-photon counter, or be imaged onto an EMCCD camera to resolve the spatial structure of superradiant emission. An absorption imaging laser for imaging BEC density travels alongŷ (not shown). (b) Two 87 Rb BECs trapped at locations x1 and x2 on opposite sides of cavity center. Images of the two BECs are created at −x1 and −x2 (not shown).
to-all in coupling [4]. The global (infinite-range) nature of these interactions necessarily implies that meanfield approximations are adequate to explain observed physics [19]. However, it has been suggested that this challenge may be met either by employing networks of single-mode cavities [20,21], using squeezed light to engineer interactions [22], or by using a single multimode cavity [13,17,23]. This work presents a realization of a multimode cavity QED system and demonstrates that such a system does indeed provide strong, tunable, and local interactions among intracavity atoms. While no beyond-meanfield physics is yet explored, we show that the crucial ingredient of local interactions is present in the system, opening the road to future investigations where beyondmean-field physics may be manifest.
We measure the interaction range versus tunable parameters by manipulating the position of Bose-Einstein condensates (BECs) within the cavity. The symmetry of our confocal cavity is exploited to measure the interaction between real BECs and their virtual images without unwanted contributions arising from the merger of real BECs. The experimental results are compared to theory, with good agreement. Furthermore, we show that the reduction in interaction range is accompanied by an increase in the effective atom-light coupling strength (g eff ) and an emergence of a continuous translational symmetry in the plane transverse to the cavity axis.
The paper is organized as follows. Section II describes in general terms how tunable-range, photon-mediated interactions arise in a transversely pumped multimode cavity QED system undergoing a superradiant, selforganization transition. Section III then describes the cavity apparatus and BEC production and manipulation. Section IV presents the experimental results while Sec. V compares these to theory. Section VI discusses in greater detail the theoretical calculation of the photon-mediated atom-atom interaction.

II. PHOTON-MEDIATED INTERACTIONS IN A MULTIMODE CAVITY
Atomic gases placed in transversely-pumped optical cavities have been shown to undergo a superradiant, selforganization transition arising from the competition between their free-particle dispersion and cavity-mediated interactions [24][25][26][27][28]. Figure 1 shows examples of transversely pumped cavities. For a pump laser red-detuned from the cavity resonance, atoms separated by a pump wavelength λ along the cavity axisẑ constructively scatter pump photons into the cavity mode, leading to a buildup of intracavity light. Conversely, scattering from atoms separated by λ/2 is suppressed. The resulting atomic light shift from the intracavity field creates an optical lattice potential that further localizes the atoms at integer-λ separations. The cavity may be interpreted to mediate a periodic, infinite-range interaction alongẑ that lowers the energy of a λ-period atomic density wave. Above a critical pump strength Ω = Ω c , the cavitymediated interaction energy of the density wave overcomes kinetic energy 2E r = h 2 /mλ 2 and the atoms selforganize into a λ-periodic pattern [29]. In doing so, the atoms spontaneously choose to localize at either the even or odd antinodes of the standing wave. Interference between the pump and cavity beams means this even/odd choice is staggered along the pump direction and leads to a 2D checkerboard lattice in the xz-plane [24,30]. Concomitantly, the atoms superradiate into the cavity. This second-order nonequilibrium phase transition is heralded both by a change in the atomic distribution [26] and by an increase in cavity emission proportional to N [25]. The momentum distribution of the atoms may be detected in time-of-flight imaging, where Bragg peaks appear at wavevectors associated with the λ-period checkerboard lattice [26,28].
In conventional Fabry-Perót cavities, i.e., those supporting a single TEM 0,0 mode near the pump frequency, the (x, y) position dependence of the interaction energy between atoms follows the Gaussian profile Ξ 0,0 (x, y) of this mode. The interaction energy vanishes at distances larger than the mode waist w 0 . We now describe the explicit form of this interaction. Atoms at position x coherently scatter pump photons into the cavity mode at a rate η = g 0 ΩΞ 0,0 (x)/∆ a , according to second-order perturbation theory, where g 0 is the single-atom atomcavity coupling rate at x = 0. This expression is valid when the atomic excited state can be adiabatically eliminated from the dynamics for sufficiently large detuning ∆ a of the pump. A virtual photon may be exchanged between atoms within a time given by the inverse of the pump-cavity detuning ∆ c . This virtual photon mediates an interaction given by (to second-order in perturbation theory) [17] U As mentioned above, this interaction energy smoothly vanishes versus distance for x−x larger than w 0 . Atoms in gases much smaller than w 0 interact with a global, allto-all coupling, up to the sinusoidal variation inẑ due to the standing-wave field modulation between the cavity mirrors. This coupling need not be global in a multimode cavity, such as a confocal Fabry-Perót resonator in which the cavity length L equals the mirrors' radius of curvature R [31]. A multimode cavity can support several Hermite-Gaussian transverse modes at the same frequency, but with orthogonal mode functions Ξ l,m (x). Ξ l,m (x) is the Hermite-Gauss function, describing the functional form of the TEM l,m mode at (x, y) position x. An atom scattering a pump photon into the cavity does so into a superposition of Ξ l,m . The weights of the superposition are given by the mode strengths at the atomic position. They are also given by any differences in detuning ∆ l,m between the (l, m) modes and the pump due to residual When δL = L − R is increased to move the system away from confocality, the Hermite-Gaussian modes of our near-confocal cavity exhibit a linear frequency dispersion with mode number: ∆ l,m = ∆ 0,0 + (l + m) . The factor S l,m , discussed in detail in Sec. VI, accounts for the overlap between the atomic density wave and the photon mode profiles along the cavity axis. Due to the nature of confocal cavities, the sum over (l, m) is restricted to l+m being either odd or even [31]. Additional dispersion, present even at δL = 0, is due to mirror aberrations [32]. Ignoring aberrations for now, we may rewrite the interaction in Eq. 2 as where the spatial dependence of the interaction is encoded in the dimensionless function D(x, x ). As we discuss in more detail below, the restriction to either odd or even modes means that this function can be thought of predominantly as containing two contributions, a direct interaction D loc (x, x ) and its mirror image, D loc (x, −x ).
As a result, we will see that the quantity appearing in Fig. 2b, evaluated at equal positions x = x contains two contributions: a broad background of self interaction providing a flat plateau from the direct term, and a sharp peak from the mirror term for small values of x 1 . The range of the cavity mediated interactions can be extracted from the width of this peak.
For an ideal cavity, supporting an infinite number of modes, there would be a delta-function interaction peak from D loc (x, −x ) because the Hermite-Gaussian polynomials form a complete basis [33], and the background contribution from D loc (x, x ) would be entirely flat and nonzero. However, real cavities support only a finite number of modes, yielding a finite-range interaction: A photon is scattered into a wavepacket localized around the atom, and only atoms with overlapping polaritonic excitations-dressed atom-photon states-can interact.
We note that there is a mirror image of the cavity mode focused through −x 1 and, for two atomic gases, through both −x 1 and −x 2 . These virtual images arise because confocal cavities only support modes of good parity at each degenerate resonance. That is, the mode content alternates between all even or all odd modes every half free spectral range [31]. The fields at the real and virtual image locations are of the same (opposite) sign for cavities tuned to even (odd) modes resonances. We employ even mode configurations in this work. The direct and mirror contributions can be seen in Fig. 1a and b, which sketch the mode-supermode-that forms around each BEC for either one or two BECs in the cavity, respectively. Each supermode is the mixture of bare cavity modes by the atomic dielectric [9,34]. The minimum waist of the supermode is as small as the width of the atomic gas if there are sufficiently many intracavity modes to create a compact superposition.
The form of D(x, x ) is set by the parameter /∆ 0,0 , which may be experimentally controlled to tune the interaction range. The length scale of the range may be tuned between that of waist w 0 for a single-mode cavity to a small fraction of w 0 for a multimode cavity. This is analogous to the phonon-mediated interaction in ion traps, where large pump-detunings from resonances in the phonon spectrum generate shorter-ranged interactions [35][36][37].
To characterize the interaction profile D(x, x ), we use local measurements of the self-organization threshold for a small BEC. Using the expression for interactions in Eq. 4, we may write a mean-field threshold condition for self-organization as where ρ TF (x) is the Thomas-Fermi density distribution of the BEC. See Ref. [17,23] for the beyond-mean-field expression. For BEC radii much smaller than w 0 , we may approximate the density by ρ TF (x) ≈ δ(x − x 1x ) to obtain an expression for D(x 1 , x 1 ) at threshold, where Ω 0 is the threshold Rabi frequency for a deltafunction-width gas localized at the center of a singlemode cavity. Measuring Ω c (x 1 , x 2 ) for a pair of atomic gases located at x 1 and x 2 allows one to determine D c (x 1 , x 2 ) via the relation where we have dropped the subscript on the D's for convenience here and below. The value of Ω c (x 1 , x 2 ) at which the superradiant, self-organization transition occurs allows us to measure the photon-mediated atomatom interaction strength versus position through a closed form expression related to Eq. 9 described in Sec. IV. When considering a pair of gases, we may exploit the mirror symmetry to cleanly measure the interactions between atoms in different gases without physically merging two real BECs. That is, we use the fact that atoms in one gas can overlap and thus interact with the image of the other gas. By avoiding overlap of the real gases, we avoid unwanted systematics due to the change in atomic density and mean-field energy from collisions as the two traps merge. From the standpoint of photon-mediated atom-atom interactions, atoms at virtual images at −x i just like dipoles near a dielectric can be thought of as interacting with their mirror images in classical electrodynamics [38].
For x 1 = x 2 , and away from x 1 = 0 where x 1 approaches −x 1 , D(x 1 , x 1 ) assumes a nearly flat distribution. D(x 1 , x 1 ) begins to decay as a Lorentzian at a distance given by w 0 √ 2M * , where (M * ) 2 is a measure of the effective number of modes coupled to the atoms, see Sec. VI. This provides a translationally invariant interaction energy over a large distance in the xy-plane. For example, this distance is ∼200 µm on either side of our near-confocal cavity, far larger than typical BEC widths. We now describe the characterization of the strength and range of cavity-photon-mediated interactions for various pump and cavity configurations.

III. EXPERIMENTAL APPARATUS
We investigate the behavior of photon-mediated interactions by trapping within an adjustable-length multi- 4. Superradiant emission into the cavity supermode above the self-organization threshold. (a)-(i) Spatial structure of superradiant emission into the modes of a near-confocal cavity (δL = 0 µm) as the BEC is translated from x1 = −1.92w0 to 1.92w0. The two peaks merge at the center of the cavity, yielding a spot size smaller than the TEM0,0 mode waist. Images (j)-(p) show superradiant emission for BECs in a single-mode cavity. The BECs are in the same locations as in the above panels. The cavity is set to δL = 65.1 µm, ≈ 60 MHz to achieve near-single-mode operation. We see that the profiles are close to the width and shape of a TEM0,0 mode. All images are plotted with identical length scales, including panels a and i, which have larger fields of view. Data are taken at a pump-cavity detuning of ∆0,0 = 20 MHz. The white bars in (a) and (e) represent the length of the waist of the TEM0,0 mode. mode optical cavity a BEC of 2.5(3) × 10 5 87 Rb atoms in the |F = 1, m F = −1 state. See Ref. [32] for BEC preparation procedure and Fig. 1 for illustration of experiment. The BEC is confined in a crossed optical dipole trap (ODT) formed by a pair of 1064-nm laser beams propagating alongx andŷ with waists of 40 µm in the xy-plane and 80 µm alongẑ. The resulting trap frequencies of (ω (1)] µm that are significantly smaller than the w 0 = 35 µm waist of the TEM 0,0 cavity mode. Acousto-optic deflectors (AODs) placed in the path of each ODT control the intensity and location of the ODTs, allowing us to translate the BEC to any point in the xyplane with an uncertainty of 0.9 µm. Some of the experiments discussed require two intracavity BECs that can be moved relative to one another. We use dynamic trap shaping techniques [39] to split the BEC into two smaller BECs of 1.0(3) × 10 5 atoms each, with an imbalance uncertainty of <10%. These BECs may be separated in x by any relative distance using the AOD; see Fig. 1(b). Absorption imaging of the BECs is performed alongŷ after a 15-ms time-of-flight (TOF) to reveal the momentum distribution of either the single or double BECs.
The cavity is operated in a near-confocal regime in which the cavity length L is set to be close to the mirrors' R = 1 cm radius of curvature. Due to astigmatism, there are two orthogonal radii of curvature that are slightly different. Because L may only be set to match one radius at a time, the cavity is never perfectly confocal. This contributes, along with spherical aberration, to the finite bandwidth (small spread) of modes seen in Fig. 2(a) for δL = 0 µm [9,32]. The mode degeneracy is maximal when L = R, as shown in Fig. 2(a) for δL = L − R = 0 µm. A slip-stick piezo attached to one of the mirrors allows us to change L in situ [32]. The frequency spacing between each family of transverse modes is controlled by δL, which provides tunability of mode density; see the transmission spectra in Fig. 2(a). By family, we mean TEM l,m modes that satisfy l + m = const. We have observed modes in cavity transmission with indices up to l + m = 300. This indicates that up to ∼10 4 modes are supported by the cavity near degeneracy.
The system with an atom at the field maximum of the TEM 0,0 mode has a single-atom cooperativity of 2.2(1), a vacuum Rabi splitting of g 0 = 2π × 1.47(3) MHz, and κ = 2π × 167(4) kHz [32]. A laser propagating alonĝ x with Rabi frequency Ω pumps the BEC-cavity system near the even modes of the confocal cavity. The pumpcavity detuning ∆ 0,0 is defined as the difference in frequency between the pump and the TEM 0,0 mode. Where unclear, e.g., at small δL, the frequency of the TEM 0,0 mode is found by measuring the resonance frequency of a TEM 0,0 mode injected using a spatial light modulator [40]. The ∆ 0,0 's employed in this work are much larger than measured dispersive shifts at the atomic detuning of ∆ a = −102 GHz. To achieve homogeneous pumping of the BEC and to minimize any perturbation to the BEC trap potentials, the transverse pump has a large waist (1/e field radius) of 500 µm. In contrast to the standing-wave pump configuration used in previous studies of cavity-induced self-organization [26,28], we employ a running-wave pump [27] in the data taken in Figs. 2, 3, 5, and 7. This is done so as not to generate a lattice potential alongx in the absence of intracavity light. The absence of such a lattice leads to a simpler dependence of threshold pump power on the cavity mediated interaction: For a standing wave, one must calculate the kinetic energy for atoms in the band-structure of the standing-wave lattice potential, and this means that pump power would appear on both sides of Eq. (5), making extraction of interaction strength less direct. We do, however, use a standing-wave pump for the cavity output and atomic density images presented in Figs. 4, 5, and 6 to avoid distortions due to atomic motion excited by running-wave pump.
ϭ͘Ϭ Ϭ͘ϴ Ϭ͘ϲ Ϭ͘ϰ Ϭ͘Ϯ Ϭ͘Ϭ Ϭ͘Ϯ Ϭ͘ϰ Ϭ͘ϲ Ϭ͘ϴ ϭ͘Ϭ The local interactions versus position between a real BEC at x1 and the virtual BEC at −x2 of a different real BEC at x2. The data were taken in the confocal configuration (δL = 0 µm) with the two BECs located on opposite sides of the cavity; see

IV. MEASUREMENTS OF CAVITY-INDUCED INTERACTIONS
We first measure the interaction energy of a single BEC as a function of its location inx. With the BEC trapped at a location x 1 , the transverse pump power is linearly increased in time while the cavity emission is monitored on a single-photon counter. A sharp increase in emission heralds the onset of the superradiant, self-organization transition and allows us to measure Ω c (x 1 ), and consequently, using Eq. (6), extract the interaction strength D(x 1 , x 1 ). Figure 2a shows the transmission spectra of the even modes in the five near-confocal cavities studied. For large values of δL, and consequently large , the interaction strength would follow a single Gaussian decay as the BEC is moved further away from the center of the cavity. This is because D(x 1 , x 1 ) is following the mode profile of the TEM 0,0 mode in this near-single-mode cavity. By contrast, as is reduced by shrinking δL, we observe a form with two components. There is a flatter, more translationally-invariant background, falling off over a length scale w coming from the self interaction of the gas. Furthermore, we observe the emergence of a prominent peak at the cavity center that decays over a much shorter range ξ due to the interaction of the cloud with its mirror image. As becomes smaller, the scale w grows and ξ shrinks. A similar behavior is observed for holding fixed but varying ∆ 0,0 , as presented in Fig. 3.
The length scale w of the background component re-flects an overall envelope of the interactions, while the scale ξ of the sharp peak reflects the interaction range. Both the growth of w and shrinking of ξ can be understood as the result of superposing ever-larger high-order Hermite-Gaussian polynomials, allowing the interaction to both extend to larger distances and resolve finer features. That we can measure the short-range interaction with only a single, compactly localized BEC is a consequence of the mirror symmetry inherent to confocal cavities. The interaction energy increases as the real atoms come near to their virtual images, even though there is only one real BEC present. Viewed equivalently, as the two spots of the supermode begin to overlap, the intracavity field magnitude increases, leading to a lower Ω c .
Images of the supermode can be directly observed in superradiant cavity emission patterns. Figures 4a-i show examples in which the superradiant spots pass through each other. One cannot differentiate the spots from the picture alone, though from the orientation of the camera and apparatus, we know that the lower (upper) spot is the real image in panels a-d (f-i). The waists of the spots are smaller than that of a TEM 0,0 mode and their size at the object plane are similar to the BEC Thomas-Fermi radius, as expected. Their small size stands in stark contrast to the single-mode cavity's size shown in Fig. 4(j)-(p): The superradiant emission pattern maintains its TEM 0,0 structure as the BEC is moved over the same distance inx, only dimming as the gas nears the edge of the single Gaussian mode [41].
We now present similar measurements of two identical intracavity BECs. The BECs are located approximately 45 µm from either side of the confocal cavity center, at x 1 and x 2 as illustrated in Fig. 1(b). Each real BEC can then overlap and interact with the nearby virtual image of the other BEC. This is accomplished by moving x 1 and x 2 by the same amount inx while keeping x 1 − x 2 fixed. Again, this allows us to probe the behavior of the photon-mediated interaction while avoiding any energy shifts due to density changes and s-wave collisions between the two BECs. As shown in Figure 5, we observe four distinct spots in the superradiant emission pattern-two at the BEC locations (x 1 , x 2 ) and two at the locations of their virtual images (−x 1 , −x 2 ). At x 1 + x 2 = 0, the BECs merge with each other's virtual images, and we again observe a peak in D(x 1 , x 2 ) arising from a photon-mediated local interaction. The sequence of BEC momenta observed by time-of-flight imaging in Fig. 6 further demonstrates how the interaction energy of two nearby BECs can push a system above threshold: A smaller BEC can undergo self-organization at a lower threshold power when it is near a larger BEC than when it is far away. The dependence of the interaction range ξ/w0 = 1/ √ 2M * on ∆0,0 inferred from the data in panel (a). We measure an interaction range of ξ/w0 = 0.09(1) at the largest value of ∆0,0 studied in the confocal configuration. This is over an order-of-magnitude smaller than the TEM0,0 waist w0. Solid lines are the same fits to the theoretical interaction profile as above. The dashed lines are extensions of the fitted curve outside the regions of validity; i.e., where M * < 1. Inset: Agreement between interaction ranges extracted from the single-cloud (unfilled squares) and two-cloud (filled circles) datasets for δL = 0 µm. The solid line is a fit to the single-cloud data. All error bars represent standard error.
where D loc (x, x ) is a local interaction between two atoms and D loc (x, −x ) represents its corresponding atom-image interaction. The third term D non (x, x ) is a weaker, non-local oscillatory interaction which will be discussed later.
The local terms have the form where K 0 is the modified Bessel function of the second kind and falls off as an exponential for large δx = x − x , the separation between the atoms. The center-of-mass coordinate of the two atoms is X cm = (x + x )/2. The strength and range of this interaction are controlled by the parameter M * = ∆ 0,0 / . The quantity (M * ) 2 may loosely be associated with the effective number of cavity modes that maximally couple in 2D to the BEC to form the supermode (see Sec. VI D); M * is the number of modes that couple in 1D. We stress that the value of M * depends on the pump detuning and any aberration of the mirrors. Therefore one should not equate (M * ) 2 with the number of modes supported by the cavity near degeneracy, which is ∼10 4 . Examining Eq. 9, we note that the interaction range ξ = w 0 / √ 2M * decreases with increasing M * . A second length scale w = √ 2M * w 0 controls the strength of this interaction as the pair of atoms is moved far from the cavity center. Small M * dilutes the strength of this interaction versus distance from the cavity center: w may be interpreted as a measure of the degree to which the system is translationally symmetric. w diverges in an ideal confocal cavity as → 0, resulting in translationally-invariant interactions determined only by atomic separation δx, with no dependence on absolute position.
We characterize the range of the local interactions in our cavity by fitting our data in Fig. 2 and 3 to the theoretical model in Eq. 8, while neglecting the weak non-local term D non . To account for the finite size of the BEC, D(x 1 , x 1 ) is evaluated by numerically integrating over the BEC's Thomas-Fermi distribution ρ TF instead of a δ-function: We fit the above expression to our data using M * and an overall scale-factor as free fit parameters. Details of how this integral may be efficiently evaluated are given in Sec. VI C. The results of these fits are shown as solid lines in Figs. 2(b), 3(b), and 5. Extracted values of M * and the interaction range ξ are presented in Fig. 7(a) and (b) respectively, for several values of δL and ∆ 0,0 using the single BEC configuration. Large values of ∆ 0,0 / result in a more uniform coupling to transverse modes of the cavity and a shorter-ranged interaction. With this control parameter, we can tune the interaction range to be as low as ξ = 3.4(4) µm. This is over an order-ofmagnitude shorter than the range set by w 0 for a singlemode cavity. Moreover, this close agreement between the data and fits for values of M * 1 highlights the validity of the theoretical model presented in Sec. VI. We note that we do not reliably infer M * for M * < 1 because the closed-form expression in Eq. 9 is only valid for ∆ 0,0 . We now turn our attention to the non-local interaction term in Eq. 8, which displays oscillatory behavior of the form As discussed below, the form of this term can be traced to the Gouy phase shifts of the bare-cavity modes [31]. While we cannot resolve the effects of this term in our interaction measurements, we do observe a weak signal in our images of superradiant cavity emission shown in Fig. 8. The cavity emission is recorded by imaging the plane containing the atoms onto the camera, so this emission records the light profile at the atom plane. For a single BEC at x 1 , the image at position x corresponds to D(x, x 1 ), and so the non-local term creates fringes in the cavity emission with a wavelength that becomes shorter as x 1 is increased, as shown in Fig. 8b-d. The oscillatory behavior can most easily be understood by considering the "hourglass" structure of confocal cavity modes. While familiar ray-tracing representations of these modes depict the parallel and diagonal "arms" of the closed hourglass path [31], they do not account for interference of the paths. A full calculation of the field of a confocal cavity supermode is shown in Fig. 8a: The parallel arms of the hourglass path create two spots at x 1 and −x 1 , while the diagonal arms interfere with each other to create fringes along x. A calculation of the superradiant intracavity field pattern shown in Fig. 8e-g, using the theory presented in the next section, reveals a similar structure and is in qualitative agreement with our data.

A. Hamiltonian and equations of motion
To derive the atom-atom interaction, we start from a model of N atoms in a condensate wavefunction Ψ(r) interacting with cavity modesâ µ by the Hamiltonian (12) where for compactness we use µ = (l, m) to index the transverse modes of our cavity. The first term is the Hamiltonian of the cavity modes with detuning ∆ µ = ∆ 0,0 + (l + m) . The remainder is the standard Hamiltonian for a weakly interacting BEC with contact interactions of strength U in an external trap V (r) and with a Stark shift proportional to 1 ∆a due to the light in the cavity. This light fieldφ consists of the running-wave pump and a sum over all cavity modes with their transverse and longitudinal spatial dependencê (13) where Ω is the pump Rabi frequency, Ξ µ (r) is a Hermite-Gauss mode of the cavity and θ µ contains other contributions to the phase which vary slowly compared to kz. In particular, its dependence on µ = (l, m) is due to the Gouy phase term (l + m)[π/4 + arctan(z/z R )], where z R = L/2 is the Rayleigh range, and this formula assumes z is measured from the center of the cavity. This term accounts for the fact that in order to have equal frequencies, a mode with higher order transverse structure must have a slower rate of change of longitudinal phase [31]. The form of Eq. (13) results in a spatially varying single-photon Rabi frequency g 0 Ξ µ (r)/Ξ 00 (0) for the mode µ [42].
To study the location of threshold, we assume that most of the condensate is in the ground state, with a small fraction having a momentum kick from either scattering a photon from the pump into the cavity or viceversa. Hence we write where Z is an envelope function which describes the confinement of the gas inẑ, ψ 0 (r) is the wavefunction of the ground state of the gas in the transverse plane, and ψ F (B) is the wavefunction of the gas which has been scattered forward (backward) by scattering between the pump beam and the cavity modes. Due to scattering into the cavity modes, these functions ψ F,B have a sinusoidal variation kz along the cavity. However, because of the different Gouy phase terms of different cavity modes, it is not a priori clear what phase the atomic density wave should take. To allow the possibility of coupling to any cavity mode we further decompose the scattered atomic wavefunctions into two out-of-phase density waves (15) and similarly for ψ B . Here, x = (x, y) is the transverse coordinate vector, and ψ F (1,2) (x) are now slowly varying envelope functions. As we see below, different cavity modes couple preferentially to ψ F 1 or ψ F 2 . For convenience, the phase offset θ 0,0 (z 0 ) corresponding to the Gouy phase of the (0, 0) mode at the position of the atomic gas is introduced. We can now use Eq. 12 to find the mean-field equations of motion for ψ 0,F,B and α µ ≡ â µ . As the threshold is where the normal state α µ , ψ F,B = 0 becomes unstable, we need only do this to leading order in these fields.
To write equations in terms of only the transverse coordinates x we must perform the z integral in Eq. 12. This can be done straightforwardly in the limit where we assume Z(z − z 0 ) has a width σ z and that λ σ z z R . The first inequality allows us to drop any terms oscillating at wavevector k; this imposes momentum conservation so that recoiling atoms pick up the difference of pump and cavity momenta. The second condition means we can evaluate the slowly varying phase terms as being effectively constant over the width of the gas: we can approximate θ µ (z) θ µ (z 0 ) ≡ θ 0,0 (z 0 ) + (l + m)θ 0 with θ 0 ≡ π/4 + arctan(z/z R ). In the linearized treatment, all relevant z integrals involve the cross pump-cavity term causing scattering between at-rest atoms ψ 0 (x) and the functions ψ F,B (r). We then find the z integrals yield two possible values, For the equations of motion, we find where we have included photon loss κ and ω r is the recoil momentum k 2 /2m. The ground state condensate has no linear perturbations, so at leading order we have: while ψ Bi (x) obeys an identical equation to Eq. 18 with F ↔ B and α µ → α * µ . The ground state density profile is that of a Thomas-Fermi gas ρ(x) = ρ 0 [1 − (x/x 0 ) 2 − (y/y 0 ) 2 ], and so we have taken ψ 0 (x) to be real.

B. Calculation of effective interaction D(x, x )
We wish to study the effective photon-mediated atomatom interaction. Since we expect the cavity field to reach a steady state on a timescale much faster than the atomic motion, we adiabatically eliminate the photons by setting the time derivative in Eq. 17 to zero and solving for α µ . We also neglect the corrections to the bare cavity modes caused by the ground state atomic gas; i.e., the term proportional to the integral of |ψ 0 (x)| 2 in Eq. 17 is set to zero. Substituting this back into the equation of motion of the atomic condensate gives where we defined an atomic Hamiltonian, and the cavity-mediated interaction takes the form: To simplify further, we assume that the atoms are close enough to the cavity center that θ(z 0 ) ≈ π/4. In this case one may see that as long as l + m is even, either O 1 µ = 0 or O 2 µ = 0, so the interaction becomes diagonal, x ) . Furthermore, using standard trigonometric identities we can reduce the expression to: In writing this, we introduced the factor 1 + (−1) l+m into S l,m so that the sum in Eq. (23) is now over all modes. This extra factor serves to cancel odd modes. This rewriting will enable us below to make use of known expressions for sums of Gauss-Hermite functions multiplied by phase factors, exp[iϕ(l + m)]. As a reminder, the detuning in the denominator takes the form ∆ l,m = ∆ 0,0 + (l + m) This term D i (x, x ) is the expression given in Eq. (4), the interaction between atoms at different points x and x due to the cavity modes (except that in Eq. (4) we neglected cavity loss). Again, an identical equation to Eq. 20 holds for ψ B (x) with F ↔ B and D(x, x ) replaced with its complex conjugate.

C. Analytic forms of interaction near confocality
In this section, we discuss those cases where it is possible to extract an analytic closed form for the interaction term, Eq. (23). We are able to find expressions for both an ideally confocal system ( = 0) and a near-confocal cavity with = 0. Moreover, we show that restricting the number of modes contributing to the interaction and including deviations from confocality affect the interaction profile similarly. This connection allows us to identify, in Sec. VI D, an effective number of modes M * that couple to the atoms.
If we first consider the ideal confocal case, = 0, the denominator in Eq. (23) becomes a constant, independent of l, m. In this case, we can make use of the har-monic oscillator Green's function, In terms of this, the interaction can be written as whereκ = κ/∆ 0,0 and we have made use of the relation G(x, x , α − iπ) = G(x, −x , α). If we then take the limit of Eq. 26 for α → 0, we find the simple expression consisting of a local interaction between atoms, a local interaction between atoms and virtual atoms at their mirror image, and a non-translationally invariant oscillatory interaction.
We can extend this result at confocality to find the interaction function in the limit of near confocality, where ∆ 0,0 . Defining˜ = /∆ 0,0 we rewrite the l, m dependence of the denominator as an integral: We may group the terms together as discussed in Eq. (8) to write The non-local contribution comes from the last two terms in Eq. (28). By using the identities sinh(θ−iπ/2) = −i cosh(θ), cosh(θ − iπ/2) = −i sinh(θ), we can write: Because the first exponential suppresses contributions where τ 1, we may consider the small˜ behavior by making a small˜ τ expansion, tanh(˜ τ ) ˜ τ and cosh(˜ τ ) 1 along with 1 + e −2˜ τ 2. The τ integral then becomes straightforward, yielding: For the local terms, a similar expansion for small˜ τ is possible, however here we must note the prefactor involves 1 − e −2˜ τ 2˜ τ . We thus find: The τ integral here can be shown to produce a modified Bessel function of the second kind, i.e.
Because the Bessel function diverges at zero argument, it is crucial to consider the smoothed version of this function when comparing D eff (x 1 , x 1 ) in Eq. (10) to experimental results. In doing this, we may note that the two terms in the argument of the Bessel function have very different dependence on coordinates. The first term depends strongly on the separation, with a characteristic length scale w 0 ˜ /2, while the second term (inside the square root) has a much weaker dependence, with a characteristic length scale w 0 2/˜ w 0 . In the smoothed function D eff (x 1 , x 1 ), we integrate over Thomas-Fermi distributions of the atom cloud. Assuming the cloud width is small compared to the length scale w 0 2/˜ , we may neglect any difference between x, x and x 1 when evaluating the term in the square root. This leads to the expression: Assuming symmetric clouds, this can further be simplified by suitable changes of variables to put it into the form of a convolution, where ρ 2 (z) = dyρ TF (y)ρ TF (z − y). This is the procedure used in fitting Figs. 2 and 3. Note that the first Bessel function in Eq. (33) describes the "self" interaction of the cloud, and its only dependence on x 1 is via the square root in the Bessel function, which ultimately leads to a slow fall off with length scale w ≡ w 0 2/˜ w 0 . The second term is the mirror interaction, and falls of exponentially with x 1 with a length scale ξ ≡ w 0 ˜ /2 w 0 . To see this behavior more clearly, we can consider the analytic expressions that result if we replace ρ TF (x) by a Gaussian of width σ. In this case ρ 2 (z) is a Gaussian with width √ 2σ. For the first term, which we denote D eff,self (x 1 ), we may use the result: which comes from an integral representation of the Bessel function and defining A = (1/2ξ) 1 + iκ + x 2 1 /w 2 . In this expression, the quantity Aσ > (σ/2ξ) 1, and thus for the typical values of τ that dominate the integral, we have 4A 2 σ 2 τ . We thus find the first part of Eq. (33) has the form: There is no such simple closed form for the image term. However, using the same approach as above we can write the expression in the form where now A = (1/2ξ) √ 1 + iκ. We still have that Aσ 1, however the extra terms in the exponent means it is no longer always true that the integral is dominated by terms for which τ 1. At large x 1 , the saddle point of the integral occurs when τ 2Ax 1 , and so for large enough x 1 we have that the dominant contribution comes from values for which τ A 2 σ 2 . The crossover occurs when x 1 σ 2 /ξ. We thus have two asymptotic limits: D. Relating the ratio of mode dispersion to mode detuning˜ to the effective number of coupled modes (M * ) 2 In order to make precise the sense in which we regard a non-zero˜ = /∆ 0,0 as corresponding to a finite mode cutoff, we discuss here the results for such a cutoff. For simplicity we consider a "square" cutoff, where we remove all modes Ξ l,m (x) with either l, m > M . This means we may write expressions in terms of the 1D Green's functions. Neglecting non-local terms we have The 1D Green's functions with finite cutoff can be written in terms of the Christoffel-Darboux identity to give: Here and throughout this section we measure all lengths in units of the cavity beam waist, i.e., w 0 ≡ 1. Using the 1D Green's functions we want to evaluate: In the following, we will find approximate forms for these terms at large M . For simplicity, we assume M is even; similar results occur for odd M , but with various sign changes in intermediate formulae. To consider the behavior at large M , we make use of the Wenzel-Kramers-Brillouin (WKB) approximation for a Gauss-Hermite function: The mirror term has a simple form as we may write The 1/x factor means we need only focus on behavior at small x. This means we can approximate the phase function S M (x) √ 2M + 1x and the envelope function as E M (x) 1/ 4 √ 2M + 1, and so we find Thus, we find that the mirror term describes a sharp peak with a width that scales as 1/ √ 2M + 1. Comparing this to the results at non-zero , we can identify an "effective" mode number, M * = ∆ 0,0 / , as parameterizing this finite peak width.
The self interaction term is more subtle. We can first rewrite the Green's function at x → x in terms of derivatives and then use the recurrence relation on Gauss-Hermite functions: One may now use that for large M , we can neglect differences between the envelope functions, E M (x) E M ±1 (x), and approximate √ 2M 2(M + 1) √ 2M + 1 in the prefactors to write: If we consider that: one may readily see that S M +1 (x)+S M −1 (x) = 2S M (x)+ O(M −2 ), and so to leading order in 1/M we have: Finally, using double angle formulae gives the result This shows the self interaction term gives a broad semicircular function. Its algebraic form does not match the finite˜ result, but we can again identify the width of this function, √ 2M + 1 with the width of the non-zero self interaction, to again give the identification M * = ∆ 0,0 / .

VII. CONCLUDING DISCUSSION
The ability to engineer tunable-range interactions among intracavity atoms and, equivalently, high-M * systems, opens several research directions. We conclude with a discussion of three such directions, one involving exotic spatial organization of superfluid atoms and two involving spin organization.
The nature of the superradiant, self-organization phase transition can differ in a multimode cavity. In the singlemode cavity, the behavior is governed by mean-field theory, due to the all-to-all coupling. In contrast, the multimode cavity allows transverse variations of phase across the cavity. The atomic gases studied in this paper are (purposefully) too small to allow such variations, but by combining much larger intracavity BECs with the confocal cavity, transverse phase variation becomes possible. This has a number of consequences.
An immediate consequence of transverse phase variation is the possibility of topological defects and phase textures. This is because atoms are no longer constrained to organize with respect to the profile of a single mode, but may fluctuate between the Hermite-Gaussian profiles of the multiple degenerate modes. Similar to classical systems like diblock copolymers and fluids undergoing Rayleigh-Bénard convection, the organization should exhibit wandering stripe-like (smectic) patterns of atoms [17,23]. The interaction length scale ξ controls the minimum size of a patch of stripes pointing in the same direction with the same spatial period, while the envelope of the interactions w controls the maximum size of an atomic gas that can fully couple to the cavity. As a result, the number of patches in the 2D transverse profile is ∼(M * ) 2 . The fact that we can engineer systems with M * 1 means that such complex, superfluid smectic states are within reach [17,23]. This opens the door to exploring analogs of the quantum liquid crystals found in strongly correlated electronic materials, such as cuprate and iron-based high-T c superconductors [43]. Controllability of quenched disorder and dimensionality using external optical dipole trap beams and speckle would provide unique ways to investigate the intertwined nature of the order-crystalline, superfluid, and even magnetic (see below)-found in these systems [44].
A second consequence of transverse degrees of freedom is their effect on the universality class of the phase transition. For a single-mode cavity, the all-to-all coupling means the phase transition-analogous to the Hepp-Lieb-Dicke transition [4]-falls within the mean-field-Ising universality class. This second-order mean-field phase transition is expected to become weakly first-order as the number of degenerate modes increases [17,23]. This occurs in a scenario akin to that of a quantum version of the Brazovskii transition known from classical liquid crystal physics [45,46]. As one approaches the critical pump strength for the second-order transition, soft modes emerge corresponding to long-wavelength transverse fluctuations. The additional, beyond-mean-field contribution to the effective action that arises from these fluctuations drives the transition first order.
In addition to modifying the universality class of the phase transition, the presence of soft transverse modes can be seen in other ways. The dispersion relation of such modes could be measured through established methods for observing dynamical susceptibilities [47]. Just like xray diffraction patterns of classical liquid crystals are arcshaped [48], signatures of this quantum liquid crystalline state might appear as arc-like Bragg diffraction peaks in time-of-flight measurements. Because of the small size of the atomic gases, no such patterns are seen in Fig. 6, but may become apparent by expanding the size of the intracavity BEC. This may easily be accomplished by lowering the optical dipole trap frequencies.
In the current configuration, the cavity mediates interactions between atomic density-wave excitations. One can also consider cavity-mediated interactions between atomic spins. These can be engineered if the transverse pump lasers drive a Raman transition between atomic Zeeman states representing a pseudospin-1/2 system [13,49,50]. If the atoms are trapped at random positions inside the cavity to realize quenched disorder, then the multiple modes of the cavity can in principle mediate frustrated spinful interactions resulting in a spin glasslike state [13,14]. However, there is a subtlety regarding the effects of summing many cavity modes: in some geometries, the sum over cavity modes may yield a short range interaction, in which case the degenerate limit produces a short-range spin model. However, as we have shown in this paper, the Gouy phase naturally present for a confocal cavity also induces a long-range sign-changing interaction D non (x, x ) ∼ cos(x·x /w 2 0 ). Such an RKKYlike sign-changing interaction is exactly the ingredient needed to enable glassy physics [12]. The ability to tune the relative strengths between this long-range interaction and the short-range interaction D loc (x, x ) provides a unique means (outside of numerical simulation) to ex-perimentally compare the dynamics of infinite-range spin glasses to those with short-range interactions. While the former has an order known to be described by mean-field replica-symmetry breaking, the latter's order defies explication despite many decades of investigation [51]. Direct spin-state detection combined with repeatable atomic disorder from shot-to-shot will allow us to create, observe, and compare system replicas. This may provide sufficient experimental information to discriminate among various theories of short-range spin glass order.
Spin glasses may serve as models for neural networks. Realizing spin glasses would provide the means to create a neural network comprised of atomic spins serving as neurons, cavity modes serving as synapses, and photons within the modes serving as action potentials [13,52]. Wiring the network to implement a particular graphical combinatorial optimization problem simply involves placing the atoms in specific locations within the cavity modes. This may be possible with optical tweezer arrays [53,54]. The combination of local and non-local interactions demonstrated here will enable the construction of a wide variety of graphical combinatorial optimization problems, not just those of a complete graph. In this way, Hopfield associative memories [13,52,55,56] and coherent Ising machines [57,58] may be implemented in the presence of quantum effects like spin entanglement and quantum criticality, providing a new route to quantum neuromorphic computation.