PHYSICAL REVIEW X 8, 021021 (2018)
Featured in Physics

A Rapidly Expanding Bose-Einstein Condensate: An Expanding Universe in the Lab
S. Eckel,1 A. Kumar,1 T. Jacobson,2 I. B. Spielman,1 and G. K. Campbell1,*
1

Joint Quantum Institute, National Institute of Standards and Technology and University of Maryland,
Gaithersburg, Maryland 20899, USA
2
Department of Physics, University of Maryland, College Park, Maryland 20742, USA(Received 16 October 2017; revised manuscript received 3 February 2018; published 19 April 2018)

We study the dynamics of a supersonically expanding, ring-shaped Bose-Einstein condensate bothexperimentally 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.
DOI: 10.1103/PhysRevX.8.021021

Subject Areas: Atomic and Molecular Physics,
Cosmology, Quantum Physics

I. INTRODUCTION
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 theredshifting of photons, which provided the first evidence foran 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 transferredfrom a homogeneous radial mode into azimuthal modes by anonlinear, turbulent cascade, reminiscent of that seen in somemodels of preheating at the end of cosmological inflation[3–9]. Experiments such as these can thus emulate bothlinear 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, theequation for long-wavelength phonons (in the hydrodynamic limit [11]) takes on a covariant form with a curved
spacetime metric [12–14]. Previous studies with ultracold
atoms illuminated different aspects of this phonon metric.
*

Corresponding author.
gcampbe1@umd.edu
Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article’s title, journal citation,
and DOI.

