A rapidly expanding Bose-Einstein condensate: an expanding universe in the lab

We study the dynamics of a supersonically expanding ring-shaped Bose-Einstein condensate both experimentally and theoretically. The expansion redshifts long-wavelength excitations, as in an expanding universe. After expansion, energy in the radial mode leads to the production of bulk topological excitations -- solitons and vortices -- driving the production of a large number of azimuthal phonons and, at late times, causing stochastic persistent currents. These complex nonlinear dynamics, fueled by the energy stored coherently in one mode, are reminiscent of a type of"preheating"that may have taken place at the end of inflation.

Cosmological expansion is central to our understanding of the universe. Here, we experimentally create a system where fields expand in a similar way as in the universe: an expanding, ring-shaped atomic Bose-Einstein condensate (BEC). Our laboratory test bed demonstrates three effects associated with the expanding universe. First, we conclusively demonstrate a redshifting of phonons analogous to the redshifting of photons, which provided the first evidence for an expanding universe [1]. Second, we observe hints of "Hubble friction" that damps the redshifted fields [2]. Third, we observe a process in which energy is rapidly transferred from a homogeneous radial mode into azimuthal modes by a nonlinear, turbulent cascade, reminiscent of that seen in some models of preheating at the end of cosmological inflation [3][4][5][6][7][8][9]. Experiments such as these can thus emulate both linear and nonlinear field theoretic aspects of cosmology.
A zero-temperature BEC is a vacuum for phonons [10], just as an empty universe is a vacuum for quantum fields, like photons. In this analogy, the speed of light is replaced by the speed of sound, c, in the BEC. Evoking general relativity, the equation for long-wavelength phonons (in the hydrodynamic limit [11]) takes on a covariant form with a curved spacetime metric [12][13][14]. Previous studies with ultra-cold atoms illuminated different aspects of this phonon metric. For example, an interface between regions of sub-sonic and super-sonic fluid flow forms a "sonic event horizon" that exhibits effects such as Hawking radiation [12,[14][15][16][17][18][19][20][21][22]. By changing the interaction strength or density, one can simulate cosmological phenomena such as pair production [23][24][25][26][27][28][29], Sakharov oscillations [30], or the dynamical Casimir effect [31], the latter two having been recently observed experimentally [32,33]. Beyond cold atoms, experimental studies have realized analog event horizons in other settings, for example in optical systems [34][35][36] and in classical fluids [37][38][39]. (For a recent review, see Ref. [40].) The expansion of our BEC-universe is forced by dynamically increasing the radius of our nearly-flat bottomed ringshaped potential [41], as opposed to being governed by an analog of the Einstein equation (see Appendix A). Figure 1 shows our BEC during a t exp = 23.1 ms inflation. The radial velocity of the trapping potential (defined as the rate of change of the mean radius, R) is directly controlled, and can be made comparable to the speed of sound. For the expansion shown in Fig. 1, the maximum velocity is v p = dR/dt ≈ 1.3c, implying that points separated by an angle π/4 recede faster than c. The condensates used in this work are well-described by mean field theory; thus, we compare our measurements to numerical simulations using the stochastic-projected Gross-Pitaevskii equation (SPGPE, see Appendix B), which accurately captures BEC dynamics with thermal fluctuations [42,43]. Images from this simulation are in excellent agreement with the corresponding experimental images. Phonon redshift -To study the red shifting of phonons, we first imprint a standing wave phonon excitation on the background BEC. During expansion, these effectively onedimensional azimuthal phonons are redshifted, i.e., their wavelength grows as shown in Fig. 2a for both experiment and theory. These images show the oscillation of a standingwave phonon, constructed by perturbing the condensate with a potential of the form sin(mθ), where m is the integer azimuthal mode number of the phonon. The (approximate) axisymmetry implies that m is conserved, in analogy with conservation of the comoving wavevector in cosmology. The phonon wavelength is therefore stretched by a factor a = R f /R i , the ratio of the geometrical radii of the expanding ring. This is related to the usual redshift parameter z through a = z + 1. Figure 2b shows the measured phonon amplitude δn vs. time for various a and m and clearly shows a shift in the frequency. (In this paper, we measure frequency and time in the laboratory frame, as opposed to using the comoving proper time as defined by the effective metric, Eq. C7.) To measure the frequency shift, we fit the oscillation before and after expansion to extract ω m,i /ω m, f , shown in Fig. 2c. At any given time, the phonon oscillation frequency is ω(t) = c θ (t)m/R(t), where c θ (t) is the azimuthal speed of sound at time t. As the ring expands, both the atomic density and c θ decrease. For the combination of harmonic confinement in the vertical direction and roughly quartic confinement in the radial direction, we find c θ ∝ R −2/7 . The solid, black curve shows the resulting ω m,i /ω m, f = a 9/7 scaling; a full Bogoliubov calculation, with the azimuthally averaged potential, is shown as the solid, colored curves.
We understand the phonon's behavior during the expansion epoch in terms of a 1D equation for the phonon amplitude χ m , where δn = ( /U 0 )∂χ m /∂t, U 0 = 4π 2 a s /M, a s is the s-wave scattering length, and M is the mass of an atom. (See Appendix C for the derivation.) There are two contributions to the damping of the amplitude. The first damping term, γ m , is phenomenological, but independently measured [45]. The second,Ṙ/R, is analogous to the "Hubble friction" in cosmology, which damps fields with frequencies ω ȧ/a. In the present case, the Hubble friction has the largest impact when for supersonic expansion, i.e., when ω Ṙ /R or mc θ Ṙ .
For our expansions, we expect that the Hubble friction will play a role, particularly for the a = 4.1, m = 1 expansion whereṘ/R 1.5ω. (At maximum velocity,Ṙ/R ≈ 3γ m for m = 2 andṘ/R 20γ m for m = 1, but this occurs only during the short expansion epoch.) The Hubble friction term changes the phase and amplitude of the phonon oscillation after expansion. However, because the observed density difference δn is proportional to ∂χ m /∂t (see Appendix C), the predominant difference in observed amplitude before and after expansion results from the change in ω. To search for the Hubble friction term, we fit all the data simultaneously to Eq. 1, takinġ R/R → γ HṘ /R, where γ H is a tunable parameter. While the best-fit value γ H = 0.55 (21) indicates the presence of Hubble friction, the deviation from unity suggests that other effects like azimuthal asymmetry and non-zero annular thickness also affect the phonon amplitude [44]. For GPE simulations of the expansion of an azimuthally symmetric, thin annulus ring with a potential of a similar functional form, Eq. 1 is an accurate description of the phonon evolution (see Appendix C).
Radial dynamics -The preceding 1D discussion (based on Eq. 1), rested on the assumption that the background BEC contained no transverse dynamics. Perhaps the first indication of additional dynamics is visible in evolution of the ring-BEC's radius, shown by the red symbols in Fig. 3. As indicated by the oscillations around the trap's mean radius (black curves), the BEC is excited after the potential has reached its final value. The amplitude of the oscillation can be estimated based on a simple harmonic oscillator model, where the oscillator is the first radial phonon mode and forces applied are due to the expansion of the confining potential. These oscillations decay rapidly, typically within a few oscillation periods for all scale factors and expansion velocities studied. If the trap were perfectly harmonic, this center-of-mass oscillation should be long-lived. In reality, our trap is more flat-bottomed, is anharmonic, and is not axially symmetric.
To understand this rapid decay, we show the atomic density and phase of a simulated Bose-Einstein condensate without an imprinted phonon during the first few oscillations after expansion in Fig. 4a and b, respectively. At t = 10 ms, the condensate reaches the far end of the potential and begins to turn around. At t = 11.5 ms, the condensate phase is approximately flat radially, with the exception of a discontinuity of ≈ π in the center of the annulus. This standing wave has nodes in the atomic density with corresponding π phase jumps, effectively imprinting a dark soliton onto the BEC [46,47]. This process is analogous to the creation of solitons upon Bragg reflection in an optical lattice [48] or reflection of a condensate off of a tunnel barrier [49]. (Unfortunately, due to imaging limitations, we are unable to resolve solitons or other similarly-sized structures in the experiment.) The number of solitons N s created from the decay of the radial mode can be estimated by comparing the energy per particles contained the radial excitation to the energy per particle of a soliton ( s ≈ 4 c/3R T , where R T is the annular width of the ring). The amplitude of the radial excitation χ r , while calculable analytically, is a complicated function that depends exponentially on the adiabaticity of the expansion relative to the frequency of the radial mode ω r . (Assuming a box-like potential in the radial direction implies ω r ≈ πc/R T .) The adiabatic condition then demands, in our system,Ṙ must nearly be supersonic to produce solitons, i.e., v p 0.8c [50].
Turbulence and reheating-Dark solitons are unstable in condensates of more than one dimension. They suffer from a "snake instability" causing the soliton to first undulate and then fragment into vortex dipoles [51][52][53]. As shown by our numerics in Fig. 4a, the undulation is underway by 12.5 ms and the fragmentation into vortices is mostly complete by 14 ms. Theoretical estimates for a single soliton in a harmonically confined BEC suggest that the snake instability will result in N vd,1 ≈ 2πR/8ξ vortex pairs, where ξ = 2 /2Mµ is the local healing length within the bulk of the condensate and µ is the chemical potential [54]. For the present case, this corresponds to N vd,1 ≈ 50 vortex pairs. At t = 13 ms in Fig. 4b, the single soliton has decayed into ≈ 6 pairs over an angle ≈ 45 • near the top of the ring. This corresponds to roughly 48 vortex pairs around the full ring. These vortex pairs then form a highly turbulent state.
We experimentally observed the fingerprints of this process through the structure factor S (k θ ), a measure of the spatially structured density fluctuations (i.e., azimuthal phonons) with wavevector k θ = m/R. For both experiment and theory we extracted S (k θ ) by first evaluating the one-dimensional density n 1D (θ) around the ring to obtain the density fluctuations δn 1D (θ) = n 1D (θ) − n 1D (θ) , where · · · denotes the average over many realizations. The structure factor is Theoretical structure factors are shown in the top row of Fig. 4c; experimental structure factors are shown in the bottom row Fig. 4c. The colors in Fig. 4c identify the times at which the structure factors were evaluated. The density obtained from experiment has limited spatial resolution, is impacted by imaging aberrations, and has additional noise from the partial transfer absorption imaging process [55]. For these reasons, we first corrected for imaging aberrations (see Appendix A) and identified the detection threshold (shown by the horizontal lines). We used the numerical simulations (which include the same aberrations) to verify the correspondence between the corrected value of S (k θ ) based on simulated imaging to the value of S (k θ ) calculated from the simulated atomic density. These agree for values of S (k θ ) above the detection threshold.
As shown by the S (k θ ) snapshots, the structure factor starts at our detection threshold [56]. After expansion and during the soliton's initial formation (t = 12 ms), S (k θ ) maintains this value, indicating that this state does not differ significantly between realizations. When the soliton begins to break apart at t = 13 ms, a small peak, still below our detection threshold, appears in the simulations near k ≈ 1 µm −1 (not shown). This corresponds roughly to the wavenumber of the snake instability, k ≈ 2π/8ξ ≈ 1.3 µm −1 . As the turbulent state develops, this peak grows and shifts to lower k, becoming detectable at 18.5 ms and becoming larger at 22.5 ms. The shift to lower k θ is expected because of the inverse cascade that occurs in two-dimensional turbulence [57].
Stochastic persistent currents -While most of the vortex dipoles recombine and produce lower energy phonons, some of the vortex dipoles manage to break apart and become free vortices. If one of the free vortices slides into the center of the ring and one leaves the outside of the ring, then the overall phase of the ring slips by 2π and the winding number , quantifying the persistent current state of the ring, changes by one [58]. Indeed, we observe stochastic persistent currents in the ring after expansion in both the experiment and simulation. Figure 5a shows the resulting distributions of winding numbers for various speeds of expansion for a = 1.4(1) Evidence for this process can be found by studying the width of the winding number distributions for expansions with different a and t exp . The number of vortex dipoles produced from N s solitons would be N vd ≈ N s (2πR f /8ξ). The measured distribution widths collapse reasonably well when plotted versus N 1/4 vd , as shown in Fig. 5b. The 1/4 may result from some combination of the stochastic nature of dipole dissociation and recombination, the interaction-driven dynamics of dipoles and free vortices in a turbulent fluid, and the random phase-slip process.
One may question whether the appearance of the winding number might involve another cosmological phenomenon: the presence of sonic horizons. If we assume the speed of sound sets a limit on the speed at which information can travel through the condensate, the rapid supersonic expansion should create regions of condensate that are causally disconnected. The typical horizon distance established during the expansion would be given by, where c 0 is the initial speed of sound [59]. This leads to N R 2πR f /R hor ≈ πR f /c 0 t exp disconnected regions. If these regions' phases evolve at different rates and become sufficiently randomized, then when the regions recombine, they can form a topological excitation in the form of a persistent current [60,61]. The probability for a given persistent current is then given by the geodesic rule [61][62][63]. The red Gaussians in Fig. 5a show the expected distributions resulting from this horizon model, which disagree with the experiment. Moreover, simple estimates for the phase fluctuations present in our condensate are a factor of 25 too low to sufficiently randomize the phase during expansion. Future studies using condensates of lower density could see this effect, as the phase fluctuations will be larger. Discussion and Outlook -In this work, we explored the physics of a rapidly expanding Bose-Einstein condensate. We observed the redshifting of phonons during this rapid expansion, which has clear analogs in cosmological physics. After expansion stops, the condensate reheats through the creation and subsequent destruction of dark solitons, producing a highly turbulent state. This process leads to the creation of global topological defects (i.e., persistent currents), which at first might be thought to arise due to the presence of cosmological horizons, but actually result from the vortices produced when the solitons break apart.
While we see evidence for Hubble friction in our system, future studies should be able to more precisely measure its influence during the expansion of the phonon modes. In particular, by varyingṘ/R, one could more easily distinguish between the Hubble friction and other damping effects. One could also contract the ring rather than expand it. Because the Hubble friction is not dissipative and is reversible, such a contraction should cause amplification of the phonon mode amplitude.
The process of expansion, which presumably cools the azimuthal degrees of freedom of the condensate, followed by the increase in azimuthal excitations (Fig. 4c-d) as the radial mode decays, is reminiscent of the reheating process in the early universe. At the end of inflation in the universe, the energy contained in the homogeneous mode of the quantum field that drove inflation, the inflaton, decayed into inhomogeneous excitations. It is not known how this occurred. In the simplest model, the inflaton oscillated around the minimum of its potential, decaying into lower energy particles [3], whereby the radial mode couples directly to lower energy azimuthal phonon modes. However, the decay of the radial mode through this process is expected to be much slower (≈ 1 s −1 , using a calculation similar to that found in Ref. [64]) compared to the observed decay of the radial mode through soliton and vortex excitations (≈ 100 s −1 ). Future studies using a ring with stronger radial confinement should suppress the non-linear excitations and enhance the direct coupling. Other models are non-perturbative and include selfinteractions in the inflaton field that can lead to turbulent cascading [4][5][6][7][8][9], much like the turbulence we observe here.
Perhaps surprisingly, the long-wavelength azimuthal phonon mode is redshifted in simple way (Fig. 2), despite the complex dynamics occurring in the underlying BEC state. This survival has a direct analogy in inflationary cosmology. During inflation, vacuum fluctuations were redshifted to large length scales and amplified. The subsequent preheating and thermalization processes took place on shorter length scales, yet the resulting thermal state was modulated by the longwavelength amplified vacuum fluctuations. This process gave rise to the large-scale structure we observe today in the universe.
In addition to the possibilities described above, we anticipate that with new developments, other interesting cosmological phenomena might be realized with expanding condensates. First, with improved imaging that captures the initial (quantum and/or thermal) fluctuations, one could observe effects related to the scaling of the vacuum. In particular, one could observe cosmological particle production [23][24][25][26][27][28][29]. Second, a ring with stronger radial confinement will suppress transverse excitations, revealing the physics arising from the recombination of causally-disconnected regions. Given these possibilities, we believe an expanding ring BEC could provide an interesting laboratory test bed for cosmological physics.

ACKNOWLEDGMENTS
The authors thank J. Ho for initial discussions and a careful reading of the manuscript. We thank W.D. Phillips, E. Goldschmidt, M. Edwards and N. Proukakis for useful discussions. We thank the anonymous referees, whose comments greatly improved the manuscript. This work was partially supported by ONR and the NSF through the PFC at the JQI. TJ was supported in part by NSF grants PHY-1407744 and PHY-1708139. IBS was partially supported by the AFOSR's Quantum Matter MURI and NIST.

Appendix A: Experimental Details
Our experimental setup consists of a BEC of 23 Na atoms in an optical dipole trap (ODT). Our BECs are created using standard laser cooling techniques, followed by evaporation in first magnetic then optical dipole traps. In this experiment, we work with BECs with between 1 × 10 5 and 4 × 10 5 atoms. For measurement, we use partial transfer absorption imaging (PTAI) [55].
The final stage of evaporation begins when thermal atoms are loaded into a combination of the vertical trap and dimple trap. Vertical confinement is created using a blue-detuned (532 nm), TEM 01 beam, tightly focused to create two parallel sheets of light with a dark region in between. The dimple trap is a red-detuned (1064 nm) Gaussian beam with 1/e 2 diameter ≈ 50 µm that provides the initial confinement in the horizontal plane. Forced dipole evaporation occurs by lowering the intensity of both the blue-detuned vertical confinement beam and the red-detuned Gaussian beam until the condensate reaches a condensate fraction > 95 %. We estimate the initial temperature to be of the order of 50 nK by extrapolation of the evaporation process [65]. The final vertical trapping frequency is 650(4) Hz. The atoms are then adiabatically transferred to the initial trap for the experiment.
To create the initial ring (or target) trap, we use a direct intensity masking technique to create the blue-detuned (532 nm) trap in any shape. Details of this technique can be found in Refs. [66,67]. Briefly, this approach images the face of a digital micromirror device (DMD) that is illuminated by a blue-detuned Gaussian beam and imaged onto the atoms. The pattern written onto the DMD is then transferred onto the potential experienced by the atoms. Using this technique, we can form fully-dynamic potentials in the shape of rings (with radii between 10 µm and 45 µm) and target shaped traps (for measuring the persistent current state of the ring). The 1/e 2 radius of the Gaussian beam that illuminates the DMD is 130(10) µm in the plane of the atoms.
Nominally, the pattern written on the DMD is given by where Θ is the Heaviside step function, R T is the ring's width, and ρ is the radial coordinate. For rings thinner than R T < 10 µm, we apply corrections by changing R T with angle θ to make the measured n 1D (θ) density of the condensate more uniform.
To expand the ring, we apply a time-dependent potential using our DMD. To minimize spurious effects related to jerk, we used a smoothly varying function of the form where erf is the error function and β is a parameter that minimizes the jerk at t = 0 and t = t exp . For the data reported in this paper, β = 0.175, which implies that at t = 0 and t exp the radius suddenly jumps by ≈ 3 × 10 −5 (R f − R i ). The DMD is pre-programmed with individual frames with ring radii calculated using Eq. A2. We use approximately 30 frames spaced ≈ 300 µs apart to encode the expansion of the ring. Given our typical chemical potentials of ≈ 1 kHz, this update rate is faster than all other timescales in the system. Moreover, we checked that our results are independent of the number of frames used. During the expansion, we increase the intensity of the trapping light to maintain constant intensity locally near the ring (compensating for the Gaussian profile of the beam illuminating the DMD). We tune the increase in the trapping light to keep the frequency of the first radial Bogoliubov mode constant with radius.
To imprint a phonon of mode number m, we instantaneously change this pattern to Here λ = 0.6 is a parameter that describes the size of the perturbation relative to overall potential depth. One cannot generate the necessary values between 0 and 1 to produce the potential described by Eq. A3 with a binary DMD device. To get the necessary grayscale to create the potential, the DMD is demagnified in order to make its pixel size (≈ 0.5 µm in the plane of the atoms) be much smaller than the aberrated point spread function (≈ 4 µm 1/e 2 full-width) of our imaging system. We then use halftoning to create the necessary grayscale effect. Ref. [67] contains more details about this imprinting process.
To measure the normalized phonon amplitude after imprinting, we first measure the 2D density in situ with (n 2D (ρ, θ)) and without (n 2D,0 (ρ, θ)) the phonon imprinted. We then integrate over the radial dimension to obtain the 1D density around the ring, e.g., n 1D (θ) = n 2D (ρ, θ)dρ. To obtain the normalized 1D density, we compute n 1D (θ)/n 1D,0 (θ). The data are fit to a m sin(m(θ + θ 0 ) at each time to extract the normalized amplitude of the phonon a m (t). The offset angle θ 0 is set by the imprinting process. Finally, we turn a m (t) into real phonon amplitude δn(t) by multiplying by the total number of atoms and dividing by the estimated Thomas-Fermi volume of the condensate V T F . Here, we have made two implicit assumptions. First, we have assumed that the phonon's amplitude is independent of ρ and z, which is valid when the thickness of the annulus is small compared to its radius. (See Appendix C for details.) Second, we have assumed that the Thomas-Fermi volume scales in the experiment according to how it would in a potential that is quartic in ρ and harmonic in z. With the same assumptions on the potential, the predicted frequency shift scales as a 9/7 , which agrees rather well with the experiment (Fig. 2c). We also note that an incorrect estimate of the original Thomas-Fermi volume would lead to a common scaling of the phonon amplitude at all later times (before and after expansion), which would not lead to any change in either the fitted frequency shift or Hubble friction.
Calibration of the aberrations in our imaging system is necessary in order to accurately measure the correlation function S (k). Conveniently, PTAI allows us to accurately calibrate our imaging system's sensitivity to density structures with wavevector k. When the transfer fraction f is low ( f 1), quantum shot noise is added and dwarfs the thermal and quantum fluctuations inherent to the condensate. This additional noise is white over all k, thus allowing for accurate calibration. To calibrate, we measure n 1D (θ) as described above and then construct S (k) as described in the main text. To compensate for the our imaging system's degraded performance at larger k θ , we minimize the functional ( S (k) C(k) − 1 f ) 2 using the tunable parameters k 1 , p 1 , k 2 and k 3 contained within the correction function: The experimentally determined parameters are k 1 = 0.34(2) µm −1 , p 1 = 3.4(2), k 2 = 1.50(4) µm −1 and p 2 = 15(6).
To measure the persistent current state, we form a trap with a ring and a concentric, central disc (i.e., a target symbol) and use the interference between the two in time-of-flight to determine the winding number [68,69]. To produce acceptable interference fringes for readout, the disc must also be expanded. This is done adiabatically over 25 ms with 40 frames. Expansion of the ring produces a host of excitations, including phonons, vortices in the annulus, and persistent currents. To accurately measure the persistent current with the least amount of interference from other excitations, we let the ring equilibrate for about 5 s. During this period, the intensity of light is ramped to ≈ 60 % of its value at the end of expansion to force evaporation of high energy excitations.

Appendix B: Stochastic-Projected Gross-Pitaevskii Calculations
To explore the behavior of our system numerically, we conducted simulations of the stochastic projected Gross-Pitaevskii equation [42,43]. This numerical framework extends the ordinary Gross-Pitaevskii equation to non-zero temperature, adding on fluctuations to the BEC ground state. While described in detail in the aforementioned references, we will briefly describe the technique here. In this formalism, the wavefunction of the BEC with fluctuations evolves in a "coherent" region -defined as the region of Hilbert space spanned by the state vectors that impact the dynamics of the BEC coherently. The BEC wavefunction in this C-region evolves as where (S ) denotes Stratonovich integration and Here, L = H sp + U 0 |ψ| 2 is the driver of Hamiltonian evolution and H sp = p 2 /2M + V is the single particle Hamiltonian. The equation for dψ G represents growth of population in the C-region from particles colliding in the incoherent (I) region.
Here, G(r) is a coefficient that sets the strength of both terms in Eq. B3, where the first term is the damping term and the second is the growth term where dW G describes a random noise seeded according to dW * G (r , t )dW G (r, t) = 2G(r)δ(r −r)dt. For this work, we neglect terms where there is an exchange of energy and momentum between the C and I region without exchange of particles [42]. Finally, the projector operator P continually projects the wavefunction into the C region.
From an implementation perspective, this involves taking a Gross-Pitaevksii equation solver and adding a noise term, and appropriately calculating the damping factor G(r), which is assumed to be constant. Our calculations are done in a Cartesian coordinate system. We apply the projection operator in momentum space, with a cutoff k c ≈ π/δx, where δx is the spacing between points in the grid.
To accurately capture the potential, we simulate the imaging process that is used to make the potential. We reproduce the image that is patterned on the DMD and simulate imaging using Fourier imaging techniques. The aperture function of the imaging system that relays the image from the DMD to the atoms is crucial in order to accurately replicate the potential at the atoms. In the experiment, the same imaging system that is used for making potentials is also used for imaging of the atoms. By measuring density-density correlations in a simple-connected thermal gas with noise dominated by quantum shot noise (by using f 1), we can extract the even (symmetric under parity reversal) aberrations [70]. To extract the odd aberrations, we use a less precise means. A second DMD in the Fourier plane of the imaging system can be used to measure the geometric spot diagram, yielding another, independent means of obtaining the aperture function. The two methods are in agreement. We use the even aberrations from the correlations and the odd aberrations for the spot diagram technique to construct the aperture function. Finally, we neglect variations in the intensity of the trapping light caused by unwanted scattering along the optical path (i.e., speckle). Because the atoms seek the darkest part of the imaged potential (the trap is blue detuned), the effect is minimized. Speckle from the surface of the DMD is eliminated by imaging.
We combine the aperture function with the Gaussian beam. We assume the beam is perfectly Gaussian and is centered on the DMD. Because the beam portion of the potential tends toward zero as r → ∞, we establish a low energy potential floor at large radius. This cutoff is determined by the minimum value of the imaged and aberrated potential between R + 3 2 R t < ρ < R + 1.1 × 3 2 R t . This prevents spurious effects like the appearance of additional BEC components out at large radius.
The resulting potential is complicated and not easily expressible in an analytic form. However, when azimuthally averaged, the potential has the form V = 1 2 Mω 2 r (ρ − R(t)) 2 + λ(ρ − R(t)) 4 , with ω r ≈ 2π × 100 Hz and λ/h ≈ 0.8 Hz µm −4 . Because most (≈ 90 %) of the confinement comes from the quartic term, it is generally acceptable to neglect the quadratic term for the purposes of calculating static properties like the initial and final µ and c.
Given that some atoms are lost during the expansion, we must also include absorbing boundary conditions in the simulation. We do this by including a potential where R c is a radial cutoff at which the potential turns on, V a is the amplitude of the potential, and w a is a parameter that controls the width. The function is a non-analytic, continuously differentiable function that minimizes the reflections from the absorbing boundary. We chose w a ≈ 25 µm and V a /h ≈ 1 kHz. This generally results in the best absorption and the least reflection. With all of these components, the simulations proceed as follows. We first find the equilibrium state by evolving the SPGPE (without V a ) for approximately 50 ms to 100 ms using a growth and decay term that are 100 times that of the value specified by the temperature (this allows for faster equilibration times). Second, we expand the ring according to that the same profile as seen in the experiment. Approximately halfway through the experiment, we turn on V a to ensure that the decay of the atom number is appropriately captured. After evolving for a total of approximately 35 ms (20 ms additional after the end of the expansion), we turn off the stochastic growth term in the SPGPE and turn on significant damping to determine whether or not a winding number is present in the condensate. We do this approximately 25 independent times to gather statistics.
We then use the same data analysis tools used on the experimental data to extract the winding number distributions, structure factor as a function of time, and radius of the ring as a function of time. The structure factor, as was done in the experiment, is measured relative to the mean density around the ring. As a result, the structure factor is determined solely by the differences in density between a given simulation and the mean of all the simulations. This density extends out to the contour in theρ-z plane where the numerator vanishes. The chemical potential µ drops as the ring expands, so that the total number of atoms remains constant.
In the experiment, we first apply a perturbation to a stationary condensate to excite an eigenmode of the wave equation. An eigenmode analysis based on the methods of Ref. [71] will be detailed in a forthcoming paper; the essential details are presented here. Assuming azimuthal symmetry, the eigenmodes for a thin ring have the form φ 1 = χ klm η klm (ρ, z; R)e i(ωt−mθ) , where η klm (ρ, z; R) is a function that describes the radial (k) and vertical excitations (l) of the Bogoliubov mode when the ring has radius R, and χ klm is its amplitude. We denote the corresponding eigenfrequencies as ω klm .
While the system may begin with only a k = l = 0 eigenmode excited, the expansion of the ring can produce transitions into other modes. The solution at all times takes the general form φ 1 (t,ρ, θ, z) = klm χ klm (t)η klm (ρ, z)e −imθ , with all χ klm (t = 0) = 0 except for k, l = 0 and our excited mode of interest m. Azimuthal symmetry precludes coupling between modes with different values of m. Furthermore, in the thin ring limit, the coupling between different k modes tends towards zero. We therefore focus here exclusively on modes that are excited only in the azimuthal direction. (The radial excitation which occurs when the ring expansion stops and is not relevant to the redshift is discussed in the main text.) When m ω 100 /(c θ /R) and for a thin ring, η 00m (ρ, z; R) is constant. (Henceforth, we will drop the k and l subscripts when they are both equal to zero.) In this limit, the equation for modes with k, l = 0 involves just t and θ derivatives, We can thus reduce the wave equation for these azimuthal phonon modes to a 1+1 dimensional wave equation, with an effective sound speed c θ . As in Ref. [71], c 2 θ is given by an average over the cross section of the ring. For a thin ring this takes the form where the integral is over the cross section of the Thomas-Fermi wavefunction of area A. For V = 1 2 Mω 2 z z 2 + λρ 4 this yields where c 2 = µ/M is the peak local sound speed. By normalizing the Thomas-Fermi solution to the number of atoms N one finds that µ ∝ R −4/7 , and therefore c θ ∝ c ∝ R −2/7 . The wave equation satisfied by our modes of interest, i.e., φ 1 = χ m (t)e imθ , is determined by the effective inverse metric density obtained from Eq. C8 by dropping theρ and z components and replacing c by c θ . In the thin ring limit this gives The resulting mode equation is where ω m := mc θ /R. This is the equation of a damped harmonic oscillator, with time-dependent frequency and damping rate. We note that this particular equation does not result from the wave equation for any 1+1 dimensional metric, since there exists no metric g 2ab for which f ab 2 = √ −gg ab . The reason is that the determinant of Eq. C12 is −c 2 θ , whereas the determinant of √ −gg ab is equal to −1 for any two-dimensional metric.
As the ring expands, the azimuthal wavenumber m is conserved, so the physical wavelength redshifts as R −1 , in analogy with the cosmological redshift. Unlike in cosmology, the sound speed is also changing, so the frequency ω m redshifts as R −9/7 . In the cosmological setting, the damping term in Eq. C13 is called the "Hubble friction" term, and would be multiplied by 3 in three spatial dimensions. The Hubble damping is not actually dissipative; in fact, Eq. C13 can be obtained from the Lagrangian L = 1 2 Rχ 2 m − 1 2 (m 2 c 2 θ /R)χ 2 m , which has the adiabatic invariant Rω m χ 2 m . To obtain Eq. 1 in the text, we add the phenomenological damping γ m observed in the experiment.
In the experiment, we measure the density variation n 1 , not the phonon velocity potential φ 1 . The relation between these quantities is given by Eq. C4. Since ∇φ 1 is azimuthal and ∇φ 0 is radial, ∇φ 0 · ∇φ 1 = 0, so we have Hence, in the experiment, we measure the time derivative of the phonon amplitude. We can verify that a phonon excitation does indeed obey Eq.C13 in a thin ring, by simulating a BEC in this regime. Figure 6 shows such a 2D simulation of BEC in a radially quartic potential, expanding from 10 to 40 µm in ≈ 15 ms with 2 × 10 5 atoms. There is no damping in this simulation; therefore, the γ m (t) in Eq. 1 is identically zero. We choose the strength of the potential to make the initial Thomas-Fermi width be 2 µm. As can be seen from the figure, Eq. C13 accurately reproduces the behavior of the redshifted phonon, but only when the Hubble friction term is included. Unlike the experiment, the adiabatic limit is satisfied (∂ω m /∂t ω 2 m ) and the final amplitude is accurately predicted using the adiabatic invariant Rω m χ 2 m . While this simulation shows the thin ring limit, we generally find that as we relax this constraint and increase the width of the annulus, the best-fit Hubble friction becomes less than unity, as might be expected from the experimental result.
For the experiment, we attempt to tease out the Hubble friction by fitting it along with the other parameters of Eq. 1. These parameters are the initial amplitudes, frequencies, and phases for each of the four expansions, the quality factor for the two m modes, Q m = ω m /2γ m , the scaling of the frequency with radius [expected to be 9/7 ≈ 1.2875, the best fit value is 1.19 (2)]. This fit therefore contains 15 parameters and 160 degrees of freedom. Fig. 7 shows the results of the reducedχ 2 fit; it shows the value of χ 2 vs. both Q m=1 and Q m=2 in the vicinity of their best fit values for three values of γ H , including the best fit value. There are several interesting features. First, the reduced-χ 2 > 1, most likely because we do not have a good estimate of the statistical uncertainties (each point represents only four realizations of the experiment) and our model does not properly account for all the relevant effects (for example, the azimuthal asymmetry and non-zero annular thickness may play a non-negligible role in determining the phonon dynamics). Second, γ H = 1 produces a better fit than γ H = 0, but both are improved slightly by taking γ H = 0.55. Third, the smallness of the change in the minimum of χ 2 with γ H indicates our uncertainty it γ H . (Part of this insensitivity comes our choice, at the time of the experiment, to have n 1 ∝φ ≈ 0 during the fastest part of the expansion, thereby inadvertently minimizing the effect of the Hubble friction [72].) Taken together, the evidence is consistent with γ H = 1 but is not conclusive.