2160-3308=18=8(2)=021021(13)For example, an interface between regions of subsonic andsupersonic fluid flow forms a “sonic event horizon” thatexhibits effects such as Hawking radiation [12,14–22].
By changing the interaction strength or density, one cansimulate cosmological phenomena such as pair production[23–29], Sakharov oscillations [30], or the dynamical
Casimir effect [31]—the latter two were recently observed
experimentally [32,33]. Beyond cold atoms, experimental
studies have realized analog event horizons in other settings,
for example, in optical systems [34–36] and in classicalfluids [37–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 ananalog of the Einstein equation (see Appendix A). Figure 1
shows our BEC during a texp ¼ 23.1 ms inflation. The radial
velocity of the trapping potential (defined as the rate ofchange of the mean radius, R) is directly controlled and canbe made comparable to the speed of sound. For the expansion
shown in Fig. 1, the maximum velocity is vp ¼ 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 thestochastic-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.
II. PHONON REDSHIFTTo study the redshifting of phonons, we first imprint astanding-wave phonon excitation on the background
BEC. During expansion, these effectively one-dimensional

021021-1

Published by the American Physical Society

S. ECKEL et al.

PHYS. REV. X 8, 021021 (2018)

FIG. 1. Measured (top) and simulated (bottom) supersonic expansion of the ring with scale factor a ¼ Rf =Ri ¼ 4.1ð3Þ, whereRf ¼ 46.4ð1.4Þ μm [Ri ¼ 11.3ð4Þ μm] is the final (initial) radius [44]. An initial distance d transforms into a larger distance d0 .
The time elapsed in the figure is approximately 15 ms.

azimuthal phonons are redshifted; i.e., their wavelength
grows as shown in Fig. 2(a) for both experiment and theory.
These images show the oscillation of a standing-wave
phonon, constructed by perturbing the condensate with apotential 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 withconservation of the comoving wave vector in cosmology. The
phonon wavelength is therefore stretched by a factor
a ¼ Rf =Ri , the ratio of the geometrical radii of the expanding
ring. This is related to the usual redshift parameter z
through a ¼ z þ 1.
Figure 2(b) shows the measured phonon amplitude δn vstime for various a and m, and it clearly shows a shift in thefrequency. [In this paper, we measure frequency and time inthe laboratory frame, as opposed to using the comoving(a)proper time as defined by the effective metric, Eq. (C7).] To
measure the frequency shift, we fit the oscillation before andafter expansion to extract ωm;i =ωm;f , shown in Fig. 2(c). At
any given time, the phonon oscillation frequency isωðtÞ ¼ cθ ðtÞm=RðtÞ, where cθ ðtÞ is the azimuthal speedof 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 ¼ a9=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 ,

(b)

(c)

FIG. 2. Redshift of long-wave excitations. (a) Atomic density difference δn at various times for both experiment and simulation for amode number m ¼ 2 and scale factor a ¼ 4.1. The density scale of images after expansion is multiplied by 10. (b) Phonon amplitude vstime for various a and m. The grey bands indicate the time during which the BEC is inflated; their intensity denotes the expansion
velocity relative to that expansion’s maximum. Vertical dashed lines in the m ¼ 2, a ¼ 4.1 panel indicate the times shown in (a).
(c) Ratio of initial to final frequency vs scale factor a. Red circles indicate m ¼ 1 modes; cyan squares, m ¼ 2. The solid black curve isthe a9=7 expectation, and the colored curves (with mode numbers matching points) are the result of full Bogoliubov calculation.

021021-2

A RAPIDLY EXPANDING BOSE-EINSTEIN CONDENSATE: …


∂ 2χm
R_ ∂χ mþ2γðtÞþ
þ ½ωðtÞ2 χ m ¼ 0;
m
R ∂t
∂t2

PHYS. REV. X 8, 021021 (2018)

ð1Þ

where δn ¼ ðℏ=U0 Þ∂χ m =∂t, U 0 ¼ 4πℏ2 as =M, as is thes-wave scattering length, and M is the mass of an atom.
(See Appendix C for the derivation.) There are twocontributions to the damping of the amplitude. The first
damping term, γ m , is phenomenological, but independently_
measured [45]. The second, R=R,
is analogous to the“Hubble friction” in cosmology, which damps fields with_frequencies ω ≲ a=a.
In the present case, the Hubble
friction has the largest impact for supersonic expansion,__i.e., when ω ≲ R=R
or mcθ ≲ R.
For our expansions, we expect that the Hubble friction
will play a role, particularly for the a ¼ 4.1, m ¼ 1_
expansion where R=R≳ 1.5ω. (At maximum velocity,_R=R ≈ 3γ m for m ¼ 2 and R=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 inthe observed amplitude before and after expansion resultsfrom the change in ω. To search for the Hubble friction
term, we fit all the data simultaneously to Eq. (1), taking__
R=R→ γ H R=R,
where γ H is a tunable parameter. While thebest-fit value γ H ¼ 0.55ð21Þ indicates the presence of
Hubble friction, the deviation from unity suggests that
other effects, like azimuthal asymmetry and nonzero
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 similarfunctional form, Eq. (1) is an accurate description of thephonon evolution (see Appendix C).
III. RADIAL DYNAMICS
The preceding 1D discussion [based on Eq. (1)] rested on
the assumption that the background BEC contained notransverse 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 andforces applied are due to the expansion of the confining
potential. These oscillations decay rapidly, typically withina few oscillation periods for all scale factors and expansion
velocities studied. If the trap were perfectly harmonic, thiscenter-of-mass oscillation should be long-lived. In reality,
our trap is more flat bottomed and anharmonic, and it is notaxially symmetric.

FIG. 3. Mean radius of the ring vs time for select expansions.
Black solid curves show the radius of the potential, red circles
show experimental data, and orange curves show simulation
results. The vertical lines for a ¼ 1.9ð1Þ correspond to times
shown in Fig. 4.

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 Figs. 4(a) and 4(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 approximately π in the centerof the annulus. This standing wave has nodes in the atomic
density with corresponding π phase jumps, effectively
imprinting a dark soliton onto the BEC [47,48]. This process
is analogous to the creation of solitons upon Bragg reflection
in an optical lattice [49] or reflection of a condensate off of atunnel barrier [50]. (Unfortunately, because of imaging
limitations, we are unable to resolve solitons or othersimilarly sized structures in the experiment.)
The number of solitons N s created from the decayof the radial mode can be estimated by comparing theenergy per particle contained in the radial excitation to theenergy per particle of a soliton (ϵs ≈ 4ℏc=3RT , where RT
is the annular width of the ring). The amplitude of theradial excitation χ r , while calculable analytically, is acomplicated function that depends exponentially on the
adiabaticity of the expansion relative to the frequency of
the radial mode ωr . (Assuming a boxlike potential in theradial direction implies ωr ≈ πc=RT .) The adiabatic condition demands that, in our system, R_ must be nearlysupersonic to produce solitons, i.e., vp ≳ 0.8c [51].

021021-3

S. ECKEL et al.

PHYS. REV. X 8, 021021 (2018)

(a)

(b)

(c)

FIG. 4. Dynamics in the radial direction and generation of azimuthal excitations for a scale factor of a ¼ 1.9ð1Þ. (a,b) A single
realization of the simulated BEC wave function’s phase (a) and magnitude jψj2 (b) vs time. (Left inset) Cross section along ŷ of thephase and density at t ¼ 11.5 ms. (Right inset) Zoomed section of the phase profile at t ¼ 13 ms, with the location of vortices
highlighted. (c) Azimuthal static structure factor. The color indicates the time and matches the vertical lines in Fig. 3 and the rectangles
in (a) and (b). The horizontal lines indicate the imaging detection threshold (see text).

IV. 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 [52–54]. As shown by our numerics in
Fig. 4(a), the undulation is underway by 12.5 ms, and thefragmentation into vortices is mostly complete by 14 ms.
Theoretical estimates for a single soliton in a harmonicallyconfined BEC suggest that the snake instability
will result
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ﬃ
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 [55]. For the present case,this corresponds to N vd;1 ≈ 50 vortex pairs. At t ¼ 13 ms in
Fig. 4(b), the single soliton has decayed into approximately
6 pairs over an angle of approximately 45° near the top of
the ring. This corresponds to roughly 48 vortex pairsaround the full ring. These vortex pairs then form a highly
turbulent state.
We experimentally observe the fingerprints of thisprocess through the structure factor Sðkθ Þ, a measure of
the spatially structured density fluctuations (i.e., azimuthal

phonons) with wave vector kθ ¼ m=R. For both experimentand theory, we extracted Sðkθ Þ by first evaluating the onedimensional density n1D ðθÞ around the ring to obtain thedensity fluctuations δn1D ðθÞ ¼ n1D ðθÞ − hn1D ðθÞi, whereh…i denotes the average over many realizations. The
structure factor is
2 
 Z


−ikθ RθSðkθ Þ ¼  δn1D ðθÞe
ð2ÞRdθ :
Theoretical structure factors are shown in the top row ofFig. 4(c); experimental structure factors are shown in thebottom row of Fig. 4(c). The colors in Fig. 4(c) 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 [56]. For these reasons, we first correct for imaging
aberrations (see Appendix A) and identify the detection
threshold (shown by the horizontal lines). We use thenumerical simulations (which include the same aberrations)to verify the correspondence between the corrected value of021021-4

A RAPIDLY EXPANDING BOSE-EINSTEIN CONDENSATE: …

PHYS. REV. X 8, 021021 (2018)(a)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 [57]. After expansion andduring 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 ourdetection threshold, appears in the simulations near k ≈1 μm−1 (not shown). This roughly corresponds to the wavenumber of the snake instability, k ≈ 2π=8ξ ≈ 1.3 μm−1 . As
the turbulent state develops, this peak grows and shifts tolower k, becoming detectable at 18.5 ms and becoming
larger at 22.5 ms. The shift to lower kθ is expected becauseof the inverse cascade that occurs in two-dimensional
turbulence [58].

(b)V. 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 freevortices slides into the center of the ring and one leaves theoutside of the ring, then the overall phase of the ring slips
by 2π and the winding number l, quantifying the persistentcurrent state of the ring, changes by 1 [59]. Indeed, we
observe stochastic persistent currents in the ring after
expansion in both the experiment and simulation.
Figure 5(a) shows the resulting distributions of winding
numbers for various speeds of expansion for a ¼ 1.4ð1Þ
[Rf ¼ 35ð2Þ μm and Ri ¼ 25ð1Þ μm].
Evidence for this process can be found by studying thewidth of the winding number distributions for expansions
with different a and texp . The number of vortex dipoles
produced from N s solitons would be N vd ≈ N s ð2πRf =8ξÞ.
The measured distribution widths collapse reasonably wellwhen plotted versus N 1=4
vd , as shown in Fig. 5(b). The 1=4
may result from some combination of the stochastic natureof dipole dissociation and recombination, the interactiondriven 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 speedof 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 establishedduring the expansion would be given byZ τ
Rhor ¼ 2
cðtÞdt ≲ 2c0 texp ;ð3Þ
0where c0 is the initial speed of sound [60]. This leads to
N R ≳ 2πRf =Rhor ≈ πRf =c0 texp disconnected regions. If

FIG. 5. (a) Measured winding number distributions for a ¼
1.4ð1Þ (histograms) with Gaussian fits (black curves) for theexperiment (top) and simulations (bottom). The gray vertical line
and bar show, for comparison, the speed of sound and its
uncertainty. The red Gaussians show the expected distribution
from the horizon model (see text). (b) Width of the winding
number distribution σ vs the number of vortex dipoles predictedby the soliton model, N vd . The red circles denote a ¼ 1.31ð9Þ
[Ri ¼ 16.7ð8Þ μm, Rf ¼ 22ð1Þ μm]; yellow squares, a ¼ 1.9ð1Þ
[Ri ¼ 16.7ð8Þ μm, Rf ¼ 31ð2Þ μm]; green triangles, a ¼ 2.3ð2Þ
[Ri ¼ 16.7ð8Þ μm, Rf ¼ 39ð2Þ μm]; and blue inverted triangles,
a ¼ 1.4ð1Þ [Ri ¼ 25ð1Þ μm, Rf ¼ 35ð2Þ μm]. The purple diamonds also show a ¼ 1.4ð1Þ but with a ring that is twice as wide.
The solid line shows σ ¼ N 1=4
vd .

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 apersistent current [61,62]. The probability for a given
persistent current is then given by the geodesic rule
[62–64]. The red Gaussians in Fig. 5(a) show the expecteddistributions resulting from this horizon model, which
disagree with the experiment. Moreover, simple estimates
for the phase fluctuations present in our condensate are afactor of 25 too low to sufficiently randomize the phaseduring expansion. Future studies using condensates oflower density could see this effect, as the phase fluctuations
will be larger.
VI. DISCUSSION AND OUTLOOKIn this work, we explored the physics of a rapidly
expanding Bose-Einstein condensate. We observed theredshifting of phonons during this rapid expansion, which021021-5

S. ECKEL et al.

PHYS. REV. X 8, 021021 (2018)

has clear analogs in cosmological physics. After expansion
stops, the condensate reheats through the creation andsubsequent 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 firstmight be thought to arise from 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=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 isreversible, such a contraction should cause amplification of
the phonon mode amplitude.
The process of expansion, which presumably cools theazimuthal degrees of freedom of the condensate, followedby the increase in azimuthal excitations [Figs. 4(c)and 4(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 thehomogeneous mode of the quantum field that drove
inflation, the inflaton, decays into inhomogeneous excitations. It is not known how this occurs. In the simplest
model, the inflaton oscillates 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 radialmode through this process is expected to be much slower(approximately 1 s−1, using a calculation similar to thatfound in Ref. [46]) compared to the observed decay of
the radial mode through soliton and vortex excitations(approximately 100 s−1 ). Future studies using a ring withstronger radial confinement should suppress the nonlinear
excitations and enhance the direct coupling. Other models
are nonperturbative and include self-interactions in the
inflaton field that can lead to turbulent cascading [4–9],
much like the turbulence we observe here.
Perhaps surprisingly, the long-wavelength azimuthal
phonon mode is redshifted in a 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 wereredshifted 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 long-wavelength amplified
vacuum fluctuations. This process gave rise to the largescale structure we observe today in the Universe.
In addition to the possibilities described above, we
anticipate that with new developments, other interestingcosmological phenomena might be realized with expanding
condensates. First, with improved imaging that captures theinitial (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–29]. Second, a ring with stronger radialconfinement 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 interestinglaboratory test bed for cosmological physics.
ACKNOWLEDGMENTS
The authors thank J. Ho for initial discussions and acareful 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
Physics Frontier Center at the Joint Quantum Institute. T. J.
was supported in part by NSF Grants No. PHY-1407744
and No. PHY-1708139. I. B. S. was partially supportedby the AFOSR’s Quantum Matter Multidisciplinary
University Research Initiative (MURI) and NIST.
APPENDIX A: EXPERIMENTAL DETAILSOur experimental setup consists of a BEC of 23Na atoms
in an optical dipole trap (ODT). Our BECs are createdusing standard laser cooling techniques, followed byevaporation first in magnetic and then in optical dipole
traps. In this experiment, we work with BECs with between
1 × 105 and 4 × 105 atoms. For measurements, we use
partial transfer absorption imaging (PTAI) [56].
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), TEM01 beam, tightly focused to create twoparallel sheets of light with a dark region in between. The
dimple trap is a red-detuned (1064-nm) Gaussian beam with
1=e2 diameter of approximately 50 μm, which provides theinitial confinement in the horizontal plane. Forced dipole
evaporation occurs by lowering the intensity of both theblue-detuned vertical confinement beam and the reddetuned Gaussian beam until the condensate reaches acondensate fraction greater than 95%. We estimate theinitial temperature to be of the order of 50 nK by extrapolation of the evaporation process [65]. The final verticaltrapping frequency is 650(4) Hz. The atoms are thenadiabatically 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 befound in Refs. [66,67]. Briefly, this approach images theface of a digital micromirror device (DMD) that isilluminated by a blue-detuned Gaussian beam onto theatoms. The pattern written onto the DMD is then021021-6

A RAPIDLY EXPANDING BOSE-EINSTEIN CONDENSATE: …
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 persistentcurrent state of the ring). The 1=e2 radius of theGaussian 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
V DMD ðrÞ ¼ Θð½RðtÞ − RT =2 − ρÞ þ Θðρ − ½RðtÞ þ RT =2Þ;ðA1Þwhere Θ is the Heaviside step function, RT is the ring’s
width, and ρ is the radial coordinate. For rings thinner thanRT < 10 μm, we apply corrections by changing RT with
angle θ to make the measured n1D ðθÞ density of thecondensate more uniform.
To expand the ring, we apply a time-dependent potential
using our DMD. To minimize spurious effects related to thejerk, we used a smoothly varying function of the form
8
>
>
> Ri
<

n 
o
t
− 12
RðtÞ ¼ Ri þ 12 ðRf − Ri Þ 1 þ erf β1 texp
>
>
>
:
Rf

t≤0
0 < t ≤ texpt > texp ;ðA2Þ

where erf is the error function and β is a parameter thatminimizes the jerk at t ¼ 0 and t ¼ texp . For the datareported in this paper, β ¼ 0.175, which implies that at
t ¼ 0 and texp the radius suddenly jumps by approximately 3 × 10−5 ðRf − Ri Þ. The DMD is preprogrammedwith 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 typicalchemical potentials of approximately 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 frequencyof the first radial Bogoliubov mode constant with radius.
To imprint a phonon of mode number m, we instantaneously change this pattern to




Rt
Rt
V DMD ðr; θÞ ¼ Θ RðtÞ −− ρ þ Θ ρ − RðtÞ þ
2
2




λ
R
R
þ Θ ρ − RðtÞ − t Θ RðtÞ þ t − ρ2
2
2
× ½1 þ sinðmθÞ:

ðA3Þ

PHYS. REV. X 8, 021021 (2018)

Here, λ ¼ 0.6 is a parameter that describes the size of theperturbation 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 gray scale to create thepotential, the DMD is demagnified in order to make its
pixel size (approximately 0.5 μm in the plane of the atoms)
much smaller than the aberrated point-spread function(approximately 4 μm 1=e2 full width) of our imaging
system. We then use halftoning to create the necessary
gray scale effect. Reference [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[n2D ðρ; θÞ] and without [n2D;0 ðρ; θÞ] the phonon imprinted.
We then integrate over the radial dimensionR to obtain the
1D density around the ring, e.g., n1D ðθÞ ¼ n2D ðρ; θÞdρ.
To obtain the normalized 1D density, we computen1D ðθÞ=n1D;0 ðθÞ. The data are fit to am sinðmðθ þ θ0 Þ ateach time to extract the normalized amplitude of the phononam ðtÞ. The offset angle θ0 is set by the imprinting process.
Finally, we turn am ðtÞ into real phonon amplitude δnðtÞ bymultiplying by the total number of atoms and dividing by theestimated Thomas-Fermi volume of the condensate, V TF .
Here, we made two implicit assumptions. First, we 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
assumed that the Thomas-Fermi volume scales in theexperiment according to how it would in a potential thatis quartic in ρ and harmonic in z. With the same assumptions
on the potential, the predicted frequency shift scales as a9=7 ,which agrees rather well with the experiment [Fig. 2(c)]. We
also note that an incorrect estimate of the original ThomasFermi volume would lead to a common scaling of thephonon 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 isnecessary 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 wave vector k. When the transfer fraction f islow (f ≪ 1), quantum shot noise is added, and it dwarfs thethermal and quantum fluctuations inherent to the condensate. This additional noise is white over all k, thus allowing
for accurate calibration. To calibrate, we measure n1D ðθÞ as
described above and then construct SðkÞ as described in themain text. To compensate for our imaging system’s
degraded performance at larger kθ , we minimize the functional f½SðkÞ=CðkÞ − ð1=fÞg2 using the tunable parameters k1 , p1 , k2 , and k3 contained within the correction
function:

021021-7

S. ECKEL et al.
1
1
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ :CðkÞ ¼ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
p1
1 þ jk=k1 j
1 þ jk=k2 jp2

PHYS. REV. X 8, 021021 (2018)
ðA4Þ

The experimentally determined parameters are k1 ¼
0.34ð2Þ μm−1 , p1 ¼ 3.4ð2Þ, k2 ¼ 1.50ð4Þ μm−1 , andp2 ¼ 15ð6Þ.
To measure the persistent current state, we form a trapwith a ring and a concentric central disc (i.e., a target
symbol) and use the interference between the two in time offlight to determine the winding number [68,69]. To produce
acceptable interference fringes for readout, the disc mustalso be expanded. This is done adiabatically over 25 mswith 40 frames. Expansion of the ring produces a host of
excitations, including phonons, vortices in the annulus, andpersistent currents. To accurately measure the persistentcurrent with the least amount of interference from otherexcitations, we let the ring equilibrate for about 5 s. During
this period, the intensity of light is ramped to approximately
60% of its value at the end of expansion to force
evaporation of high-energy excitations.
APPENDIX B: STOCHASTIC-PROJECTEDGROSS-PITAEVSKII CALCULATIONS
To explore the behavior of our system numerically, we
conducted simulations of the stochastic projected GrossPitaevskii equation [42,43]. This numerical frameworkextends the ordinary Gross-Pitaevskii equation to nonzero
temperature, adding fluctuations to the BEC ground state.
While described in detail in the aforementioned references,
we briefly describe the technique here. In this formalism,
the wave function of the BEC with fluctuations evolves in a
“coherent” region—defined as the region of Hilbert spacespanned by the state vectors that impact the dynamics of the
BEC coherently. The BEC wave function in this C region
evolves asðSÞdψ ¼ dψ H þ dψ G þ ðSÞdψ M ;ðB1Þ

where (S) denotes Stratonovich integration andidψ H ¼ P − Lψdt ;ℏdψ G ¼ P

GðrÞðμ − LÞψdt þ dW G ðr; tÞ :
kB T

ðB2Þ
ðB3Þ

Here, L ¼ Hsp þ U 0 jψj2 is the driver of Hamiltonian
evolution, and Hsp ¼ p2 =2M þ V is the single-particle
Hamiltonian. The equation for dψ G represents growth of
population in the C region from particles colliding in theincoherent (I) region. Here, GðrÞ is a coefficient that sets
the strength of both terms in Eq. (B3), where the first termis the damping term and the second is the growth term, withdW G describing a random noise seeded according to

hdW G ðr0 ; t0 ÞdW G ðr; tÞi ¼ 2GðrÞδðr0 − rÞdt. For this work,
we neglect terms where there is an exchange of energy andmomentum between the C and I regions without exchangeof particles [42]. Finally, the projector operator P continually projects the wave function into the C region.
From an implementation perspective, this involves
taking a Gross-Pitaevskii equation solver and adding anoise term, and appropriately calculating the damping
factor GðrÞ, which is assumed to be constant. Our calculations are performed in a Cartesian coordinate system. We
apply the projection operator in momentum space, with acutoff kc ≈ π=δx, where δx is the spacing between points in
the grid.
To accurately capture the potential, we simulate theimaging process that is used to make the potential. We
reproduce the image that is patterned on the DMD andsimulate imaging using Fourier imaging techniques. The
aperture function of the imaging system that relays theimage from the DMD to the atoms is crucial in order toaccurately replicate the potential at the atoms. In theexperiment, the same imaging system that is used formaking potentials is also used for imaging of the atoms. By
measuring density-density correlations in a simplyconnected 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 aperturefunction. The two methods are in agreement. We use theeven aberrations from the correlations and the odd aberrations for the spot diagram technique to construct the aperturefunction. Finally, we neglect variations in the intensity ofthe trapping light caused by unwanted scattering along theoptical path (i.e., speckle). Because the atoms seek thedarkest part of the imaged potential (the trap is blue
detuned), the effect is minimized. Speckle from the surfaceof the DMD is eliminated by imaging.
We combine the aperture function with the Gaussianbeam. We assume the beam is perfectly Gaussian and iscentered 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 andaberrated potential in the range R þ 32 Rt < ρ < Rþ
1.1 × 32 Rt . This prevents spurious effects like the appearance of additional BEC components at a 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 ¼ 12 Mω2r (ρ−RðtÞ)2 þ
4˜λ(ρ−RðtÞ), with ωr ≈ 2π × 100 Hz and λ=h ≈ 0.8 Hzμm−4 .
Because most (approximately 90%) of the confinement
comes from the quartic term, it is generally acceptable to021021-8

A RAPIDLY EXPANDING BOSE-EINSTEIN CONDENSATE: …
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, wemust also include absorbing boundary conditions in thesimulation. We do this by including a potential
Va ¼

0

ρ < Rc−wa =ðρ−Rc Þ−iV a eρ > Rc ;ðB4Þ

where Rc is a radial cutoff at which the potential turns on,
V a is the amplitude of the potential, and wa is a parameterthat controls the width. The function is a nonanalytic,
continuously differentiable function that minimizes thereflections from the absorbing boundary. We choose wa ≈25 μm and V a =h ≈ 1 kHz. This choice generally results inthe best absorption and the least reflection.
With all of these components, the simulations proceed asfollows. We first find the equilibrium state by evolving the
SPGPE (without V a ) for approximately 50 ms to 100 ms,
using growth and decay terms that are 100 times that of thevalue specified by the temperature (this allows forfaster equilibration times). Second, we expand the ringaccording to 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 isappropriately 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 thecondensate. We do this approximately 25 independent
times to gather statistics.
We then use the same data analysis tools as on theexperimental data to extract the winding number distributions, the structure factor as a function of time, and theradius of the ring as a function of time. The structure factor,
as was done in the experiment, is measured relative to themean 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 thesimulations.
APPENDIX C: EVOLUTION
OF AZIMUTHAL PHONONSIn this appendix, we derive the wave equation satisfied
by the phonon field and explain the origin of the redshiftingand Hubble friction. We start by noting that the GrossPitaevskii equation, written in termsof thepﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ﬃ density n andphase ϕ defined through ψ ¼ nðr; tÞeiϕðr;tÞ , can be
expressed as an ideal fluid with an equation for continuity,−∂n ℏ¼ ∇ · ðn∇ϕÞ;
∂t m

ðC1Þ

−ℏ

PHYS. REV. X 8, 021021 (2018)

pﬃﬃﬃ ℏ2∂ϕℏ2¼ − pﬃﬃﬃ ∇2 n þð∇ϕÞ2 þ V þ U0 n;
∂t
2m
2m n

ðC2Þ

where V is the potential for the atoms, U0 ¼ 4πℏ2 as =M isthe interaction constant, M is the mass of an atom, and as is
the s-wave scattering length. In this treatment, we neglect
the quantum pressure term [the first term on the right-handside of Eq. (C2)]. By linearizing the equations about
the background solution n0 and ϕ0 , i.e., n ¼ n0 þ n1
and ϕ ¼ ϕ0 þ ϕ1 , one obtains the coupled differential
equations

−ℏ

∂n1
ℏ
¼ − ∇ · ½n0 ∇ϕ1 þ n1 ∇ϕ0 
m
∂t

ðC3Þ

∂ϕ1 ℏ2¼ ∇ϕ0 · ∇ϕ1 þ U0 n1 :
∂t
m

ðC4Þ

(In the main text, we use δn ¼ n1 .) Solving the second
equation for n1 and inserting it into the first one yields awave equation for ϕ1, which, expressed in covariant
form, is1
pﬃﬃﬃﬃﬃﬃ
pﬃﬃﬃﬃﬃﬃ ∂ μ ð −ggμν ∂ ν ϕ1 Þ ¼ 0;
−g

ðC5Þ

where gμν is the inverse metric, g is the metric’s determinant, and ϕ1 is the phonon’s velocity potential field. The
metric, in its most general form, is given by the line elementds2 ¼ gμν dxμ dxν
¼ c0 ½−c20 dt2 þ ðdx − v0 dtÞ · ðdx − v0 dtÞ:

ðC6Þ

field of the backHere, v0 ¼ ðℏ=mÞ∇ϕ0 is the velocity
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ground condensate, and c0 ¼ U0 n0 =M is the speedof sound.
In the expanding ring experiment, the potential is(approximately) axisymmetric and is translated radiallyas a function of time. In terms of cylindrical coordinatesðρ; θ; zÞ, the velocity v0 thus has only ρ and z components,
with its dominant component being radial. The central
radius of the ring is given by a function ρ ¼ RðtÞ, so,assuming the velocity is purely radial, the line element[Eq. (C6)] takes the form
ds2 ¼ c0 ½−c20 dt2 þ ðR þ ρ̃Þ2 dθ2 þ dρ̃2 þ dz2 ;ðC7Þ

where ρ̃ ¼ ρ − RðtÞ is a comoving radial coordinate. The
tensor density that directly enters the wave equation[Eq. (C5)] is then, in ðt; θ; ρ̃; zÞ coordinates,
f μν ≔

pﬃﬃﬃﬃﬃﬃ μν
−gg ¼ ðR þ ρ̃Þdiag½−1; c20 =ðR þ ρ̃Þ2 ; c20 ; c20 :
ðC8Þ

and an equation analogous to the Euler equation,

021021-9

S. ECKEL et al.

PHYS. REV. X 8, 021021 (2018)

The speed c0 is determined by the background condensate
density, which, to a first approximation, follows theinstantaneous ground-state Thomas-Fermi distribution at
all times during the expansion,n0 ¼

μ − Vðρ̃; zÞ
:
U0

ðC9Þ

This density extends out to the contour in the ρ̃ − z planewhere the numerator vanishes. The chemical potential μdrops as the ring expands, so the total number of atoms
remains constant.
In the experiment, we first apply a perturbation to astationary condensate to excite an eigenmode of the wave
equation. An eigenmode analysis based on the methods ofRef. [71] will be detailed in a forthcoming paper; theessential details are presented here. Assuming azimuthal
symmetry, the eigenmodes for a thin ring have the form
ϕ1 ¼ χ klm ηklm ðρ̃; z; RÞeiðωt−mθÞ , where ηklm ðρ̃; z; RÞ is afunction 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Þ ¼
P
−imθ
, with all χ klm ðt ¼ 0Þ ¼ 0
klm χ klm ðtÞηklm ðρ̃; zÞe
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 that occurs when the ring expansion stops, which 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 drop the k andl 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 theseazimuthal phonon modes to a 1 þ 1-dimensional wave
equation, with an effective sound speed cθ . As in Ref. [71],c2θ is given by an average over the cross section of the ring.
For a thin ring, this takes the form
Z
1
2
cθ ¼
½μ − Vðρ̃; zÞdρ̃dz;
ðC10Þ
AM
where the integral is over the cross section of the ThomasFermi wave function of area A. For V ¼ 12 Mω2z z2 þ λρ̃4,this yields4c2θ ¼ c2 ;
7

ðC11Þ

where c2 ¼ μ=M is the peak local sound speed. By
normalizing the Thomas-Fermi solution to the numberof atoms N, one finds that μ ∝ R−4=7 , and thereforecθ ∝ c ∝ R−2=7 .
The wave equation satisfied by our modes of interest,i.e., ϕ1 ¼ χ m ðtÞeimθ , 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
2
f ab
2 ¼ diag½−R; cθ =R:

ðC12Þ

The resulting mode equation isχ̈ m þ

R_
χ_ þ ω2m χ m ¼ 0;
R m

ðC13Þ

where ωm ≔ mcθ =R. This is the equation of a damped
harmonic oscillator, with time-dependent frequency anddamping 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 g2ab for which
pﬃﬃﬃﬃﬃﬃ ab
f ab−gg . The reason for this is that the determinant
2 ¼
pﬃﬃﬃﬃﬃﬃ
of Eq. (C12) is −c2θ , whereas the determinant of −ggab isequal to −1 for any two-dimensional metric.
As the ring expands, the azimuthal wavenumber m isconserved, so the physical wavelength redshifts as R−1 , inanalogy with the cosmological redshift. Unlike in cosmology, the sound speed is also changing, so the frequency ωmredshifts as R−9=7 . In the cosmological setting, the damping
term in Eq. (C13) is called the Hubble friction term andwould 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 ¼ 12 R_χ 2m − 12 ðm2 c2θ =RÞχ 2m , which has the adiabatic invariant Rωm χ 2m . To obtain Eq. (1) in the text, we add thephenomenological damping γ m observed in the experiment.
In the experiment, we measure the density variation n1 ,not the phonon velocity potential ϕ1 . The relation between
these quantities is given by Eq. (C4). Since ∇ϕ1 isazimuthal and ∇ϕ0 is radial, ∇ϕ0 · ∇ϕ1 ¼ 0, so we haven1 ¼ −ℏ _ℏ
ϕ1 ¼ − χ_ m eimθ :
U0
U0

ðC14Þ

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 thisregime. Figure 6 shows such a 2D simulation of BEC in aradially quartic potential, expanding from 10 to 40 μm in
approximately 15 ms with 2 × 105 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 inorder to make the initial Thomas-Fermi width 2 μm. As can

021021-10

A RAPIDLY EXPANDING BOSE-EINSTEIN CONDENSATE: …

FIG. 6. Phonon evolution in the thin ring limit. Blue circles
show 2D GPE simulation data; the red (green) curve is the_solution to Eq. (C13) with (without) the R=R“Hubble friction”
term included. The gray band indicates the time during which the
BEC is inflated; its intensity denotes the expansion velocity
relative to the maximum.

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, theadiabatic limit is satisfied (∂ωm =∂t ≪ ω2m ), and the final
amplitude is accurately predicted using the adiabatic
invariant Rωm χ 2m . While this simulation shows the thinring 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 fromthe 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, Qm ¼ ωm =2γ m , andthe 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.
Figure 7 shows the results of the reduced-χ 2 fit; it shows thevalue of χ 2 vs both Qm¼1 and Qm¼2 in the vicinity of theirbest-fit values for three values of γ H , including the best-fit
value. There are several interesting features. First, thereduced χ 2 is greater than 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 andnonzero annular thickness may play a non-negligible rolein 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

PHYS. REV. X 8, 021021 (2018)

FIG. 7. Goodness of fit χ 2 =ν, where ν is the number of degrees
of freedom, versus quality factors for the m ¼ 1, 2 modes forvarious values of γ H .

change in the minimum of χ 2 with γ H indicates our
uncertainty in γ H . (Part of this insensitivity comes fromour choice, at the time of the experiment, to have n1 ∝ϕ_ ≈ 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.

[1] P. J. E. Peebles, Principles of Physical Cosmology,
Princeton Series in Physics (Princeton University Press,
Princeton, NJ, 1993).
[2] D. Baumann, Inflation, in Physics of the Large and the
Small (World Scientific, Singapore, 2012), pp. 523–686.
[3] L. Kofman, A. Linde, and A. A. Starobinsky, Reheating
After Inflation, Phys. Rev. Lett. 73, 3195 (1994).
[4] G. Felder and L. Kofman, Development of Equilibrium After
Preheating, Phys. Rev. D 63, 103503 (2001).
[5] L. Kofman, Tachyonic Preheating, in Proceedings of theEighth International Conference, Particles, Strings, andCosmology, edited by P. Frampton and J. Ng (Rinton Press,
Paramus, 2001), p. 167.
[6] E. J. Copeland, S. Pascoli, and A. Rajantie, Dynamics of
Tachyonic Preheating After Hybrid Inflation, Phys. Rev. D
65, 103517 (2002).
[7] M. Gleiser, N. Graham, and N. Stamatopoulos, Long-Lived
Time-Dependent Remnants During Cosmological Symmetry
Breaking: From Inflation to the Electroweak Scale, Phys.
Rev. D 82, 043517 (2010).
[8] M. A. Amin, R. Easther, H. Finkel, R. Flauger, and M. P.
Hertzberg, Oscillons After Inflation, Phys. Rev. Lett. 108,241302 (2012).

021021-11

S. ECKEL et al.

PHYS. REV. X 8, 021021 (2018)[9] M. A. Amin, M. P. Hertzberg, D. I. Kaiser, and J. Karouby,
Nonperturbative Dynamics of Reheating After Inflation: A
Review, Int. J. Mod. Phys. D 24, 1530003 (2015).
[10] L. P. Pitaevskii and S. Stringari, Bose-Einstein Condensation, International Series of Monographs on Physics(Clarendon Press, Oxford, 2003).
[11] That is, neglecting the quantum pressure term.
[12] W. G. Unruh, Experimental Black-Hole Evaporation? Phys.
Rev. Lett. 46, 1351 (1981).
[13] U. R. Fischer and M. Visser, Riemannian Geometry of
Irrotational Vortex Acoustics, Phys. Rev. Lett. 88,
110201 (2002).
[14] C. Barceló, S. Liberati, and M. Visser, Analogue Gravity,
Living Rev. Relativity 8, 12 (2005).
[15] L. J. Garay, J. R. Anglin, J. I. Cirac, and P. Zoller, Sonic
Analog of Gravitational Black Holes in Bose-Einstein
Condensates, Phys. Rev. Lett. 85, 4643 (2000).
[16] L. J. Garay, J. R. Anglin, J. I. Cirac, and P. Zoller, Sonic
Black Holes in Dilute Bose-Einstein Condensates, Phys.
Rev. A 63, 023611 (2001).
[17] R. Schützhold, Effective Horizons in the Laboratory, in
Quantum Analogues: From Phase Transitions to Black
Holes and Cosmology, Lecture Notes in Physics, Vol. 718,
edited by R. Schützhold and W. G. Unruh (Springer, Berlin,
Heidelberg, 2007), pp. 5–28.
[18] P. Jain, A. S. Bradley, and C. W. Gardiner, Quantum deLaval Nozzle: Stability and Quantum Dynamics of Sonic
Horizons in a Toroidally Trapped Bose Gas Containing a
Superflow, Phys. Rev. A 76, 023617 (2007).
[19] J. Macher and R. Parentani, Black-Hole Radiation in
Bose-Einstein Condensates, Phys. Rev. A 80, 043601(2009).
[20] O. Lahav, A. Itah, A. Blumkin, C. Gordon, S. Rinott,
A. Zayats, and J. Steinhauer, Realization of a Sonic Black
Hole Analog in a Bose-Einstein Condensate, Phys. Rev.
Lett. 105, 240401 (2010).
[21] J. Steinhauer, Observation of Self-Amplifying Hawking
Radiation in an Analogue Black-Hole Laser, Nat. Phys.
10, 864 (2014).
[22] J. Steinhauer, Observation of Quantum Hawking Radiation
and Its Entanglement in an Analogue Black Hole, Nat. Phys.
12, 959 (2016).
[23] P. O. Fedichev and U. R. Fischer, Gibbons-Hawking Effect in
the Sonic de Sitter Space-Time of an Expanding BoseEinstein-Condensed Gas, Phys. Rev. Lett. 91, 240407 (2003).
[24] C. Barceló, S. Liberati, and M. Visser, Probing Semiclassical Analog Gravity in Bose-Einstein Condensateswith Widely Tunable Interactions, Phys. Rev. A 68,
053613 (2003).
[25] P. O. Fedichev and U. R. Fischer, “Cosmological” Quasiparticle Production in Harmonically Trapped Superfluid
Gases, Phys. Rev. A 69, 033602 (2004).
[26] U. R. Fischer and R. Schützhold, Quantum Simulation of
Cosmic Inflation in Two-Component Bose-Einstein Condensates, Phys. Rev. A 70, 063615 (2004).
[27] E. A. Calzetta and B. L. Hu, Early Universe Quantum
Processes in BEC Collapse Experiments, Int. J. Theor.
Phys. 44, 1691 (2005).
[28] P. Jain, S. Weinfurtner, M. Visser, and C. W. Gardiner,
Analog Model of a Friedmann-Robertson-Walker Universe

in Bose-Einstein Condensates: Application of the ClassicalField Method, Phys. Rev. A 76, 033616 (2007).
[29] A. Prain, S. Fagnocchi, and S. Liberati, Analogue Cosmological Particle Creation: Quantum Correlations in
Expanding Bose-Einstein Condensates, Phys. Rev. D 82,
105018 (2010).
[30] A. D. Sakharov, The Initial Stage of an Expanding Universe
and the Appearance of a Nonuniform Distribution of
Matter, Zh. Eksp. Theor. Fiz. 49, 345 (1965) [English
translation by J. G. Adashko, Sov. Phys. JETP 22, 241(1966)].
[31] G. T. Moore, Quantum Theory of the Electromagnetic Field
in a Variable-Length One-Dimensional Cavity, J. Math.
Phys. (N.Y.) 11, 2679 (1970).
[32] C.-L. Hung, V. Gurarie, and C. Chin, From Cosmology
to Cold Atoms: Observation of Sakharov Oscillations in
a Quenched Atomic Superfluid, Science 341, 1213
(2013).
[33] J.-C. Jaskula, G. B. Partridge, M. Bonneau, R. Lopes, J.
Ruaudel, D. Boiron, and C. I. Westbrook, Acoustic Analog
to the Dynamical Casimir Effect in a Bose-Einstein
Condensate, Phys. Rev. Lett. 109, 220401 (2012).
[34] F. Belgiorno, S. L. Cacciatori, M. Clerici, V. Gorini, G.
Ortenzi, L. Rizzi, E. Rubino, V. G. Sala, and D. Faccio,
Hawking Radiation from Ultrashort Laser Pulse Filaments,
Phys. Rev. Lett. 105, 203901 (2010).
[35] E. Rubino, J. McLenaghan, S. C. Kehr, F. Belgiorno,
D. Townsend, S. Rohr, C. E. Kuklewicz, U. Leonhardt,
F. König, and D. Faccio, Negative-Frequency Resonant
Radiation, Phys. Rev. Lett. 108, 253901 (2012).
[36] D. Vocke, C. Maitland, A. Prain, F. Biancalana, F. Marino,
E. M. Wright, and D. Faccio, Rotating Black Hole
Geometries in a Two-Dimensional Photon Superfluid,
arXiv:1709.04293.
[37] S. Weinfurtner, E. W. Tedford, M. C. J. Penrice, W. G.
Unruh, and G. A. Lawrence, Measurement of Stimulated
Hawking Emission in an Analogue System, Phys. Rev. Lett.
106, 021302 (2011).
[38] L.-P. Euvé, F. Michel, R. Parentani, T. G. Philbin, and G.
Rousseaux, Observation of Noise Correlated by the Hawking
Effect in a Water Tank, Phys. Rev. Lett. 117, 121301(2016).
[39] T. Torres, S. Patrick, A. Coutant, M. Richartz, E. W.
Tedford, and S. Weinfurtner, Rotational SuperradiantScattering in a Vortex Flow, Nat. Phys. 13, 833 (2017).
[40] D. Faccio, F. Belgiorno, S. Cacciatori, V. Gorini, S. Liberati,
and U. Moschella, Analogue Gravity Phenomenology:
Analogue Spacetimes and Horizons, from Theory to Experiment, Lecture Notes in Physics (Springer International
Publishing, New York, 2013).
[41] With perfect imaging, our potential would be flat bottomed,
but a combination of finite resolution and aberrations causes
our potential to have both curvature and broken axial
symmetry.
[42] S. J. Rooney, P. B. Blakie, and A. S. Bradley, Stochastic
Projected Gross-Pitaevskii Equation, Phys. Rev. A 86,
053634 (2012).
[43] A. S. Bradley, S. J. Rooney, and R. G. McDonald,
Low-Dimensional Stochastic Projected Gross-Pitaevskii
Equation, Phys. Rev. A 92, 033631 (2015).

021021-12

A RAPIDLY EXPANDING BOSE-EINSTEIN CONDENSATE: …
[44] Unless stated otherwise, all uncertainties and error bars inthis paper are the uncorrelated combination of 1σ statisticaland systematic uncertainties.
[45] We measure Q ¼ γ m =ωm ≈ 6 for all m > 1 and Q ≳ 20 for
m ¼ 1, indicating that the intrinsic damping γ m is most
likely Landau damping [46]. Note that Q decreases slightly
for increasing m and R. For fitting a given m, we take Q to
be constant with time.
[46] L. P. Pitaevskii and S. Stringari, Landau Damping in Dilute
Bose Gases, Phys. Lett. A 235, 398 (1997).
[47] C. J. Pethick and H. Smith, Bose-Einstein Condensation in
Dilute Gases (Cambridge University Press, Cambridge,
England, 2002).
[48] D. J. Frantzeskakis, Dark Solitons in Atomic Bose-Einstein
Condensates: From Theory to Experiments, J. Phys. A 43,
213001 (2010).
[49] R. G. Scott, A. M. Martin, T. M. Fromhold, S. Bujkiewicz,
F. W. Sheard, and M. Leadbeater, Creation of Solitons andVortices by Bragg Reflection of Bose-Einstein Condensatesin an Optical Lattice, Phys. Rev. Lett. 90, 110404 (2003).
[50] A. M. Martin, R. G. Scott, and T. M. Fromhold,Transmission and Reflection of Bose-Einstein Condensates
Incident on a Gaussian Tunnel Barrier, Phys. Rev. A 75,
065602 (2007).
[51] For the expansion described by Eq. (A2), χ r ¼ ðRf − Ri Þ×
exp½−ðβωr texp Þ2 =4. Thus, we must compare the energy ofthe soliton to that of the energy in the radial mode, 12 mω2r χ 2r .
The full adiabatic condition then becomes vp ≳cðRf −Ri Þ=
pﬃﬃﬃ
RT × π ½log(3π 2 mcðRf −Ri Þ2 =ℏRT )−1=2 .
[52] E. A. Kuznetsov and S. K. Turitsyn, Instability and Collapseof Solitons in Media with a Defocusing Nonlinearity,
Zh. Eksp. Theor. Fiz. 94, 119 (1988) [English translation
by A. Tybulewicz, Sov. Phys. JETP 67, 1583 (1988)].
[53] A. M. Kamchatnov and L. P. Pitaevskii, Stabilization of
Solitons Generated by a Supersonic Flow of Bose-Einstein
Condensate Past an Obstacle, Phys. Rev. Lett. 100, 160402(2008).
[54] A. Cetoli, J. Brand, R. G. Scott, F. Dalfovo, and L. P.
Pitaevskii, Snake Instability of Dark Solitons in FermionicSuperfluids, Phys. Rev. A 88, 043639 (2013).
[55] L. A. Toikka and K. A. Suominen, Snake Instability ofRing Dark Solitons in Toroidally Trapped Bose-Einstein
Condensates, Phys. Rev. A 87, 043601 (2013).
[56] A. Ramanathan, S. R. Muniz, K. C. Wright, R. P. Anderson,
W. D. Phillips, K. Helmerson, and G. K. Campbell, PartialTransfer Absorption Imaging: A Versatile Technique for
Optimal Imaging of Ultracold Gases, Rev. Sci. Instrum. 83,
083119 (2012).
[57] In the simulations, Sðkθ Þ as computed from the density(as opposed to simulated, aberrated images) is approximately 1 at these early times. This corresponds to the

PHYS. REV. X 8, 021021 (2018)

equilibrium value for a condensate where kB T ≈ μ, T is thetemperature, and kB is Boltzmann’s constant.
[58] R. H. Kraichnan, Inertial Ranges in Two-Dimensional
Turbulence, Phys. Fluids 10, 1417 (1967).
[59] P. W. Anderson, Considerations on the Flow of Superfluid
Helium, Rev. Mod. Phys. 38, 298 (1966).
[60] The approximation of constant speed of sound is good towithin 10% because of the weak dependence of c on R.
[61] W. H. Zurek, Cosmological Experiments in Superfluid
Helium?, Nature (London) 317, 505 (1985).
[62] D. R. Scherer, C. N. Weiler, T. W. Neely, and B. P.
Anderson, Vortex Formation by Merging of Multiple
Trapped Bose-Einstein Condensates, Phys. Rev. Lett. 98,
110402 (2007).
[63] M. J. Bowick, L. Chandar, E. A. Schiff, and A. M.
Srivastava, The Cosmological Kibble Mechanism in the
Laboratory: String Formation in Liquid Crystals. Science
263, 943 (1994).
[64] We find, through Monte Carlo simulations, that the probabilitiespare
ﬃﬃﬃﬃﬃﬃﬃ well described by normal distributions withσ ≈ 0.3 N R .
[65] A. Kumar, S. Eckel, F. Jendrzejewski, and G. K. Campbell,Temperature-Induced Decay of Persistent Currents in a
Superfluid Ultracold Gas, Phys. Rev. A 95, 021602(2017).
[66] J. G. Lee and W. T. Hill, Spatial Shaping for Generating
Arbitrary Optical Dipole Traps for Ultracold Degenerate
Gases, Rev. Sci. Instrum. 85, 103106 (2014).
[67] A. Kumar, N. Anderson, W. D. Phillips, S. Eckel,
G. K. Campbell, and S. Stringari, Minimally Destructive,
Doppler Measurement of a Quantized Flow in a RingShaped Bose-Einstein Condensate, New J. Phys. 18,025001 (2016).
[68] S. Eckel, J. G. Lee, F. Jendrzejewski, N. Murray,
C. W. Clark, C. J. Lobb, W. D. Phillips, M. Edwards,
and G. K. Campbell, Hysteresis in a Quantized Superfluid “Atomtronic” Circuit, Nature (London) 506, 200(2014).
[69] L. Corman, L. Chomaz, T. Bienaimé, R. Desbuquois, C.
Weitenberg, S. Nascimbène, J. Dalibard, and J. Beugnon,
Quench-Induced Supercurrents in an Annular Bose Gas,
Phys. Rev. Lett. 113, 135302 (2014).
[70] C.-L. Hung, X. Zhang, L.-C. Ha, S.-K. Tung, N. Gemelke,
and C. Chin, Extracting Density-Density Correlations fromIn Situ Images of Atomic Quantum Gases, New J. Phys. 13,075019 (2011).
[71] E. Zaremba, Sound Propagation in a Cylindrical BoseCondensed Gas, Phys. Rev. A 57, 518 (1998).
[72] Our experiments were performed before we realized thesignificance of the Hubble friction effect, which was
only later pointed out by our more enlightened theory
co-author.

021021-13