# On the survival of the quantum depletion of a condensate after release from
a magnetic trap

J. A. Ross  Research School of Physics, Australian National University,
Canberra 0200, Australia  P. Deuar  Institute of Physics, Polish Academy of
Sciences, Aleja Lotników 32/46, 02-688 Warsaw, Poland  D. K. Shin  Research
School of Physics, Australian National University, Canberra 0200, Australia
K. F. Thomas  Research School of Physics, Australian National University,
Canberra 0200, Australia  B. M. Henson  Research School of Physics, Australian
National University, Canberra 0200, Australia  S. S. Hodgman  Research School
of Physics, Australian National University, Canberra 0200, Australia  A. G.
Truscott

###### Abstract

We present observations of the high momentum tail in expanding Bose-Einstein
condensates of metastable Helium atoms released from a harmonic trap. The far-
field density profile exhibits features that support identification of the
tails of the momentum distribution as originating in the in-situ quantum
depletion prior to release. Thus, we corroborate recent observations of
slowly-decaying tails in the far-field beyond the thermal component. This
observation is in conflict with the hydrodynamic theory, which predicts that
the in-situ depletion does not survive when atoms are released from a trap.
Indeed, the depleted tails even appear stronger in the far-field than expected
before release, and we discuss the challenges of interpreting this in terms of
the Tan contact in the trapped gas. In complement to these observations, full
quantum simulations of the experiment show that, under the right conditions,
the depletion can persist into the far field after expansion. Moreover, the
simulations provide mechanisms for survival and for the the large-momentum
tails to appear stronger after expansion due to an acceleration of the
depleted atoms by the mean-field potential. However, while in qualitative
agreement, the final depletion observed in the experiment is much larger than
in the simulation.

## Introduction

In the Bogoliubov description of an ultracold interacting superfluid, the
ground state is composed of a macroscopically-occupied condensate and
correlated particle pairs due to s-wave interactions between constituent
particles [1, 2]. A consequence of these pairs is that excited single-particle
modes are populated even at zero temperature. This is the _quantum depletion_
of the condensate and presents as an occupation of single particle modes,
which at large momentum $p$ decays [3, 4] like $p^{-4}$.

Since the realization of atomic Bose-Einstein condensates (BECs) there has
been considerable experimental [5, 6, 7, 8, 9, 10, 2, 11, 12, 13, 14, 15, 16,
17, 18, 19] and theoretical [20, 21, 4, 22, 23, 24, 25, 26, 27, 28, 29, 30,
31, 32, 33] interest in the Bogoliubov theory [2, 34, 17, 32, 33] (and quantum
depletion specifically [12, 7, 10, 11, 13, 35]). In contrast to the case of
liquid helium, where the depleted fraction is large (of order 93% of the fluid
[36, 37, 38]) due to the strong interparticle interactions, the depletion is
generally very small (less than 1% [12, 7]) in weakly-interacting dilute
gases. The intimately related thermodynamic (Tan’s) contact has also received
growing attention, [5, 39, 40, 41, 24, 25, 20, 8, 9, 4, 22, 7, 23, 6, 42, 26,
27, 22, 14, 15, 16, 18, 19, 28, 29, 30, 31, 35], in part due to Tan’s proof
that the contact is directly related to the amplitude of the $p^{-4}$ tail
[39].

Experiments examining the large-momentum tails have typically employed
Feshbach resonances to enhance interactions in ultracold gases and produce a
depleted fraction visible in the far-field with standard imaging techniques,
but the power-law tails have proven elusive [8, 9] in this regime. A handful
of theories have emerged [21, 20, 22] which elucidate the role played by many-
body interactions in modifying the momentum distribution during the evolution
following a quench to a large scattering length. A very recent experiment [43]
was able to detect pairs of atoms with anticorrelated momenta in the far-field
by use of an optical lattice to create a BEC in a high-density, strongly-
interacting regime. However, measurements in the weakly-interacting regime
have returned unexpected results. A previous experiment reported the presence
of power-law-like tails in the far-field distribution after releasing a BEC of
metastable helium from a harmonic optical trap [7]. This was surprising
because conventional wisdom argues that the density decreases adiabatically
during expansion (even when the trap release is non-adiabatic), motivating a
hydrodynamic approximation wherein the tails are predicted to vanish [23, 10].
Moreover, the tails were reported to be approximately six times heavier than
predicted by Bogoliubov theory. It is important to verify the anomaly and
understand its origin because far-field measurements play a central role in
the study of ultracold gases. The prospect of extracting correlated depleted
pairs from a zero-temperature ground state is also conceptually, and possibly
technologically, interesting in itself.

To these ends, we measure the far-field momentum distribution of a BEC of
metastable helium (He${}^{*}\leavevmode\nobreak\ $) expanding from a harmonic
trap. We observe tails in the large-momentum part of the (far-field)
condensate wavefunction whose amplitude depends nonlinearly on the condensate
population, and whose density profile is consistent with a $p^{-4}$-like power
law decay, in a manner consistent with the predictions of the Tan and
Bogoliubov theory. Specifically, the amplitude of the far-field momentum tails
is shown to have a linear relationship with the product of the condensate
population and peak density, as predicted by both theories. However, there is
a quantitative difference in amplitude between the predicted and measured
values. Our measurements are complemented by numerical simulations of the
dynamics of the momentum distribution after the trap release using a
Stochastic Time-Adaptive Bogoliubov (STAB) method in the positive-P framework
[33, 44]. These demonstrate a mechanism for survival associated with the non-
adiabatic release of the trap, and suggest that the depleted particles acquire
additional kinetic energy from the mean-field energy of the condensate during
the subsequent adiabatic expansion. These factors result in an amplification
of the density of the far-field momentum tails relative to the in-situ values
by a factor of up to about two, and are absent from the hydrodynamic
approximation. However, even taking these effects into account, the amplitude
of the measured tails is still significantly larger than expected from the
simulations.

### Quantum depletion of the condensate by contact interactions

Before presenting our results, let us introduce the central theoretical
assumptions and predictions relevant for this work. The Hamiltonian of a
homogeneous system of interacting bosons can be written in terms of plane-wave
field operators $\hat{a_{\textbf{k}}}$, labeled by the wavevector
$\textbf{k}=\textbf{p}/\hbar$, and diagonalized by the Bogoliubov
transformation to a free Bose gas of collective excitations through the
operator transformation
$\hat{b}_{\textbf{k}}^{\dagger}=u_{k}\hat{a}_{\textbf{k}}^{\dagger}+v_{k}\hat{a}_{-\textbf{k}}$
[1, 45]. The collective excitations are superpositions of particles with
opposite momenta [2], and the $u_{k}$ and $v_{k}$ coefficients are given by

$\displaystyle u_{k}$ $\displaystyle=\cosh\theta_{k},\qquad
v_{k}=\sinh\theta_{k}$ (1) $\displaystyle\theta_{k}$
$\displaystyle=\frac{1}{2}\log\frac{\hbar^{2}k^{2}/2m}{\epsilon(k)}<0$ (2)

where the denominator is the quasiparticle dispersion

$\epsilon(k)=\sqrt{\left(\frac{\hbar^{2}k^{2}}{2m}\right)^{2}+gn\frac{\hbar^{2}k^{2}}{m}}.$
(3)

determined by the particle density $n$, the atomic mass $m$, and the effective
interaction strength $g=4\pi\hbar^{2}a/m$, where $a$ is the s-wave scattering
length [3, 45]. In the non-interacting ($a\rightarrow 0$) limit, $u_{k}=1$ and
$v_{k}=0$, so the transformation reduces to the identity and the dispersion is
that of free particles. The occupation of single-particle momentum modes can
be found using the inverse transformation and is given by

$\displaystyle\rho(\textbf{k})$
$\displaystyle=\langle\hat{a}_{\textbf{k}}^{\dagger}\hat{a}_{\textbf{k}}\rangle$
(4)
$\displaystyle=\left(u_{k}^{2}+v_{k}^{2}\right)\langle\hat{b}_{\textbf{k}}^{\dagger}\hat{b}_{\textbf{k}}\rangle+v_{k}^{2},$
(5)

wherein the quasiparticle population statistics follow the canonical ensemble
as [3, 7]
$\langle\hat{b}^{\dagger}_{\textbf{k}}\hat{b}_{\textbf{k}}\rangle=(\exp[\epsilon(k)/k_{B}T]-1)^{-1}$.
At finite temperatures, quasiparticle modes are thermally populated and
deplete the condensate. Even at zero temperature, when the thermal fraction
vanishes, the $v_{k}^{2}$ term in Eqn (5) persists giving a zero-temperature
population of excited particles [46, 4, 7] which decays as [45, 3, 7]
$\lim_{k\rightarrow\infty}\rho(\textbf{k})\propto k^{-4}$. Bogoliubov’s theory
makes accurate predictions of the total depleted population in ultracold
atomic Bose-Einstein condensates (BECs) [10, 12] and exciton-polariton
condensates in solid substrates [11].

In the case of a harmonically trapped gas, one can employ the local-density
approximation (LDA) to compute the amplitude of the $k^{-4}$ tail by
integrating $v_{k}^{2}$ across a Thomas-Fermi distribution [7]. One can also
compute the expected amplitude of the tails using thermodynamic relations
between the condensate mean-field energy and momentum distribution: The
amplitude of the tails was shown by Tan to be exactly the quantity called the
_contact_ , which is proportional to the derivative of the energy with respect
to the s-wave scattering length [39, 25]. For a Bose gas at equilibrium in a
harmonic trap, the tail amplitude can be calculated using Tan’s original
theorems. The two-body _contact intensity_ is defined by [39, 25]

$C=\lim_{k\rightarrow\infty}k^{4}\rho(k),$ (6)

which is related to the total contact (or just _contact_) $\mathcal{C}=\int
C(\textbf{r})d^{3}\textbf{r}$. The contact can be derived from the total
energy $E$ through the _adiabatic sweep theorem_ [40],

$\mathcal{C}=\frac{8\pi ma^{2}}{\hbar^{2}}\frac{\partial E}{\partial a}.$ (7)

Applying this to the Thomas-Fermi energy of a harmonically trapped condensate,

$\frac{E}{N_{0}}=\frac{5}{7}\mu=\frac{5}{7}\frac{\hbar\bar{\omega}}{2}\left(\frac{15N_{0}a}{a_{\textrm{HO}}}\right)^{2/5},$
(8)

where $a_{\textrm{HO}}=\sqrt{\hbar/(m\bar{\omega})}$ is the harmonic
oscillator length and $\bar{\omega}=\sqrt[\scriptstyle
3]{\omega_{x}\omega_{y}\omega_{z}}$ is the geometric trapping frequency [3,
45], leads to the expression

$\mathcal{C}=\frac{8\pi}{7}\left(15^{2}(aN_{0})^{7}\left(\frac{m\bar{\omega}}{\hbar}\right)^{6}\right)^{1/5},$
(9)

for the contact and thus the asymptotic momentum (density) distribution $n(k)$
of the in-situ condensate is,

$\lim_{k\rightarrow\infty}n(k)={\frac{\mathcal{C}}{k^{4}}}=\frac{64\pi^{2}a^{2}}{7}\frac{N_{0}n_{0}}{k^{4}},$
(10)

which depends on the peak density of the harmonically trapped condensate, in
turn given by

$n_{0}=\frac{1}{8\pi}\left((15N_{0})^{2}\left(\frac{m\bar{\omega}}{\hbar\sqrt{a}}\right)^{6}\right)^{1/5}.$
(11)

Note that hereon we refer to the momentum distribution $n(\textbf{k})$ rather
than the occupation numbers $\rho(\textbf{k})=n(k)d^{3}\textbf{k}/(2\pi)^{3}$,
and that the total number of atoms in this normalisation is
$N=\frac{1}{(2\pi)^{3}}\int d^{3}\textbf{k}\,n(\textbf{k})$.

## Results

### Experimental measurements

Our experimental sequence began with BECs consisting of between $2\times
10^{5}$ and $5\times 10^{5}$ 4He∗ atoms, spin-polarized in the
$2^{3\\!}S_{1}(m_{J}=1)$ state and cooled to $\sim$ 300 nK by forced
evaporative cooling in a harmonic magnetic trap generated by field coils in a
Bi-planar Quadrupole Ioffe configuration [47]. The trap was then switched off
with a $1/e$ time of ${\tau_{\rm release}}\approx 38\mu$s. The condensates
were allowed to expand for 2ms before we transferred about one quarter of the
initial $m_{J}=1$ condensate into the magnetically insensitive $m_{J}=0$ state
with a radio-frequency (RF) Landau-Zener sweep to preserve it against
distortion by stray magnetic fields during the free fall to the detector. We
deflected the $m_{J}=\pm 1$ clouds away from the detector with a Stern-Gerlach
scheme immediately after the RF pulse by switching on a magnetic field. The
centre of mass of the cloud then impacts on the detector after a $\tau=417$ms
time of flight following the trap switch-off.

Investigations of the quantum depletion in He${}^{*}\leavevmode\nobreak\ $are
challenging because the absence of a known Feshbach resonance precludes
control over the contact
$\mathcal{C}\propto((aN_{0})^{7}\bar{\omega}^{6})^{1/5}$ via the scattering
length $a$. Given the small fixed $a=7.512$nm [48], we test the validity of
Eqn. (10) for describing the far-field by varying the density of the gas,
$n\propto\left(N_{0}\bar{\omega}^{3}\right)^{2/5}$ (c.f. Eqn (10)). To achieve
this we used two trap configurations with
$(\omega_{x},\omega_{y},\omega_{z})\approx 2\pi(45,425,425)$ Hz (geometric
mean $\bar{\omega}=2\pi\cdot 201$ Hz) and $\approx 2\pi(71,902,895)$ Hz
($\bar{\omega}=2\pi\cdot 393$ Hz) where the frequencies are known within 1%
the (weak) axis of symmetry is horizontal. We varied the endpoint of the
evaporative cooling ramp to adjust the number of atoms in the condensate.

Our experiment uses single-particle detection with multichannel-plate and
delay-line detector (MCP-DLD) stacks [49] after a long time of flight (hence
in the far-field regime) enabled by the large (19.8eV) internal energy [50] of
the metastable $2^{3\\!}S_{1}$ state, He∗. The unique capabilities of such
setups have permitted the observation of many-body momentum correlations [51,
52] and the Hanbury Brown-Twiss effect in both condensed [53, 54, 49, 55, 56,
57] and quantum depleted atoms [13, 35]. We are thus able to reconstruct the
full single-atom momentum distribution in three dimensions and examine the
dilute far-field momentum tails of the $m_{J}=0$ clouds in detail.

In Figure 1 we show the empirical far-field density $n(k)$ for two data
collection runs at the extreme values of $n_{0}$ we used. The black (purple)
correspond to condensates with an average of $3.5\times 10^{5}(4.5\times
10^{5})$ atoms and a thermal fraction of 9% (10%). The (geometric) trap
frequencies were $2\pi\cdot$201 and $2\pi\cdot 393$ Hz, and the healing length
$\xi=\hbar/\sqrt{2mgn_{0}}$ at the center of these clouds were
$56\leavevmode\nobreak\ \mu\textrm{m}$ and $36\leavevmode\nobreak\
\mu\textrm{m}$, respectively. The three regimes of the condensate, thermal
depletion, and quantum depletion span over five orders of magnitude in
density. The thermal part of the distribution is well fitted by the momentum
distribution of an ideal Bose gas [58]

$\frac{n_{T}(k)}{{(2\pi)^{3}}}=\frac{N_{T}}{\zeta(3)}\leavevmode\nobreak\
\left(\frac{\lambda_{dB}}{2\pi}\right)^{3}g_{3/2}\left(\exp\left(-\frac{k^{2}\lambda_{dB}^{2}}{4\pi}\right)\right)$
(12)

wherein the thermal de Broglie wavelength
$\lambda_{dB}=\sqrt{2\pi\hbar^{2}/(mk_{B}T)}$ yields an estimate of the
temperature $T$ which ranges from 100 to 320 nK in our experiments. Here,
$g_{3/2}(\cdot)$ is the standard Bose integral, $\zeta(\cdot)$ is the Riemann
zeta function, and $N_{T}$ is the number of atoms in the thermal component.
Note that for a non-interacting gas in the thermodynamic limit, the number of
thermal atoms is simply ${N_{T}^{\rm
id}}=\zeta(3)(k_{B}T/\hbar\bar{\omega})^{3}{=\eta_{T}N}$, but for our
condensates the critical temperature is reduced by $\approx 20\%$ by
interactions [3, 45]. We account for this and the approximately twofold
increase in the thermal fraction $\eta_{T}$ (relative to the non-interacting
case) by explicitly using $N_{T}$ as a fit parameter. At larger values of
momentum, where the thermal component makes a negligible contribution, there
appears a slow decay which we identify as the quantum depletion.

Figure 1: The measured far-field density of particle momenta from two trap
configurations (black and magenta). Three regions are shown: At low $k$ the
parabolic distribution of the BEC dominates. For larger $k$, the thermal parts
(fits shown by dashed lines) decay super-exponentially as $e^{-k^{2}}$. For
even larger $k$, these give way to the surmised quantum depletion region. A
combined fit of the form $n_{T}(k)+C_{4}/k^{4}$ (green dot-dash lines) yields
temperatures consistent with the thermal fit and also an amplitude $C_{4}$ of
the depleted tail. The grey dotted line is a guide to the eye showing a
$k^{-4}$ decay. Due to constraints of the detector geometry (see supplementary
materials for details), these profiles were integrated over two spherical
segments, each subtending an angle of $\pi/6$ radians with the $\pm z$ axes
(see also Figure 5). The detector shows signs of saturation for low
$k(\lesssim 1.5\leavevmode\nobreak\ \mu\textrm{m}^{-1}$). These factors imply
the total area under the curves is less than the average number of trapped
atoms.

### Analysis of experimental results

A standard approach to analysing the empirical momentum density would be to
proceed with a routine fit of the $k$-space histogram with an additional term
of the form $C_{\alpha}/k^{\alpha}$ to estimate the parameters of the
purported quantum-depleted tail. If we augment the thermal fit function (Eqn.
(12)) with a power-law term as per

${n(k)}=n_{T}(k)+\frac{C_{\alpha}}{k^{\alpha}}$ (13)

and leave $\alpha$ as a free parameter, the average exponent over all runs is
4.2(4). For comparison, the prior work [7] reported power-law tails with an
exponent 4.2(2). At first glance, one could simply determine the amplitude of
the tails by fixing the exponent to 4, and if we do so, we find an average
$C_{\alpha=4}$ which is approximately 8(2) times greater than the coefficient
predicted by Eqn. (10), and in general agreement with Ref. [7]. However, as we
detail in the supplementary materials, the covariance of the fit parameters
$C$ and $\alpha$, coupled with the exponential relationship to the independent
variable $k$, means that this gives a significant underestimate of the
uncertainty in $C_{\alpha}$. In general, fitting power laws to data is known
to be prone to return biased estimates of parameters and to drastically under-
report uncertainties,[59, 60] especially when data is available over less than
a couple of decades of dynamic range.

Below we present a number of lines of evidence that support the identification
of these tails as originating in the quantum depletion, but we also argue
there is not sufficient reason to assume that the fit with a fixed $\alpha=4$
is appropriate. The main reason for the latter is that the far-field momentum
distribution is known to be a modification of the in-situ distribution due to
the dispersal of the condensate mean-field energy into kinetic energy. Even
neglecting this effect, it is not a given that the far-field distribution
could be modified in such a way as to simply increase the amplitude of the
tails without otherwise altering the functional form (i.e. the exponent in
this case). As the authors of Ref. [59] note, “In practice, we can rarely, if
ever, be certain that an observed quantity is drawn from a power-law
distribution. The most we can say is that our observations are consistent with
the hypothesis that $x$ is drawn from […] a power law". Indeed, this analysis
does show that the far-field momentum distribution is consistent with a power-
law exponent $3.8\leq\alpha\leq 4.6$, but the data at hand cannot precisely
determine the exponent $\alpha$ (nor $C_{\alpha}$), as detailed in the
supplement.

Figure 2: Population of momentum tails, including excess compared to Tan-
Bogoliubov theory. (a) The product $N_{0}n_{0}$ is a linear predictor of the
number of counts within the region $(k_{\textrm{min}}=6\leavevmode\nobreak\
\mu\textrm{m}^{-1},k_{\textrm{max}}=10\leavevmode\nobreak\
\mu\textrm{m}^{-1})$, consistent with Eqn (10) (solid orange line, dashed
lines 95% CI). The gradient $\Lambda$ in Eqn. (15) can be predicted using Eqn
(10) ($\Lambda_{\rm pred}$ solid purple line) but this disagrees with the
experiment by a factor of about 8. Our simulations (dashed line, CE in Fig.
3a) show an increase in counts after release but by less than in the
experiment. In (b,c) linear fits to the experimental data yield $\Lambda_{\rm
fit}$ (points) which vary with the choice of $k$ bounds (fixing
$k_{\textrm{max}}=10\mu\textrm{m}^{-1}$ in (b) and
$k_{\textrm{min}}=6\mu\textrm{m}^{-1}$ in (c))w. For comparison, we show
predictions of $\Lambda$ based directly on Eqn (10) ($\Lambda_{\rm pred}$,
blue, $n(k)=\mathcal{C}/k^{4}$), along with the predictions from Eqn. (15)
using a density function $n(k){=\mathcal{AC}/k^{4}}$ that has an additional
prefactor $\mathcal{A}=$8(3) (green) and one that has a modified exponent of
$\alpha=3.86(2)$ via $n(k)=\mathcal{C}/k^{\alpha}$ (yellow). A log-normal
distribution produces nearly identical predictions (red, offset vertically for
visibility). Quoted error estimates correspond to 95% CI of the fit
parameters. In (b), the deviation from the predictions at
$k_{\textrm{min}}\lesssim 6\leavevmode\nobreak\ \mu\textrm{m}^{-1}$ is because
the collection area starts to overlap with the thermal cloud.

Rather than explicitly enforce a power law decay assumption, we focus on
another observable which can be readily measured and predicted: The number of
atoms whose wavevector has a modulus in the interval
$k\in(k_{\textrm{min}},k_{\textrm{max}})$,

$N_{k_{\textrm{min}},k_{\textrm{max}}}=\frac{\mathcal{C}}{2\pi^{2}}\left(\frac{1}{k_{\textrm{min}}}-\frac{1}{k_{\textrm{max}}}\right)$
(14)

Note that the integral of $n(k)$ is most easily performed in spherical
coordinates and requires the Jacobian $(2\pi)^{-3}{d^{3}}\textbf{k}$ to ensure
normalization. For fixed $k_{\textrm{min}}$ and $k_{\textrm{max}}$, Eqn. (14)
has the form

$N_{k_{\textrm{min}},k_{\textrm{max}}}=\Lambda N_{0}n_{0}$ (15)

(c.f. Eqn. (10)). We can thus test Eqn. (15) directly by measuring the number
of counts detected in the interval $(k_{\textrm{min}},k_{\textrm{max}})$ after
producing a BEC of $N_{0}$ atoms with peak density $n_{0}$. A key advantage of
this method is that theoretical assumptions (such as the exponent of the power
law) are not required when analysing the experimental data, but only when
calculating the (independent) prediction, i.e. the data processing is
essentially theory-free.

Under the null hypothesis (based on the hydrodynamic theory) that the _in
situ_ depletion does not survive the expansion, $\Lambda=0$. Further, most
types of technical noise masquerading as high energy tails would be expected
to not follow the $N_{0}n_{0}$ scaling and give at best a poor correlation
with Eq. (15). As we show in Fig. 2, a linear fit of the form
$\hat{N}_{k_{\textrm{min}},k_{\textrm{max}}}=\Lambda_{\textrm{fit}}n_{0}N_{0}+\beta$
yields an intercept consistent with zero ($\beta$=-0.9, 95% CI (-3.1, 1.2))
and a good correlation ($r^{2}\approx 0.8$, $p=1\times 10^{-3}$), providing
evidence supporting the expected linear relationship with $n_{0}N_{0}$, and
against the high energy tails being due to some technical noise. The
correlation coefficient between the variables $N_{k_{\rm min},k_{\rm max}}$
and $N_{0}n_{0}\propto(N_{0}^{7}\bar{\omega}^{6})^{1/5}$ is 0.9. We conclude
that the product $N_{0}n_{0}$ is a predictor of the high energy population,
which is consistent with Eqn. (10).

For comparison, a linear fit proves that the atom number $N_{0}$ itself is a
poor predictor of the detected number ($r^{2}=0.05\leavevmode\nobreak\
,p=0.54$), as is the central density $n_{0}$ alone
($r^{2}=0.4\leavevmode\nobreak\ ,p=0.04$). Accordingly, the particular
nonlinear scaling of detected counts with the predictor $N_{0}n_{0}$ is
concordant with the tails’ originating in the quantum depletion, and
inconsistent with any technical noise that we know of.

The gradient $\Lambda_{\textrm{fit}}$ is of particular interest because it can
be predicted using Eqn. (14). Given a region of interest (ROI) over which we
count atoms, one can calculate $\Lambda_{\textrm{pred}}=32\epsilon
a^{2}(k_{\textrm{min}}^{-1}-k_{\textrm{max}}^{-1})/7$, where $\epsilon$ is the
total detection efficiency. In our experiment, $\epsilon\approx 0.23(5)\%$
(see Methods for details). In comparison with the predicted value
$\Lambda_{\mathrm{pred}}=2.7(6)\times 10^{-7}$ (units of
$\mu\textrm{m}^{3}$/atom), we find that the the empirical fit disagrees with
the predicted slope by a factor of
$\mathcal{A}_{\mathrm{exp}}=\Lambda_{\textrm{fit}}/\Lambda_{\textrm{pred}}=8.3$,
95% CI $(5.5,11)$, which rules out the null hypothesis.

While this result may appear to restate the previously-mentioned fitting
approach, which gave an increase of the $C_{4}$ coefficient by a factor of
8(2), it in fact complements it. In this case the overpopulation of the tails
is directly measured without any recourse to assumptions about power-law
behaviour in the data itself. The direct comparison of the populations in a
given $k$-interval allows for an independent comparison between the prediction
and measured result and seeks to simply answer the question does the data
satisfy the most general model of quantum depletion proposed by Eq. (15).

In summary, there are three robust conclusions that can be drawn from the
data. _First_ , the population in the high-momentum tails depends linearly on
the product $n_{0}N_{0}$, which is a prediction of the Tan and Bogoliubov
theories and not readily associated with any other known physical process.
_Second_ , there are some 8(3) times as many particles in the far-field, high-
momentum tails as would be expected to be found in the same interval of the
in-situ distribution. _Third_ , the data is consistent with power-laws with
exponents in the range $3.8\leq\alpha\leq 4.6$.

Interestingly, if one were to take the first observation as sufficient (and
indeed independent) evidence to identify the tails with the quantum depletion
_and_ assume a power-law decay of the form $C_{4}k^{-4}$, then one obtains a
value of $C_{4}$ that is consistent with the regression against $n_{0}N_{0}$.
While this is evidence that the $\alpha=4$ hypothesis is not inconsistent with
the data, it is essentially the same as calculating $C$ from the results of
the linear regression by assuming $\alpha=4$. In Figure 2 (b,c), as a
counterpoint to power law fits over $k_{\rm min,max}$ shown in blue, green and
yellow, we also show, in red, predictions obtained by assuming log-normally
distributed $k$ with parameters $(\mu,\sigma)\approx(1.235,0.95)$ and
normalized to the relevant amplitude. This underscores the challenge of
identifying power-law behaviour in range-limited data, because although the
log-normal distribution eventually diverges from the power law, it does so
over a much larger domain than available in either Helium experiment (here or
[7]). These fits scarcely differ in their goodness-of-fit criterion (the mean
square error) and so offer no obvious way to reconcile the expected
distribution with these divergent statistical conclusions.

### Findings from simulations

Figure 3: Simulations of release from the trap. (a) Steady-state values of the
simulated contact. Simulations of condensates released from a cigar-shaped
trap (CT) are consistent with the Tan theory (TT) before release, and show an
increase in contact after the trap release (CE). A slow relaxation of the
transverse trapping frequencies (CS) shows a decrease in line with the
predicted value of the lower density. Spherical traps (ST,SE) lack any
directions of tight confinement, wherein a longer interaction time prevents
the escape of depleted particles as seen in cigar traps. (b) the time-
dependence of the contact stabilizes after a time on the order of
$1/\omega_{x}$, several hundred $\mu$s. The expanded contact is consistently
about 1.7 times the Tan theory. For comparison, the experimental control
pulses are implemented after 2ms of expansion. When the transverse trapping
frequencies are slowly (1.2ms) reduced by half (dotted line), the in-situ
contact relaxes on a faster timescale than the ramp.

In order to understand whether the depletion could survive the expansion and
to investigate what effects are taking place during the initial release, we
performed simulations of the BEC expansion from harmonic traps using the first
principles STAB method [33, 44]. The simulations started from a cigar-shaped
trap with parameters matched to the experimental conditions. The in-trap state
before release from the trap at time $t=0$ (marked CT in Fig. 3 (a)) was
consistent with the adiabatic sweep theorem applied to the in-situ condensate.
Following expansion from the cigar trap, the simulated tail amplitude
increased and stabilized within a few hundred microseconds (CE in Fig. 3 (a)),
which is much slower than the timescale of the trap potential’s vanishing, and
implies the far-field tails stabilize in appearance much sooner than the 2ms
delay between the trap release and application of the RF and Stern-Gerlach
pulses. Fig. 3 (b) shows the time evolution of the tail amplitude $C_{\rm
sim}$ extracted from a $n(k)=C_{\rm sim}/k^{4}$ fit to the simulated density.
In this configuration the steady-state value of the momentum tails was a
factor of $C_{\rm sim}/\mathcal{C}=$1.64(9) above the predictions of Eqn.
(10). An analysis of the occupation of the tails according to (14), gives very
similar factors $\mathcal{A}_{\rm sim}$ for the increase in the strength of
the tails (relative to _in-situ_ predictions) during evolution, as shown in
Supplementary table S2.

To understand the disagreement with earlier theory [23], which predicted no
depletion survival, we also investigated the effect of adiabatic expansion on
the in-trap depletion. The characteristic healing timescale
$t_{\xi}=\hbar/gn_{0}=15-40\mu$s in the centre of the trapped cloud is
comparable to the trap release time $\tau_{\rm release}$, so a suspicion that
adiabaticity is broken in the CE trap release simulations is warranted. For
example, $t_{\xi}$ is a characteristic timescale for relaxation of density
correlatons due to depletion after a quantum quench [61]. To test the
hypothesis that the difference is due to our system breaking the adiabaticity
assumed in [23], we ran simulations in which the trap is not rapidly released,
but ramped down to half transverse strength over a much longer time period (CS
in Fig. 3). The _in situ_ expression (Eqn 9) predicts that the depletion
should reduce $\propto\bar{\omega}^{6/5}$ to about half its original value.
Indeed it was found that the in-trap contact $C_{\rm sim}$ as well as the the
tail strength $N_{k_{\rm min},k_{\rm max}}$ from (14) decreased roughly as
predicted — see the dashed line in Fig. 3 (b) and Supplementary table S2,
strongly supporting the hypothesis that adiabaticity is needed for agreement
with the results of [23].

As a check on whether we correctly identify the processes involved in
depletion survival, we compared release of atoms from the experimental
elongated clouds with spherically trapped clouds having the same central
density $n_{0}$ and particle number $N$. These clouds are labelled (ST,SE) for
initial and released clouds, respectively. We find that the survival of
depleted atoms is reduced in the spherical trap compared to the elongated
ones.

### Analysis of simulation results

Our understanding of the above dependencies in the simulations is that the
survival and tail strength behaviour are a consequence of the rapid ramp-down
of the trap and quench of density which allows escape of non-condensed
particles, as well as their acceleration by the non-uniform mean-field energy
of the condensate during the expansion.

In detail, after a quench into the untrapped regime, the condensate expands
hydrodynamically on timescales of $1/\omega$, and the equilibrium depletion
density drops in accordance with falling central density $n_{0}$ in Eqn. (10).
However, whether the actual density in k-space modes follows this equilibrium
relationship depends on the reabsorption timescale. Low momentum depletion
atoms are unable to escape the condensate before being reabsorbed and are
absorbed back into the condensate in agreement with [23]. However, if
reabsorption occurs slower than the change in density, the drop in depletion
will be incomplete. High momentum atoms have sufficient velocity to escape the
expanding cloud without being reabsorbed and thus transition to free atoms. In
our system, as seen in Fig. S4 in the supplement, this concerns particles with
wavenumber on the order of $k\gtrsim 2\leavevmode\nobreak\
\mu\textrm{m}^{-1}$, which in particular includes the high momentum tails that
are the focus of the experiment. This is the same kind of escape mechanism
seen for the appearance of halos of $k,-k$ paired atoms in supersonic BEC
collision experiments. [62, 63, 44]

Moreover, an atom inside the BEC experiences an effective force from the
gradient of the mean-field potential $\textbf{F}=-4\pi\hbar^{2}m^{-1}a\nabla
n(x,t)$. This endows escaping depleted particles with a greater momentum. This
phenomenon dubbed a “skiing effect” [64] has been observed for the thermal
part of the cloud in other experiments [65, 66], and for supersonic BEC
collision halos [63, 67, 68]. For a scale-free distribution such as the
$k^{-4}$ power law sought here, such a shift of momentum will manifest itself
as an increase of the amplitude of the tails in the far-field, thus explaining
how the observed depletion can appear stronger than in-situ. The simplest very
rough estimate of this effect can be made by adding an energy of $gn_{0}$ to
each atom during expansion, obtaining a modified density profile of the form
$n(k)\to\approx\mathcal{C}k/(k^{2}-2gn_{0}m/\hbar^{2})^{5/2}$. This leads, for
example, to a doubling of the apparent contact $C_{\rm sim}$ at $k\approx
6/\mu$m for clouds with $n_{0}=39\mu{\rm m}^{-3}$. Thus, this modification
alone is not sufficient to explain the excess counts in the detection region.

A third element is that it is much easier for depletion atoms to escape and
the acceleration is larger along the tightly-confined axes of a cigar-shaped
cloud because the distances $R_{\perp}=(1/\omega_{y,z})\sqrt{2gn_{0}/m}$ are
reduced by $\bar{\omega}/\omega_{y,z}$, whereas the initial mean depletion
velocities in situ $v\sim\sqrt{2gn_{0}/m}$ are isotropic. Indeed, the
simulations show that spherical clouds (SE) exhibit a much weaker effect than
the elongated clouds (CE) in agreement with the longer escape time. This
anisotropy effect also presents as an increase in $C_{\textrm{sim}}$ and
$\mathcal{A}_{\rm sim}$ for simulation collection regions (ROI) that include a
narrower range of angles around the tight trapping plane. Our ability to test
this experimentally was limited because atoms with momenta larger than about 5
$\mu\textrm{m}^{-1}$ in the horizontal plane expanded beyond the detector’s
active surface. Therefore, we obtain only weak evidence of such anisotropy in
the experimental data, which is discussed in the supplementary material.

The above picture is corroborated by another observation within the
simulations: During the expansion we observe a decrease in the total number of
depleted particles (reabsorption) as seen in Supplementary Table S3 by
comparing CE to CT and SE to ST values of $N_{B}$, and a simultaneous increase
of the large-k population (forcing) described by $C_{\rm sim}$. A toy model of
$k,-k$ depletion modes in a uniform gas undergoing external change in the
background density was also investigated to verify our interpretation of the
processes involved.

### Toy model of escape

The reabsorption mechanism and qualitative features of the escape of depletion
from the condensate discussed above can be seen in a toy model of two $k$ and
$-k$ Bogoliubov modes in a uniform volume of gas at zero temperature when the
background density is quenched due to external factors, as described in the
supplementary material. The simplest such “caricature”, when the density
$n_{0}$ is quenched to $n^{\prime}<n_{0}$ at $t=0$ and then stays constant,
has the occupation of each $k$ mode evolve as

$\rho(k,t)=\rho(k,0)-\frac{gn[\epsilon_{0}(k)^{2}-\epsilon(k)^{2}]}{4\epsilon(k)^{2}\epsilon_{0}(k)}\,[1-\cos
2\epsilon(k)t].$ (16)

Here, $\epsilon_{0}(k)$ and $\epsilon(k)$ are given by (3) using the initial
$n_{0}$ and later $n$ values of density, respectively. The reabsorption comes
about then via the initial dip of the Rabi oscillations seen in Fig. 4(a). The
Rabi oscillations are between the two superpositions corresponding to the
initial Bogoliubov ground state and the final one. However, in the actual
expansion it can happen that later (return) stages of the Rabi oscillations
never eventuate if the density drops faster than the oscillation frequency;
roughly $\omega_{\perp}\gtrsim 2\epsilon(k)$.

Figure 4: Toy model simulations for $k$ and $-k$ modes in a uniform gas
initially in the Bogoliubov ground state. Panel (a): behaviour after a
“caricature” quench to half-density, as per Eqn. (16) for modes with
$k\xi=0.02,0.04,\dots,0.18$. The remaining panels concern the better toy model
of Eqs. (S10,S17) with parameters like the full simulation with
$\omega=902\times 895\times 71$ Hz, $N=455852$, $n_{0}=43.66/\mu m^{3}$ peak
density, and $\tau_{\rm release}=38\mu$s, and show the mode occupation
evolution relative to the initial value (the survival rate). Panel (b): For
different initial locations in the condensate in the narrow direction:
$R_{0}=y/R_{\perp}$ where $R_{\perp}=(1/\omega_{y})\sqrt{2gn_{0}/m}$, and
initial $k=1.5/\mu$m. Panel (c): Evolution of the relative occupation for
different ramp speeds $\tau_{\rm release}$. Cycles of reabsorption are seen
for the slow ramps. Panel (d): final survival rate for the same parameters.
Panel (e): shows the dependence of the final survival rate on the trap aspect
ratio $\lambda$, when initial central density $n_{0}$ is kept constant. The
magenta dashed line indicates the experimental $\lambda=12$ (like CE
simulations), the black dashed line a spherical trap $\lambda=1$ (like SE
simulations).

Figs. 4(b-e) show the behaviour of a more careful toy model given by Eqn.
(S10) and (S17) (requiring numerical integration) in which the background
density decays in a time-dependent fashion that at least qualitatively
approximates a significant part of what happens during release. Panel (b)
concerns escape of particles that start in the outer parts of the cloud
($R_{0}=y/R_{\perp}\gtrsim 0.5)$. Panels (c-d) show the dependence of survival
on the speed of the ramp which turns off the trap, indicating that ramp speeds
$\tau_{\rm release}\lesssim 100\mu$s are mostly neutral for the effect but
slower ramps strongly suppress the escape. Panel (e) considers the survival
rate for different trap aspect ratios, including the case of $\lambda=12$ like
in the experiment (CE) and the spherical case (SE). Survival is greatly aided
by elongated traps.

There are still many effects missing from the toy model compared to the STAB
simulations (high momentum atoms on trapped trajectories temporarily outside
the condensate at time of release, multidirectional flight of the depleted
atoms, reduced skiing due to simultaneous collapse of the condensate density,
energy-momentum uncertainty, the effect of the remnant trapping potential on
depletion atoms, to name a few) and survival rate is lower than in the full 3D
simulations. It does, though, give a first qualitative underpinning for the
escape effects seen in the full many-mode simulations.

## Discussion

We find that the number of atoms in the large-$k$ tails in the far-field is
consistent with the tail amplitude scaling as a linear function of the product
$N_{0}n_{0}\propto(N_{0}^{7}\bar{\omega}^{6})^{1/5}$, in line with Tan’s
theory of the contact (Eqns. 10,14). Neither thermal populations nor a number
of technical effects (imaging, single particle processes, background noise)
can be expected to have the same $N_{0}n_{0}$ scaling. However, the effect
size is significantly different than naively expected from in-situ values by a
factor of order 8(3), and from simulated values by a factor of 5(3) which is
not accounted for by any obvious systematic effects.

The earlier experiment by Chang _et al._ [7] also noted the $N_{0}n_{0}$
scaling and an excess by a factor of about 6 that falls within our error bars.
However, recent investigations [35] found that this was correlated with the
presence of a small impurity fraction consisting of $m_{j}=0$ atoms ($\sim
1$%) in their optically trapped cloud, which was otherwise spin polarised in
the $m_{j}=1$ state. When the impurity fraction was reduced to $0.05$% tail
survival was no longer observed. In contrast, our experiment is conducted in a
magnetic trap, which is unable to confine any spin state other than the
$m_{j}=1$ state. Thus our results can not be explained by the presence of
similiar trapped impurities.

While the thermal quasiparticles in the Bogoliubov picture simply map onto the
thermal population of constituent particles of the same momentum (See, for
example, [45] Chap. 8.3. or [2]), a remnant thermal population is not a good
candidate to explain the observations. This is because it decays super-
exponentially with $k$ (Bose-Einstein distribution), and hence does not
account for the atoms we observe beyond $k\gtrsim 6\leavevmode\nobreak\
\mu\textrm{m}^{-1}$, even though it is subject to the same mean-field forcing
as the depletion [65]. We can show this with a simple calculation, noting that
on physical grounds the maximum energy that can be imparted by the “skiing
effect”[64] is $\mu=gn_{0}$, where $n_{0}$ is the initial density in the
centre of the cloud. For an atom with momentum $k=6\mu\textrm{m}^{-1}$, at the
edge of the thermal region in the densest cloud we consider
(44$\mu\textrm{m}^{-3}$), the additional energy $\mu$ imparts at most a
momentum shift of order 0.7$\mu\textrm{m}^{-1}$, which is insufficient to
account for the population detected as far out as $k=10\mu\textrm{m}^{-1}$.
The phonon/particle changeover is also not responsible for the inflections
seen at high $k$ in Figure 1 because this changeover occurs at $k\sim
1/\xi\approx 2\mu{\rm m}^{-1}$.

Another candidate explanation to consider would be large depletion produced
after release in the short-lived mixed-species condensate (where the
$m_{J}=0,1$, and -1 clouds overlap after the Landau-Zener sweep). This could
have the observed $N_{0}n_{0}$ scaling. However, the expression for the
contact in a mixed-species bosonic gas [30] can be combined with the energy of
a condensed mixture [45] to show that the contact in spin-mixed systems is
bounded from above by the contact of the same system polarized in the most
strongly-interacting state [25]. The inter-spin scattering lengths $a_{ij}$
(between He${}^{*}\leavevmode\nobreak\ $atoms in the $i$ and $j$ spin states)
are not fully characterized by experiments, but can be estimated to be
$a_{11}=a_{-1-1}=a_{01}=a_{0-1}\approx 140\leavevmode\nobreak\ a_{0}$,
$a_{00}=120\leavevmode\nobreak\ a_{0}$ and $a_{1-1}\approx
60\leavevmode\nobreak\ a_{0}$, in terms of the Bohr radius $a_{0}$[69]. Thus,
the contact in a He${}^{*}\leavevmode\nobreak\ $condensate is maximized when
the cloud is purely polarized in the $m_{J}=1$ state. Any mixture of He∗ spin
states is thus predicted to have a lower contact (and thus less-populated
tails) than the initial condensate, which appears to rule out this route to
explain our observations.

On the other hand the simulations and toy model demonstrate a route for escape
of the fast depletion atoms from the cloud, and indicate that the survival of
the quantum depletion into the far-field is possible when the release is non
adiabatic, but not as a straightforward mapping into the far-field density.
The latter is due to which imparts some acceleration to the atoms during the
early stages of the expansion [65].

We thus are led to surmise that the experimentally observed tails are indeed a
remnant of the quantum depletion (per the observed scaling with $N_{0}n_{0}$
and $\approx k^{-4}$ that matches the Tan theory, qualitative similarity in
behaviour to the simulations, and lack of convincing counter-hypotheses),
albeit subject to some physical effect during the expansion or some
nonequilibrium enhancement in the trapped state.

In conclusion, we find statistically robust evidence that the quantum
depletion can, remarkably, survive the expansion and dilution of its original
condensate under certain conditions. Our simulations also demonstrate a
mechanism by which the uncondensed quantum depleted atoms of a single species
can be visible in the far-field momentum distribution, and that the
hydrodynamic approximation does not capture sufficient short-wavelength
information to make detailed predictions about the high-momentum behaviour. We
thus find a partial explanation for the experimental deviation of the far-
field distribution from both the in-situ and the hydrodynamic pictures,
although there is an unexplained discrepancy at this time between theory and
experiment as to the amount of this growth. The results reported here expand
the growing body of data and knowledge regarding the somewhat mysterious
behaviour of the far-field quantum depletion [13, 7, 23, 35].

Since the exact mechanism responsible for the impurity effects seen in [35]
remains unrecognised, it is uncertain whether it can also be responsible for
our measurements of tail strength well in excess of the single-species
simulation. Our experiment does not involve impurities in the initial trapped
cloud.

## Methods

### Experimental setup

Our experimental sequence for measurement runs as described in the
“Experimental meaurements” section above is shown schematically in Fig. 5.

We prepared our BECs via forced evaporative cooling in a harmonic magnetic
trap with trap frequencies $\approx(45,425,425)$ Hz and a DC bias stabilized
by our auxiliary field compensation coils [47, 70]. For the tight trap we
increased the coil current after the cooling sequence to obtain trapping
frequencies $\approx(71,902,895)$ Hz, ramping the current as a sigmoid step
function to minimize in-trap oscillations. Note that the weak ($x$) axis of
the trap is horizontal, with tight vertical confinement. The RF pulse was
created by a function generator, amplified, and applied to the experiment
chamber by a coiled antenna inserted into the BiQUIC coil housing. The pulse
swept from 1.6-2.6MHz over 1ms and was centred on the resonance between the
$m_{J}$ states. The determination of the transfer efficiencies $\eta_{J}$ for
each of the $m_{J}$ states is discussed below. The sweep was $10^{6}$-fold
wider than the RF width of the BEC to ensure uniform transfer at all momenta.
Immediately after the RF sweep, the bias coils are switched off and auxiliary
push coils in the vertical (Z) and weak horizontal (X) axes are activated
using a fast MOSFET switch to implement a Stern-Gerlach deflection of the
$m_{J}=-1,$ and $+1$ atoms, such that only $m_{J}=0$ state atoms reach the
detector. The Stern-Gerlach (SG) pulse was designed by increasing the pulse
duration until the $m_{j}=\pm 1$ clouds were given sufficient velocity to
reach the edges of the detector ($\approx 10$ cm/s), and then doubling the
current passed through the field-generating coils.

We use an 80mm diameter multichannel plate and delay-line detector stack [49]
located 848mm below the trap, which registers the arrival times and positions
$(t,x,y)$ of each atom. The velocity of each atom relative to the centre of
mass of each cloud is calculated by
$(v_{x},v_{y},v_{z})=t_{i}^{-1}(x_{i}-\bar{x},y_{i}-\bar{y},\tfrac{1}{2}g_{0}(t_{cen}^{2}-t_{i}^{2}))$,
where $g_{0}$ is the local gravitational acceleration, the overbar denotes the
within-shot average and $t_{cen}$ is the time of flight of the centre of mass
of the cloud. The far-field momentum is thus obtained via
$m\textbf{v}=\hbar\textbf{k}$, noting that this cannot be identified with the
in-situ momentum (see Discussion section). The space and time resolution of
the detector are 100 $\mu$m and 3 $\mu$s, respectively [71]. Sets of ten
experimental runs were interleaved with calibration measurements to determine
the shot-to-shot variation in atom number, trapping frequencies, magnetic
state transfer efficiency, and noise contributions in the manner described in
the supplementary material.

The detector quantum efficiency (QE) of ${\eta_{Q}=}8(2)\%$ was determined
from analysis of the squeezing parameter of correlated atoms on the opposite
sides of scattering halos [72, 73, 74]. A second factor affecting the total
collection efficiency $\epsilon$ is that the $k$-space field of view is
restricted by the detector radius to $k\lesssim 5/\mu{\rm m}$ in the $(x,y)$
plane, which is only just sufficient to reach past the edge of the thermal
region. We thus face a tradeoff in the choice of $k_{\textrm{max}}$, therefore
we define the bounds of our region of interest (ROI) by the minimum elevation
angle $\phi_{c}=\pi/3$ rad above the $(x,y)$ plane and an upper bound of
$k_{\textrm{max}}=10\mu\textrm{m}^{-1}$ (beyond which the signal to noise
ratio becomes too poor). This amounts to an ROI consisting of two vertically
oriented spherical segments, each with half-angle $\pi/6$ from the $z$ axis,
encompassing a total solid angle of
${\Omega_{ROI}=4\pi(1-\sin{\phi_{c}})=}0.13\times 4\pi$ steradians.

We must also account for the state-transfer efficiency of ${\eta_{0}=}25(2)$%
during the RF sweep, and combine all these factors into the total efficiency
$\epsilon{=\eta_{Q}\eta_{0}(1-\sin{\phi_{c}})}\approx 0.23(5)\%$. The dominant
uncertainty in the collection efficiency $\epsilon$ is the 25% error in the
detector quantum efficiency (QE), whereas the other factors (cutoff angle
$\phi_{c}$ and transfer efficiency $\eta_{0}$) are more precisely known.

We performed the analysis of the depletion tails described above for a range
of $\phi_{c}$ and values of the QE $\eta$ and found that the excess of counts
(expressed as $\Lambda_{\textrm{fit}}/\Lambda_{\textrm{pred}}$) was not
significantly affected. This is summarized Supplementary Table S1.

Figure 5: Sketch of the experimental sequence. A BEC is released from a
harmonic trap (a) and expands during freefall before being split into a
superposition of the $m_{J}\in\\{-1,0,1\\}$ states (b) by an RF chirp. A
magnetic field gradient separates the clouds (c) ensuring that only the
magnetically insentitive $m_{J}=0$ cloud lands on the detector (d), from which
the momentum information is reconstructed. Due to the finite detector radius,
the collection region in momentum-space is restricted to two vertically-
oriented spherical segments (shaded region) whose boundary subtends an angle
${\phi_{c}}=\pi/3$ with the horizontal $(x,y)$ plane. The quantum depletion
lies in the dilute tails at large momentum $\gtrsim 6\mu\textrm{m}^{-1}$ (see
Fig. 1)

### Simulations

The STAB method (stochastic Time-Adaptive Bogoliubov) [33, 44] uses the
positive-P representation [75, 76] to describe Bogoliubov quasiparticles
around a dynamically evolving condensate [32]. This allows for straightforward
treatment of inhomogenous and evolving condensates with their associated
quantum depletion, without the need to diagonalise the Bogoliubov-de Gennes
equations. The systems considered here require $4-6\times 10^{6}$ modes for
simulation, so avoiding diagonalisation is very relevant. Previous use of the
STAB method[44, 33, 63, 74, 77, 78, 79, 67] has been according to the
equations described in detail in [33] which relied on a separation of the
condensate and Bogoliubov quasiparticles in k-space that arose from initial
conditions and system dynamics. Here this does not occur, and there is a
significant overlap in momentum space. The standard STAB formulation leads to
an unphysical amplification of the part of the Bogoliubov field that overlaps
with the condensate. Therefore a theory that explicitly imposes orthogonality
between condensate and Bogoliubov modes is required. We summarise our approach
here, with some technical details in the supplemental material. Details of the
derivation and proper benchmarking of the modified method will be reported in
[80].

#### Orthogonalised STAB method

In terms of operators, the Bose field of the atoms
$\widehat{\Psi}({\mathbf{x}},t)$ is written as

$\widehat{\Psi}({\mathbf{x}},t)=\phi({\mathbf{x}},t)+\widehat{\Psi}_{B}({\mathbf{x}},t),$
(17)

where $\phi({\mathbf{x}},t)$ is the condensate order parameter described in
3-dimensional space ${\mathbf{x}}$, and $\widehat{\Psi}_{B}({\mathbf{x}},t)$
is a relatively small operator fluctuation field. The smallness requirement
can be written

$\displaystyle N=$ $\displaystyle\int d^{3}{\mathbf{x}}\
|\phi({\mathbf{x}},t)|^{2}$ (18) $\displaystyle\gg$ $\displaystyle\int
d^{3}{\mathbf{x}}\left\langle\widehat{\Psi}^{\dagger}_{B}({\mathbf{x}},t)\widehat{\Psi}_{B}({\mathbf{x}},t)\right\rangle=N_{B}=N\delta_{B},$
(19)

i.e. $N_{B}$ the number of particles in the Bogoliubov field is small overall,
but locally the Bogoliubov field density need not be smaller then the
condensate – $\delta_{B}$ is the small parameter of the theory [81]. The
condition (19) allows one to discard third and higher orders of
$\widehat{\Psi}_{B}$ in the effective Hamiltonian (the Bogoliubov
approximation). A second condition, not applied in standard STAB, but present
in more precise flavours of Bogoliubov theory is

$\int d^{3}{\mathbf{x}}\
\widehat{\Psi}^{\dagger}_{B}({\mathbf{x}},t)\phi({\mathbf{x}},t)=0.$ (20)

which imposes orthogonality and prevents seeping of condensate atoms into the
fluctuation field $\widehat{\Psi}_{B}({\mathbf{x}},t)$.

The condensate order parameter $\phi({\mathbf{x}},t)$ is assumed to evolve
according to the Gross-Pitaevskii equation (correct to leading order, given
(19)):

$\displaystyle
i\hbar\frac{d\phi}{dt}=\left[-\frac{\hbar^{2}}{2m}\nabla^{2}+g|\phi|^{2}+V({\mathbf{x}},t)\right]\phi.$
(21)

and is normalised to the (conserved) total number of particles $\int
d^{3}{\mathbf{x}}\ |\phi({\mathbf{x}},t)|^{3}=N$. The
$g=4\pi\hbar^{2}a_{1,1}/m$ is the s-wave contact interaction between He∗ atoms
in the initial $m_{J}=1$ state (we take $a_{1,1}=7.51$nm), and
$V({\mathbf{x}},t)$ is the trap potential with in general time-dependent
frequency. We then represent the Bogoliubov quasiparticles using the
positive-P representation [75, 33], which leads to the following equations of
motion:

$\displaystyle i\hbar\frac{d\psi_{B}}{dt}$ $\displaystyle=$
$\displaystyle\left[-\frac{\hbar^{2}}{2m}\nabla^{2}+g|\phi|^{2}+V({\mathbf{x}},t)\right]\psi_{B}+{\mathcal{P}}_{\perp}\left\\{g|\phi|^{2}\psi_{B}+g\phi^{2}{\widetilde{\psi}}_{B}^{*}+\sqrt{-ig}\
\phi\,\xi({\mathbf{x}},t)\right\\}$ (22a) $\displaystyle
i\hbar\frac{d{\widetilde{\psi}}_{B}}{dt}$ $\displaystyle=$
$\displaystyle\left[-\frac{\hbar^{2}}{2m}\nabla^{2}+g|\phi|^{2}+V({\mathbf{x}},t)\right]{\widetilde{\psi}}_{B}+{\mathcal{P}}_{\perp}\left\\{g|\phi|^{2}{\widetilde{\psi}}_{B}+g\phi^{2}\psi_{B}^{*}+\sqrt{-ig}\
\phi\,{\widetilde{\xi}}({\mathbf{x}},t)\right\\}.$ (22b)

Here the ket $\psi_{B}({\mathbf{x}},t)$ and bra
${\widetilde{\psi}}_{B}({\mathbf{x}},t)$ amplitudes provide the positive-P
representation of the Bogoliubov field $\hat{\Psi}_{B}({\mathbf{x}},t)$ in 3D
space. We used the robust stochastic integration procedure described in [82].
The $\xi({\mathbf{x}},t)$ and ${\widetilde{\xi}}({\mathbf{x}},t)$ are
independent white Gaussian noise fields of zero mean and variance:

$\langle\xi({\mathbf{x}},t)\xi({\mathbf{x}}^{\prime},t^{\prime})\rangle=\langle{\widetilde{\xi}}({\mathbf{x}},t){\widetilde{\xi}}({\mathbf{x}}^{\prime},t^{\prime})\rangle=\delta^{3}({\mathbf{x}}-{\mathbf{x}}^{\prime})\delta(t-t^{\prime}).$
(23)

An ensemble of field trajectories with independent noise in each trajectory
and in each trajectory’s initial state is generated to represent the
Bogoliubov field. We typically used ${\mathcal{S}}=4000$ trajectories.
Notably, the equations (22) allow not only for production of additional
Bogoliubov quasiparticles quantum depleted from the condensate but also for
their reabsorption. The main additional element in (22) compared to the
standard STAB equations [79] is the projection ${\mathcal{P}}_{\perp}$ which
imposes the orthogonality requirement (20) and avoids the aforementioned
amplification of the Bogoliubov field where it overlaps with the condensate.
The projection ${\mathcal{P}}_{\perp}$ of a field $f({\mathbf{x}})$ can be
carried out efficiently by

${\mathcal{P}}_{\perp}f({\mathbf{x}})=f({\mathbf{x}})-\frac{1}{N}\left[\int
d^{3}{{\mathbf{x}}^{\prime}}\
\phi({{\mathbf{x}}^{\prime}})^{*}f({{\mathbf{x}}^{\prime}})\right]\,\phi({\mathbf{x}}).$
(24)

The kinetic part of the evolution equations (21)-(22) is also carried out
efficiently by a split-step approach which evaluates kinetic terms in k-space
and the rest in x-space, moving between k-space and x-space using a fast
Fourier transform. Calculation of observables is described in the
supplementary material.

#### Initial condition

Our simulations aim to study the evolution of the quantum depletion particles
in $\widehat{\Psi}_{B}$ after release from the trap. We use a zero temperature
initial condition, since the object is to study the behaviour of the high
momentum tails beyond the edge of the thermal cloud, in which $T>0$ effects
are negligible. The $T=0$ initial state is also more straightforward to
obtain, allowing one to use lower k-values to access the $k^{-4}$ tails, since
they are not obscured by the stronger thermal cloud at intermediate momenta.
This significantly reduces the size of the computational lattice needed. For
the low temperatures in the experiment we do not expect any significant
interaction between the behavior of the thermal cloud and the depleted atoms
because both are well approximated by the Bogoliubov approximation which
neglects interactions between excited modes. Therefore the neglect of the
thermal cloud does not significantly affect the properties of the higher k
depletion or its evolution.

However, one cannot use the standard Gross-Pitaevskii ground state since that
has 100% condensate and no quantum depletion. The task of generating a cloud
with the appropriate depletion in such a large nonuniform system turns out to
be nontrivial. Conceptually the issue is simple – diagonalise the Bogoliubov
Hamiltonian, and give the well known Bogoliubov $T=0$ occupation to each
quasiparticle mode. However, for a system with $10^{6}$ modes diagonalisation
is not a good option. Our solution to this situation is to make a calibrated
quantum quench from the Gross-Pitaevskii solution to the full Bogoliubov
equations of motion, which provides a state with an appropriate quantity of
quantum depletion. The technique is described in detail in the supplementary
material.

#### Simulation types

Several types of simulations were made, with shorthand labels as per Fig. 3,
and summarised in Supplementary Table S3:

(CE) Release of atoms from the trap, as in the experiment. Here the potential
was reduced exponentially

$V({\mathbf{x}})=\frac{m}{2}\left(\omega_{x}^{2}x^{2}+\omega_{y}^{2}y^{2}+\omega_{z}^{2}z^{2}\right)\,e^{-t/\tau_{\rm
release}},$ (25)

with time constant $\tau_{\rm release}=37.5\mu$s, matched to the experiment.
The initial trap frequencies were $\omega=425\times 425\times 45$ Hz and
$\omega=902\times 895\times 71$ Hz, and two variants of the initial state were
simulated: a low density and a high density cloud.

(CS) Slow decrease of the transverse trapping frequencies by a factor of two.
Here we ramped the trap as follows:

$V({\mathbf{x}})=\frac{m}{2}\left[\left(\omega_{x}^{2}x^{2}+\omega_{y}^{2}y^{2}\right)\left(1-\frac{t}{2t_{\rm
ramp}}\right)^{2}+\omega_{z}^{2}z^{2}\right],$ (26)

with timescales of order 1-2ms (see Supplementary Table S3). The simulations
were run up till $t=t_{\rm ramp}$ when the transverse trap frequency was half
the original one.

(ST,SE) Trap release of spherical clouds. The velocity distribution of the
depletion atoms is isotropic in situ, being given by $mv^{2}/2\approx gn_{0}$.
However, the distance to travel to escape reabsorption depends on the cloud
shape. In particular, escape is made easier in the tight trap directions (less
distance to travel), and harder in the long trap direction. Here we used
spherically trapped clouds having the same central density $n_{0}$ and
particle number $N$. These clouds had isotropic trapping frequency
${\overline{\omega}}=(\omega_{x}\omega_{y}\omega_{z})^{1/3}$ and are labelled
(ST). Trap release (SE) followed (25) as before.

## Acknowledgements

We Would like to thank David Clement, Jean Dalibard, Karen Kherunstyan, and
Raphael Lopes for their helpful discussions. This work was supported by
Australian Research Council (ARC) Discovery Project Grants No. DP160102337 and
No. DP190103021. S. S. H was supported by DECRA DE150100315, J.A.R., D. K. S.
by the Australian Postgraduate Award (APA), and K.F.T. by the Australian
Government Research Training Program (RTP) Scholarship. The simulations by P.
D. were supported by National Science Centre (Poland) grants No.
2018/31/B/ST2/01871 and 2012/07/E/ST2/01389.

## Author contributions statement

The data was collected and analysed by J. A. Ross, D. K. Shin, K. F. Thomas,
and B. M. Henson under the supervision of S. S. Hodgman and A. G. Truscott. P.
Deuar conducted the simulations. All authors contributed to the interpretation
of results. J. A. Ross and P. Deuar wrote the paper with input from all
authors.

## Additional information

The authors declare no competing interests. The data generated and analysed
during the current study are available from the corresponding author on
reasonable request.

## References

  * [1] Bogoliubov, N.  On the theory of superfluidity.  _Journal of Physics USSR_ XI, 23 (1947).
  * [2] Vogels, J. M., Xu, K., Raman, C., Abo-Shaeer, J. R. & Ketterle, W.  Experimental observation of the Bogoliubov transformation for a Bose-Einstein condensed gas.  _Phys. Rev. Lett._ 88, 060402, DOI: 10.1103/PhysRevLett.88.060402 (2002).
  * [3] Pitaevskiǐ, L. P. & Stringari, S.  _Bose-Einstein condensation and superfluidity_.  No. 164 in International series of monographs on physics (Oxford University Press, Oxford, United Kingdom, 2016), first edition edn.  OCLC: ocn919186901.
  * [4] Decamp, J., Albert, M. & Vignolo, P.  Tan’s contact in a cigar-shaped dilute Bose gas.  _Phys. Rev. A_ 97, 033611, DOI: 10.1103/PhysRevA.97.033611 (2018).
  * [5] Stewart, J. T., Gaebler, J. P., Drake, T. E. & Jin, D. S.  Verification of universal relations in a strongly interacting Fermi gas.  _Phys. Rev. Lett._ 104, 235301, DOI: 10.1103/PhysRevLett.104.235301 (2010).
  * [6] Wild, R. J., Makotyn, P., Pino, J. M., Cornell, E. A. & Jin, D. S.  Measurements of Tan’s contact in an atomic Bose-Einstein condensate.  _Phys. Rev. Lett._ 108, 145305, DOI: 10.1103/PhysRevLett.108.145305 (2012).
  * [7] Chang, R. _et al._ Momentum-resolved observation of thermal and quantum depletion in a Bose gas.  _Phys. Rev. Lett._ 117, 235303, DOI: 10.1103/PhysRevLett.117.235303 (2016).
  * [8] Makotyn, P., Klauss, C. E., Goldberger, D. L., Cornell, E. A. & Jin, D. S.  Universal dynamics of a degenerate unitary Bose gas.  _Nature Physics_ 10, 116–119, DOI: 10.1038/nphys2850 (2014).
  * [9] Eigen, C. _et al._ Universal prethermal dynamics of Bose gases quenched to unitarity.  _Nature_ 221, DOI: https://doi.org/10.1038/s41586-018-0674-1 (2018).
  * [10] Xu, K. _et al._ Observation of strong quantum depletion in a gaseous Bose-Einstein condensate.  _Phys. Rev. Lett._ 96, 180405, DOI: 10.1103/PhysRevLett.96.180405 (2006).
  * [11] Pieczarka, M. _et al._ Observation of quantum depletion in a non-equilibrium exciton–polariton condensate.  _Nature Communications_ 11, 429, DOI: 10.1038/s41467-019-14243-6 (2020).
  * [12] Lopes, R. _et al._ Quantum depletion of a homogeneous Bose-Einstein condensate.  _Phys. Rev. Lett._ 119, 190404, DOI: 10.1103/PhysRevLett.119.190404 (2017).
  * [13] Cayla, H. _et al._ Hanbury Brown and Twiss bunching of phonons and of the quantum depletion in an interacting Bose gas.  _Phys. Rev. Lett._ 125, 165301, DOI: 10.1103/PhysRevLett.125.165301 (2020).
  * [14] Kuhnle, E. D. _et al._ Temperature dependence of the universal contact parameter in a unitary Fermi gas.  _Phys. Rev. Lett._ 106, 170402, DOI: 10.1103/PhysRevLett.106.170402 (2011).
  * [15] Sagi, Y., Drake, T. E., Paudel, R. & Jin, D. S.  Measurement of the homogeneous contact of a unitary Fermi gas.  _Phys. Rev. Lett._ 109, 220402, DOI: 10.1103/PhysRevLett.109.220402 (2012).
  * [16] Fletcher, R. J. _et al._ Two- and three-body contacts in the unitary Bose gas.  _Science_ 355, 377–380, DOI: 10.1126/science.aai8195 (2017).
  * [17] Lopes, R. _et al._ Quasiparticle energy in a strongly interacting homogeneous Bose-Einstein condensate.  _Phys. Rev. Lett._ 118, 210401, DOI: 10.1103/PhysRevLett.118.210401 (2017).
  * [18] Mukherjee, B. _et al._ Spectral response and contact of the unitary Fermi gas.  _Phys. Rev. Lett._ 122, 203402, DOI: 10.1103/PhysRevLett.122.203402 (2019).
  * [19] Carcy, C. _et al._ Contact and sum rules in a near-uniform Fermi gas at unitarity.  _Phys. Rev. Lett._ 122, 203401, DOI: 10.1103/PhysRevLett.122.203401 (2019).
  * [20] Colussi, V. E. _et al._ Cumulant theory of the unitary Bose gas: Prethermal and Efimovian dynamics.  _Phys. Rev. A_ 102, 063314, DOI: 10.1103/PhysRevA.102.063314 (2020).
  * [21] Kira, M.  Coherent quantum depletion of an interacting atom condensate.  _Nature Communications_ 6, 6624, DOI: 10.1038/ncomms7624 (2015).
  * [22] Smith, D. H., Braaten, E., Kang, D. & Platter, L.  Two-body and three-body contacts for identical bosons near unitarity.  _Phys. Rev. Lett._ 112, 110402, DOI: 10.1103/PhysRevLett.112.110402 (2014).
  * [23] Qu, C., Pitaevskii, L. P. & Stringari, S.  Expansion of harmonically trapped interacting particles and time dependence of the contact.  _Phys. Rev. A_ 94, 063635, DOI: 10.1103/PhysRevA.94.063635 (2016).
  * [24] Braaten, E., Kang, D. & Platter, L.  Short-time operator product expansion for rf spectroscopy of a strongly interacting Fermi gas.  _Phys. Rev. Lett._ 104, 223004, DOI: 10.1103/PhysRevLett.104.223004 (2010).
  * [25] Braaten, E., Kang, D. & Platter, L.  Universal relations for identical bosons from three-body physics.  _Phys. Rev. Lett._ 106, 153005, DOI: 10.1103/PhysRevLett.106.153005 (2011).
  * [26] Rakhimov, A.  Tan’s contact as an indicator of completeness and self-consistency of a theory.  _Phys. Rev. A_ 102, 063306, DOI: 10.1103/PhysRevA.102.063306 (2020).
  * [27] Braaten, E. & Platter, L.  Exact relations for a strongly interacting Fermi gas from the operator product expansion.  _Phys. Rev. Lett._ 100, 205301, DOI: 10.1103/PhysRevLett.100.205301 (2008).
  * [28] Zhang, S. & Leggett, A. J.  Universal properties of the ultracold Fermi gas.  _Phys. Rev. A_ 79, 023601, DOI: 10.1103/physreva.79.023601 (2009).
  * [29] Combescot, R., Alzetto, F. & Leyronas, X.  Particle distribution tail and related energy formula.  _Phys. Rev. A_ 79, 053640, DOI: 10.1103/PhysRevA.79.053640 (2009).
  * [30] Werner, F. & Castin, Y.  General relations for quantum gases in two and three dimensions. II. Bosons and mixtures.  _Phys. Rev. A_ 86, 053633, DOI: 10.1103/PhysRevA.86.053633 (2012).
  * [31] Werner, F. & Castin, Y.  General relations for quantum gases in two and three dimensions: Two-component fermions.  _Phys. Rev. A_ 86, 013626, DOI: 10.1103/PhysRevA.86.013626 (2012).
  * [32] Sinatra, A., Castin, Y. & Lobo, C.  A Monte Carlo formulation of the Bogolubov theory.  _Journal of Modern Optics_ 47, 2629–2644, DOI: 10.1080/09500340008232186 (2000).
  * [33] Deuar, P., Chwedeńczuk, J., Trippenbach, M. & Ziń, P.  Bogoliubov dynamics of condensate collisions using the positive-P representation.  _Phys. Rev. A_ 83, 063625, DOI: 10.1103/PhysRevA.83.063625 (2011).
  * [34] Steinhauer, J. _et al._ Bragg spectroscopy of the multibranch Bogoliubov spectrum of elongated Bose-Einstein condensates.  _Phys. Rev. Lett._ 90, 060404, DOI: 10.1103/PhysRevLett.90.060404 (2003).
  * [35] Cayla, H. _et al._ Observation of $1/k^{4}$-tails in the asymptotic momentum distribution of bose polarons.  _arXiv: 2204.10697_ DOI: 10.48550/arXiv.2204.10697 (2022).
  * [36] Dmowski, W. _et al._ Observation of dynamic atom-atom correlation in liquid helium in real space.  _Nature Communications_ 8, 15294, DOI: doi:10.1038/ncomms15294 (2017).
  * [37] Glyde, H. R., Azuah, R. T. & Stirling, W. G.  Condensate, momentum distribution, and final-state effects in liquid ${}^{4}\mathrm{He}$.  _Phys. Rev. B_ 62, 14337–14349, DOI: 10.1103/PhysRevB.62.14337 (2000).
  * [38] Moroni, S. & Boninsegni, M.  Condensate fraction in liquid 4He.  _Journal of Low Temperature Physics_ 136, 129, DOI: https://doi.org/10.1023/B:JOLT.0000038518.10132.30 (2004).
  * [39] Tan, S.  Large momentum part of a strongly correlated Fermi gas.  _Annals of Physics_ 323, 2971–2986, DOI: 10.1016/j.aop.2008.03.005 (2008).
  * [40] Tan, S.  Energetics of a strongly correlated Fermi gas.  _Annals of Physics_ 323, 2952–2970, DOI: 10.1016/j.aop.2008.03.004 (2008).
  * [41] Tan, S.  Generalized virial theorem and pressure relation for a strongly correlated Fermi gas.  _Annals of Physics_ 323, 2987–2990, DOI: 10.1016/j.aop.2008.03.003 (2008).
  * [42] Hoinka, S. _et al._ Precise determination of the structure factor and contact in a unitary Fermi gas.  _Phys. Rev. Lett._ 110, 055305, DOI: 10.1103/PhysRevLett.110.055305 (2013).
  * [43] Tenart, A., Hercé, G., Bureik, J., Dareau, A. & Clément, D.  Observation of pairs of atoms at opposite momenta in an equilibrium interacting bose gas.  _Nature Physics_ 17, DOI: https://doi.org/10.1038/s41567-021-01381-2 (2021).
  * [44] Kheruntsyan, K. V. _et al._ Violation of the Cauchy-Schwarz inequality with matter waves.  _Phys. Rev. Lett._ 108, 260401, DOI: 10.1103/PhysRevLett.108.260401 (2012).
  * [45] Pethick, C. & Smith, H.  _Bose-Einstein condensation in dilute gases_ (Cambridge University Press, Cambridge, New York, 2008), 2nd ed edn.
  * [46] Olshanii, M. & Dunjko, V.  Short-distance correlation properties of the Lieb-Liniger system and momentum distributions of trapped one-dimensional atomic gases.  _Phys. Rev. Lett._ 91, 090401, DOI: 10.1103/PhysRevLett.91.090401 (2003).
  * [47] Dall, R. & Truscott, A.  Bose–Einstein condensation of metastable helium in a bi-planar quadrupole Ioffe configuration trap.  _Optics Communications_ 270, 255–261, DOI: 10.1016/j.optcom.2006.09.031 (2007).
  * [48] Moal, S. _et al._ Accurate determination of the scattering length of metastable helium atoms using dark resonances between atoms and exotic molecules.  _Phys. Rev. Lett._ 96, DOI: 10.1103/physrevlett.96.023203 (2006).
  * [49] Manning, A. G., Hodgman, S. S., Dall, R. G., Johnsson, M. T. & Truscott, A. G.  The Hanbury Brown-Twiss effect in a pulsed atom laser.  _Optics Express_ 18, 18712, DOI: 10.1364/oe.18.018712 (2010).
  * [50] Hodgman, S. S. _et al._ Metastable helium: A new determination of the longest atomic excited-state lifetime.  _Phys. Rev. Lett._ 103, 053002, DOI: 10.1103/physrevlett.103.053002 (2009).
  * [51] Hodgman, S. S., Dall, R. G., Manning, A. G., Baldwin, K. G. H. & Truscott, A. G.  Direct measurement of long-range third-order coherence in Bose-Einstein condensates.  _Science_ 1046–1049 (2011).
  * [52] Dall, R. G. _et al._ Ideal _n_ -body correlations with massive particles.  _Nature physics_ 341–344 (2013).
  * [53] Schellekens, M. _et al._ Hanbury Brown Twiss effect for ultracold quantum gases.  _Science_ 648–651 (2005).
  * [54] Jeltes, T. _et al._ Comparison of the Hanbury Brown–Twiss effect for bosons and fermions.  _Nature_ 402–405 (2007).
  * [55] Dall, R. G. _et al._ Observation of atomic speckle and Hanbury Brown–Twiss correlations in guided matter waves.  _Nature communications_ (2011).
  * [56] Perrin, A. _et al._ Observation of atom pairs in spontaneous four-wave mixing of two colliding Bose-Einstein condensates.  _Phys. Rev. Lett._ 99, 150405, DOI: 10.1103/PhysRevLett.99.150405 (2007).
  * [57] A, P. _et al._ Hanbury Brown and Twiss correlations across the Bose-Einstein condensation threshold.  _Nature physics_ 195–198 (2012).
  * [58] Dalfovo, F., Giorgini, S., Pitaevskii, L. P. & Stringari, S.  Theory of bose-einstein condensation in trapped gases.  _Rev. Mod. Phys._ 71, 463–512, DOI: 10.1103/RevModPhys.71.463 (1999).
  * [59] Clauset, A., Shalizi, C. R. & Newman, M. E. J.  Power-law distributions in empirical data.  _SIAM Review_ 51, 661–703, DOI: 10.1137/070710111 (2009).  ArXiv: 0706.1062.
  * [60] Virkar, Y. & Clauset, A.  Power-law distributions in binned empirical data.  _The Annals of Applied Statistics_ 8, 89–119, DOI: 10.1214/13-AOAS710 (2014).  ArXiv: 1208.3524.
  * [61] Pietraszewicz, J., Stobińska, M. & Deuar, P.  Correlation evolution in dilute Bose-Einstein condensates after quantum quenches.  _Phys. Rev. A_ 99, 023620, DOI: 10.1103/PhysRevA.99.023620 (2019).
  * [62] Vogels, J. M., Xu, K. & Ketterle, W.  Generation of macroscopic pair-correlated atomic beams by four-wave mixing in bose-einstein condensates.  _Phys. Rev. Lett._ 89, 020401, DOI: 10.1103/PhysRevLett.89.020401 (2002).
  * [63] Krachmalnicoff, V. _et al._ Spontaneous four-wave mixing of de Broglie waves: Beyond optics.  _Phys. Rev. Lett._ 104, 150402, DOI: 10.1103/PhysRevLett.104.150402 (2010).
  * [64] Deuar, P., Ziń, P., Chwedeńczuk, J. & Trippenbach, M.  Mean field effects on the scattered atoms in condensate collisions.  _The European Physical Journal D_ 65, 19–24, DOI: 10.1140/epjd/e2011-20066-7 (2011).
  * [65] Ozeri, R., Steinhauer, J., Katz, N. & Davidson, N.  Direct observation of the phonon energy in a bose-einstein condensate by tomographic imaging.  _Phys. Rev. Lett._ 88, 220401, DOI: 10.1103/PhysRevLett.88.220401 (2002).
  * [66] Simsarian, J. E. _et al._ Imaging the phase of an evolving bose-einstein condensate wave function.  _Phys. Rev. Lett._ 85, 2040–2043, DOI: 10.1103/PhysRevLett.85.2040 (2000).
  * [67] Deuar, P. _et al._ Anisotropy in s -wave Bose-Einstein condensate collisions and its relationship to superradiance.  _Phys. Rev. A_ 90, 033613 (2014).
  * [68] Hodgman, S. S., Khakimov, R. I., Lewis-Swan, R. J., Truscott, A. G. & Kheruntsyan, K. V.  Solving the quantum many-body problem via correlations measured with a momentum microscope.  _Phys. Rev. Lett._ 118, 240402, DOI: 10.1103/PhysRevLett.118.240402 (2017).
  * [69] Vassen, W., Notermans, R. P. M. J. W., Rengelink, R. J. & van der Beek, R. F. H. J.  Ultracold metastable helium: Ramsey fringes and atom interferometry.  _Applied Physics B_ 122, 289, DOI: 10.1007/s00340-016-6563-0 (2016).
  * [70] Dedman, C. J., Dall, R. G., Byron, L. J. & Truscott, A. G.  Active cancellation of stray magnetic fields in a Bose-Einstein condensation experiment.  _Review of Scientific Instruments_ 78, 024703, DOI: 10.1063/1.2472600 (2007).
  * [71] Henson, B. M. _et al._ Bogoliubov-Cherenkov radiation in an atom laser.  _Phys. Rev. A_ 97, 063601, DOI: 10.1103/physreva.97.063601 (2018).
  * [72] Shin, D. K. _et al._ Bell correlations between spatially separated pairs of atoms.  _Nature communications_ 4447, DOI: https://doi.org/10.1038/s41467-019-12192-8 (2019).
  * [73] Shin, D. K., Ross, J. A., Henson, B. M., Hodgman, S. S. & Truscott, A. G.  Entanglement-base 3D magnetic gradiometry with an ultracold atomic scattering halo.  _New Journal of Physics_ 013002, DOI: 10.1088/1367-2630/ab66de (2020).
  * [74] Jaskula, J.-C. _et al._ Sub-Poissonian number differences in four-wave mixing of matter waves.  _Phys. Rev. Lett._ 105, 190402, DOI: 10.1103/PhysRevLett.105.190402 (2010).
  * [75] Drummond, P. D. & Gardiner, C. W.  Generalised P-representations in quantum optics.  _Journal of Physics A: Mathematical and General_ 13, 2353 (1980).
  * [76] Deuar, P. & Drummond, P. D.  Correlations in a BEC collision: First-principles quantum dynamics with 150 000 atoms.  _Phys. Rev. Lett._ 98, 120402, DOI: 10.1103/PhysRevLett.98.120402 (2007).
  * [77] Lewis-Swan, R. J. & Kheruntsyan, K. V.  Proposal for demonstrating the Hong–Ou–Mandel effect with matter waves.  _Nature Commun._ 5, 3752, DOI: 10.1038/ncomms4752 (2014).
  * [78] Lewis-Swan, R. J. & Kheruntsyan, K. V.  Proposal for a motional-state Bell inequality test with ultracold atoms.  _Phys. Rev. A_ 91, 052114, DOI: 10.1103/PhysRevA.91.052114 (2015).
  * [79] Deuar, P., Wasak, T., Ziń, P., Chwedeńczuk, J. & Trippenbach, M.  Tradeoffs for number squeezing in collisions of Bose-Einstein condensates.  _Phys. Rev. A_ 88, 013617, DOI: 10.1103/PhysRevA.88.013617 (2013).
  * [80] Deuar, P., Ross, J. D. & Truscott, A. G. (2021).  In preparation.
  * [81] Castin, Y. & Dum, R.  Low-temperature Bose-Einstein condensates in time-dependent traps: Beyond the u(1) symmetry-breaking approach.  _Phys. Rev. A_ 57, 3008–3021, DOI: https://doi.org/10.1103/PhysRevA.57.3008 (1998).
  * [82] Deuar, P.  Multi-time correlations in the positive-P, Q, and doubled phase-space representations.  _Quantum_ 5, 455, DOI: https://doi.org/10.22331/q-2021-05-10-455 (2021).
  * [83] Henson, B. M. _et al._ Approaching the adiabatic timescale with machine learning.  _Proceedings of the National Academy of Sciences_ 115, 13216, DOI: https://doi.org/10.1073/pnas.1811501115 (2018).
  * [84] Goldstein, M. L., Morris, S. A. & Yen, G. G.  Problems with fitting to the power-law distribution.  _The European Physical Journal B_ 41, 255–258, DOI: 10.1140/epjb/e2004-00316-5 (2004).
  * [85] Hanel, R., Corominas-Murtra, B., Liu, B. & Thurner, S.  Fitting power-laws in empirical data with estimators that work for all exponents.  _PLOS ONE_ 12, e0170920, DOI: 10.1371/journal.pone.0170920 (2017).
  * [86] Ruostekoski, J. & Martin, A. D.  _The Truncated Wigner Method for Bose Gases_ , chap. 13, 203–214 (Imperial College Press, 2013).
  * [87] Sinatra, A., Castin, Y., Carusotto, I., Lobo, C. & Witkowska, E.  _Number-Conserving Stochastic Approaches for Equilibrium and Time-Dependent Bose Gases_ , chap. 14, 215–228 (Imperial College Press, 2013).
  * [88] Martin, A. D. & Ruostekoski, J.  Nonequilibrium quantum dynamics of atomic dark solitons.  _New Journal of Physics_ 12, 055018 (2010).
  * [89] Martin, A. D. & Ruostekoski, J.  Quantum and thermal effects of dark solitons in a one-dimensional Bose gas.  _Phys. Rev. Lett._ 104, 194102, DOI: 10.1103/PhysRevLett.104.194102 (2010).
  * [90] Sinatra, A., Lobo, C. & Castin, Y.  The truncated Wigner method for Bose-condensed gases: limits of validity and applications.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ 35, 3599 (2002).
  * [91] Norrie, A. A., Ballagh, R. J. & Gardiner, C. W.  Quantum turbulence and correlations in Bose-Einstein condensate collisions.  _Phys. Rev. A_ 73, 043617, DOI: 10.1103/PhysRevA.73.043617 (2006).
  * [92] Deuar, P.  _First-principles quantum simulations of many-mode open interacting Bose gases using stochastic gauge methods_.  Ph.D. thesis, University of Queensland, arXiv:cond-mat/0507023 (2005).
  * [93] Drummond, P. D. & Opanchuk, B.  Initial states for quantum field simulations in phase space.  _Phys. Rev. Research_ 2, 033304, DOI: 10.1103/PhysRevResearch.2.033304 (2020).
  * [94] Castin, Y. & Dum, R.  Bose-Einstein condensates in time dependent traps.  _Phys. Rev. Lett._ 77, 5315–5319, DOI: https://doi.org/10.1103/PhysRevLett.77.5315 (1996).

Supplementary information to
_On the survival of the quantum depletion of a condensate after release from a
magnetic trap_

J. A. Ross, P. Deuar, D. K. Shin, K. F. Thomas, B. M. Henson, S. S. Hodgman,
and A. G. Truscott*
* Corresponding author email<EMAIL_ADDRESS>
## Experimental details

#### Peak density calibration

The quantum depletion and contact are both predicted to depend solely on the
condensed number and trapping frequencies via the condensate density, hence it
is important to determine both quantities accurately.

The sole experimental parameters in the expression for the peak density (Eqn.
(11) ) are the geometric trap frequency
$\bar{\omega}=\left(\omega_{x}\cdot\omega_{y}\cdot\omega_{z}\right)^{1/3}$,
and $N_{0}$, the number of atoms in the condensate. We simultaneously
determine the total atom number $N$ and trap frequency $\bar{\omega}$ in a
single shot using a pulsed atom laser and use the thermal fraction $\eta_{T}$
(see below) to determine the condensed number $N_{0}=(1-\eta_{T})N$.

The pulsed atom laser consists of a series of Fourier-broadened RF pulses
centred on the minimum Zeeman splitting in the trap. The pulse transfers atoms
in the trap to the untrapped $m_{J}=0$ state with an approximately constant
transfer rate across the cloud [49, 71]. We outcouple approximately 2% of the
atoms per 100$\leavevmode\nobreak\ \mu$s pulse for $\approx$200 pulses, which
eventually depletes the entire trap. The atom laser thus prevents the detector
from saturating and allows an accurate determination of the atom number, up to
a factor of the quantum efficiency. We determine the trapping frequencies by
inducing centre-of-mass oscillations with a magnetic impulse, and find the
oscillation period from the atom laser pulses [83].

#### Determining spin transfer efficiency

To calibrate the transfer efficiencies, we applied a weaker Stern-Gerlach than
for the depletion measurement, such that each $m_{J}$ cloud hits detector at
different times (and positions), as illustrated in Figure S1. The efficiencies
$\eta_{J}$ cannot be calculated by counting the atoms in each cloud because
the detector saturates during the peak condensate flux. However, we can
compare the thermal parts. We align each cloud along the time (Z) axis and
compute the pointwise fraction of the atomic flux $\phi(t)$ accounted for by
each cloud, $\eta_{j}(t)=\phi_{j}(t)/\sum_{j}\phi_{j}(t)$, as depicted in
Figure S1 (a-c). The ratio of densities between the clouds is roughly constant
in the thermal part (Figure S1 (c)), indicating the absence of saturation
effects in the thermal part and a spin transfer efficiency that is independent
of $k$. The fraction of the original cloud transferred into each $m_{J}$ state
is determined by taking the average $\langle\eta_{j}(t)\rangle$ over the
thermal tails.

Figure S1: Determining the RF transfer efficiency. The time-of-flight profiles
of each pulse are resolved (a) by applying a weak Stern-Gerlach pulse during
the time of flight. The pulses are aligned with respect to their centre-of-
mass (b) and used to determine the pointwise fraction ((c), dotted line).
Detector saturation is evident in the peaks (dashed lines), but not in the
thermal tails (solid lines), which are used to compute the transfer
efficiency. Because of its lower flux, the $m_{J}=-1$ pulse does not show any
clear evidence of saturation (d) and is used to determine the thermal fraction
and hence $N_{0}$.

We find these efficiencies are approximately 74%, 24%, and 2% in all runs for
the $m_{J}=+1$, 0, and -1 states, respectively. While the $m_{J}=0$ and
$m_{J}=1$ clouds clearly saturate the detector, the small fraction ($\approx
2\%$) of the atoms transferred to the $m_{J}=-1$ state does not (Figure S1
(d)). A bimodal fit to the condensed and thermal parts, plus constant
background, yields the thermal fraction $\eta_{T}$ and condensed fraction
$1-\eta_{T}$.

#### Noise sources

In early tests of our measurement sequence we noticed a contamination of the
signal by spurious counts. Specifically, an extremely dilute remnant of the
$m_{J}=+1$ cloud was observed as a small peak in counts in the region
$-10\mu\textrm{m}^{-1}\lesssim k_{z}\lesssim-8.5\mu\textrm{m}^{-1}$. We
inferred these were remnant counts from the $m_{J}=+1$ cloud as they were
still visible when we ran an experimental sequence without the Landau-Zener
transfer, but including the Stern-Gerlach pulse.

While the cause of the appearance of the remnant counts is unclear, we observe
that the count density outside the region of interest is similar in both the
shots with the RF pulse (Landau-Zener transfer) and those without. We note
that only about one in a million atoms from the $m_{J}=1$ cloud are present in
this manner in a given shot. On average, this accounted for approximately
$N_{mJ=+1}=1$ atom per shot

when running the sequence without the Landau-Zener transfer, and approximately
half this amplitude after applying the pulse. As the signal of interest was
not the angular density profile but the integrated number of counts, we
corrected for this contribution by subtraction,
$N_{k_{\text{min}},k_{\text{max}}}=(N_{k_{\text{min}},k_{\text{max}}}^{\text{detected}}-\eta_{1}N_{m_{J}=+1})/(\eta_{0}\times\eta_{Q}\times\Omega_{\text{ROI}})$,
where $\Omega_{\text{ROI}}$ is the fraction of the sphere retained by the
windowing procedure described in the main text. and $\eta_{Q}$ is the detector
quantum efficiency. Such counts constitute about 10(5)% of the detection
events in the ROI and are not sufficient to account for the excess tail
strength. We hypothesize that the remnant counts are atoms transferred into
the $m_{J}=0$ state by non-ideal behaviour of the Stern-Gerlach pulses or
magnetic field switches.

#### Dependence of tail properties on the ROI

Table S1 shows the absolute and relative tail amplitude estimates
$\Lambda_{\textrm{fit}}$ and $\Lambda_{\textrm{fit}}/\Lambda_{\textrm{pred}}$,
respectively, that proceed from variations of the ROI and assumed quantum
efficiency in the analysis of the experimental data. It is seen that neither
the uncertainty in the collection efficiency nor the choice of elevation angle
cutoff $\phi_{c}$ have a significant effect on the findings.

QE | fit $r^{2}$ | $\Lambda_{\textrm{fit}}$ | $\Lambda_{\textrm{fit}}/\Lambda_{\textrm{pred}}$
---|---|---|---
0.05 | 0.83 | 0.2(0.1,0.3) | 6.8(4.5,9.1)
0.06 | 0.83 | 0.3(0.2,0.4) | 7.4(4.9,9.8)
0.07 | 0.83 | 0.3(0.2,0.4) | 7.8(5.2,10.5)
0.08 | 0.83 | 0.4(0.3,0.5) | 8.3(5.5,11)
0.09 | 0.83 | 0.5(0.3,0.6) | 8.7(5.8,11.6)
0.1 | 0.83 | 0.6(0.4,0.7) | 9.0(6.0,12.1)
0.11 | 0.83 | 0.6(0.4,0.8) | 9.4(6.2,12.5)
$\phi_{c}$ | fit $r^{2}$ | $\Lambda_{\textrm{fit}}$ | $\Lambda_{\textrm{fit}}/\Lambda_{\textrm{pred}}$
80∘ | 0.82 | 0.1(0.0,0.1) | 9.3(6.0,12.7)
70∘ | 0.86 | 0.2(0.1,0.3) | 9.2(6.4,12)
60∘ | 0.83 | 0.4(0.3,0.5) | 8.3(5.5,11)
Table S1: Tail predictions with modified ROI or assumptions about the value of
quantum efficiency (QE). A change in quantum efficiency (QE) presents by
changing the value of $\epsilon$ used in the prediction and also through the
factor of $N_{0}^{7/5}$ used to compute the condensate density $n_{0}$. These
effects partially cancel to produce a weak scaling of
$\Lambda_{\textrm{fit}}/\Lambda_{\textrm{pred}}$ with respect to $\epsilon$.
Re-running the analysis using different QE and the standard ROI with
$\phi_{c}=\pi/3$ yields fits that barely differ in the goodness-of-fit
criterion and present comparable results for $\Lambda_{\textrm{fit}}$. We also
find that the choice of collection area defined by the elevation angle cutoff
$\phi_{c}$ (while using the best QE value of 0.08) has a weak effect on the
result, but below statistical significance. Terms in brackets are the upper
and lower 95% confidence intervals. Neither the uncertainty in the collection
efficiency nor the choice of elevation angle cutoff $\phi_{c}$ have a
significant effect on the findings.

## Statistical details

In Figure S2 we illustrate how fixing the scaling exponent $\alpha$ when
fitting to a given dataset obscures the enormous covariance between $\alpha$
and the scale coefficient $C_{\alpha}$. The red diamonds in Figure S2 show the
fit exponent and amplitude obtained from a fit of Eqn. (13) to single data
sets, each representing a different BEC density. The variation in the fit
exponent is small (mean 4.2, standard deviation 0.4), and although the mean
$\alpha$ from the fits is only half a standard deviation (one standard-error
interval) different from 4, the corresponding difference in amplitude $C$
varies exponentially with $\alpha$, spanning over six orders of magnitude.

Figure S2: Illustrating the large systematic errors in a fit to the density
using a power law ansatz. If the fit to the density profile (Fig 1 in the main
text) ) is replaced by Eqn. (13) (also in main) and $\alpha$ and $C_{\alpha}$
used as fit parameters, it gives a wide variation in best-fit exponents and
scale coefficients between each data set (red diamonds). The red square shows
the mean fitted $\alpha$ and (geometric) mean $C_{\alpha}$ with standard
deviation (in $\alpha$) shown as error bars. Blue circles show the amplitude
coefficient $C_{\alpha}$ obtained from the same data sets with $\alpha$ fixed
at a number of values within the range of the free fit results. The choice of
$\alpha$ strongly determines the coefficient $C_{\alpha}$, but the error bars
(standard errors in fit parameters) are smaller than the markers in all cases.

The blue circles in Figure S2 show the amplitudes returned from fits to data
from all datasets when $\alpha$ is constrained to the set values as shown.
These set values are chosen within about one standard deviation of the mean
$\alpha$ from the two-parameter fits to separate datasets. These blue fits
scarcely differ in their goodness-of-fit criterion (the mean square error) and
so offer no obvious way to reconcile the expected distribution with these
divergent statistical conclusions. Furthermore, this is not reflected in the
error estimates in the fitting routines: The error bars representing the
uncertainty in parameters from the fixed-$\alpha$ fits are smaller than the
markers used in Figure S2. A linear fit reveals that
$d\log_{10}C_{\alpha}/d\alpha\approx 6.8$.

Ultimately, the question of whether the data is drawn from a power law at all
is not amenable to a decisive conclusion. For example, a log-normal
distribution can produce similarly accurate predictions to the power law (see
Figure 2 in the article) although there is no physical hypothesis that
predicts such a distribution. This points to one of the primary challenges
with power laws; the exponents are strongly entwined with the rate of
occurrence of rare events, which by definition are subject to large
statistical fluctuations and and thus subvert even the most meticulous
investigations. These problems with fitting power laws are ubiquitous, and
made more difficult by the small range of $k$ which are visible in the helium
experiments. In general, estimating the exponent of a purported power law is
difficult and requires data spanning several orders of magnitude in scale [84,
59, 60, 85], which are not present in either of the helium experiments to
date.

Nevertheless — the analysis of Fig. S2 does allow us to at least say that the
data is inconsistent with power laws outside the range $3.3<\alpha<4.7$.

## Simulations

#### Observables

Detailed data extracted from the simulations are given in Table. S3.

The total quantum depletion of the condensate $\delta_{B}$ is given by

$\delta_{B}=\frac{N_{B}}{N}=\frac{1}{N}\int
d^{3}{\mathbf{x}}\left\langle\real{\psi_{B}({\mathbf{x}},t){\widetilde{\psi}}^{*}_{B}({\mathbf{x}},t)}\right\rangle_{\rm
stoch.}$ (S1)

Stochastic averaging over all trajectories in the ensemble is denoted by
$\langle\cdot\rangle_{\rm stoch.}$. The density of the depleted particles is
evaluated by the standard positive-P expression,

$n_{B}({\mathbf{k}})=\left\langle\widehat{\Psi}^{\dagger}_{B}({\mathbf{k}},t)\widehat{\Psi}_{B}({\mathbf{k}},t)\right\rangle={\rm
Re}\left\langle{\widetilde{\psi}}^{*}_{B}({\mathbf{k}},t)\psi_{B}({\mathbf{k}},t)\right\rangle_{\rm
stoch.}.$ (S2)

The density of condensate

$n_{BEC}({\mathbf{k}})=\left(1-\delta_{B}\right)|\phi({\mathbf{k}},t)|^{2},$
(S3)

is augmented by the $1-\delta_{B}$ factor when calculating observables to
conserve overall particle number. The k-space fields here are normalised as

$\left[\begin{array}[]{c}\widehat{\Psi}_{B}({\mathbf{k}})\\\
\psi_{B}({\mathbf{k}})\\\
{\widetilde{\psi}}_{B}({\mathbf{k}})\end{array}\right]=\frac{1}{{(2\pi)^{3}}}\int
d^{3}{\mathbf{x}}\ e^{-i{\mathbf{k}}{\mathbf{x}}}\
\left[\begin{array}[]{c}\widehat{\Psi}_{B}({\mathbf{x}})\\\
\psi_{B}({\mathbf{x}})\\\
{\widetilde{\psi}}_{B}({\mathbf{x}})\end{array}\right]$ (S4)

Figure S3:  Details of initial state preparation, shown on the $N=455852$
case. Upper panel: Calibration data. Shown is the calculated contact $C_{\rm
sim}$ at the end of the initial state generation ($t=0$), as a function of the
quench ramp time $t_{\rm ramp}$. The horizontal dashed line shows the LDA
prediction $C=0.102$ pm-1, the blue data point the ensemble that was deemed to
agree, and was used for subsequent simulations for $t>0$. Lower panel:
Evolution and stabilisation of the contact. The blue shaded area shows the
duration of the ramp from $t_{\rm start}$ to $t_{\rm start}+t_{\rm ramp}$.
Orange shading denotes the error bars. ${\mathcal{S}}=4000$ trajectories in
all cases.

#### Initial ensemble

A procedure for generating the equilibrium state has been developed in the
Wigner representation [32, 86, 87, 88, 89]. Unfortunately we cannot use this
directly, nor a direct Wigner representation of the quasiparticles because the
condition required for correctness of the Wigner representation – that there
are more particles than modes – is very far from being met [90, 91, 33] (In
fact here we have about ${\mathcal{O}}(1000)$ particles in the depletion, and
${\mathcal{O}}(10^{6})$ modes.). It is also unclear how to translate a local
density formulation of depletion in a uniform section of gas to a positive-P
ensemble without introducing discontinuities. Instead we turned to dynamically
generating a state with the required quantum depletion.

We begin with a fully condensed ground state with
$\psi_{B}={\widetilde{\psi}}_{B}=0$ and
$\phi({\mathbf{x}})=\phi_{0}({\mathbf{x}})$, the ground state of the Gross-
Pitaevskii equation (21). The latter is obtained by imaginary time propagation
of the GPE augmented with an appropriately chosen chemical potential $\mu$
according to (8) which sets the central density.

Our first attempt to generate the equilibrated quantum depleted state thus
started with $\phi_{0}({\mathbf{x}})$ and then adiabatically ramped the
interaction from $g=0$ to the experimental value while evolving the equations
(21)–(22). This did not work for two reasons. Firstly, a very strong
collective oscillation was induced, since the width of the Thomas-Fermi ground
state depends strongly on $g$. Secondly, very long time evolutions succumb to
excessive noise in the positive-P simulation and produce a state that is too
noisy to be useful. Note that in the positive-P representation different
ensembles can represent the same state but exhibit very different noisiness
and practical usefulness [92, 93].

The second attempt began with the Gross-Pitaevskii ground state and the target
physical interaction strength $g$ in the GPE (21), while slowly ramping the
interaction strength in the Bogoliubov equations (22). This eliminates the
main oscillations in the mean-field evolution. The least unwanted nonadiabatic
disturbance occurs when $g$ is ramped only within the projected part of the
Bogoliubov evolution as per

$\displaystyle i\hbar\frac{d\phi}{dt}$ $\displaystyle=$
$\displaystyle{\mathcal{H}}(g,\phi)\phi=(\ref{GPE-eq})$ (S5a) $\displaystyle
i\hbar\frac{d\psi_{B}}{dt}$ $\displaystyle=$
$\displaystyle{\mathcal{H}}(g,\phi)\psi_{B}+{\mathcal{P}}_{\perp}\left\\{g_{B}|\phi|^{2}\psi_{B}+g_{B}\phi^{2}{\widetilde{\psi}}_{B}^{*}+\sqrt{-ig_{B}}\,\phi\,\xi({\mathbf{x}},t)\right\\}$
(S5b) $\displaystyle i\hbar\frac{d{\widetilde{\psi}}_{B}}{dt}$
$\displaystyle=$
$\displaystyle{\mathcal{H}}(g,\phi){\widetilde{\psi}}_{B}+{\mathcal{P}}_{\perp}\left\\{g_{B}|\phi|^{2}{\widetilde{\psi}}_{B}+g_{B}\phi^{2}\psi_{B}^{*}+\sqrt{-ig_{B}}\,\phi\,{\widetilde{\xi}}({\mathbf{x}},t)\right\\}.$
(S5c)

with $g_{B}(t)$. However, over adiabatic timescales, this still introduced far
too much noise in the Bogoliubov fields $\psi_{B}$, ${\widetilde{\psi}}_{B}$
to be useful.

To work around the problem, we take advantage of the fact that an
instantaneous quantum quench of a weakly interacting condensate produces
depletion with a similar time-integrated momentum profile $\int
dt\,n_{B}({\mathbf{k}},t)\propto 1/k^{4}$ to the ground state value, but with
a somewhat higher depletion, as described in detail in [61, 80]. Simulations
show that non-instantaneous but rapid quenches produce lower amounts of
quantum depletion, as shown in Figure S3 (left panel). To generate these
ensembles, we used the following ramp:

$g_{B}(t)=\left\\{\begin{array}[]{c@{\quad\text{when}\quad}l}\left(\frac{t-t{\rm
start}}{t_{\rm ramp}}\right)\,g&t<t_{\rm start}+t_{\rm ramp}\\\ g&t\geq t_{\rm
start}+t_{\rm ramp}\end{array}\right.$ (S6)

starting at negative $t_{\rm start}$, and evolving till $t=0$. Thus $t_{\rm
ramp}$ was the ramp time, and the remaining time: $|t_{\rm start}|-t_{\rm
ramp}$ an equilibration time. We carried out a calibration like that shown in
Figure S3 (left panel) for each set of cloud and trap parameters needed. Then
for the initial in situ ensemble in the trap we chose the ensemble generated
with the $t_{\rm ramp}$ that produces total quantum depletion in agreement
with the _in situ_ value (9) for the condensate ground state (like the blue
data point in Figure S3, left panel). The value of $t_{\rm start}$ was
$-101\mu$s for the weakly trapped cases ($\omega=45\times 425\times 425$ Hz
and $\bar{\omega}=201$ Hz) and $-198\mu$s for the strongly trapped cases
($\omega=71\times 902\times 895$ Hz and $\bar{\omega}=393$ Hz) detailed in
Table. S3. Figure S3 (right panel) shows an example of the evolution of the
contact in situ during this initial state generation. Properties of the
initial states are labelled (CT) and (ST) in Fig. 3 (a) (main text).

#### Determination of the contact and tail strengths

Much as in the experiment, the k-space density $n({\mathbf{k}})$ in the
simulations is very noisy in the asymptotic region of large $|{\mathbf{k}}|$,
and a lot of averaging is needed to extract the contact. In light of the
discussion on power laws, the exponent $\alpha$ was not fitted here, but we
did calculate both the fit coefficients $C_{\rm sim}$ on the function
$n_{B}(k)=C_{\rm sim}/k^{4}$ with $\alpha=4$ assumed and the count of
particles in the tails, $N_{k_{\rm min},k_{\rm max}}$ as in Table S2, in order
to have a systematic view of how correlated the two kinds of results are.

We proceed as follows: The simulations provide a density of depleted particles
$n_{B}({\mathbf{k}})$. Like with the experimental counts, we keep only density
that is far enough away ($\phi_{c}$, usually $60^{o}$) from the long axis of
the cloud. The simulated system is axially symmetric around the long axis. We
do not restrict counting to the vicinity of the vertical axis because here
there are no detector irregularities to avoid, and this allows us to improve
the signal to noise ratio.

The particle counts $N_{B}({\mathbf{k}})=n_{B}({\mathbf{k}}){/V}$ are finely
binned according to the absolute value of the momentum $k=|{\mathbf{k}}|$,
giving total bin count $N_{k}$. The field is simulated on a square lattice in
k space, in which each site corresponds to a k-space volume of $\Delta
V_{K}{=(2\pi)^{3}/V}$ with real space volume $V$. Therefore the total bin
volume $V_{k}$ is obtained by binning the site volumes, and the mean density
in each bin is $N_{k}/V_{k}$ so $n_{k}=(2\pi)^{3}N_{k}/V_{k}$. Each bin gives
an estimate of the corresponding apparent contact:

$c(k)={n_{k}k^{4}}=\frac{N_{k}}{V_{k}}\,k^{4}(2\pi)^{3}.$ (S7)

The statistical error estimate on $c(k)$, $\Delta c(k)$, is obtained by
averaging subensembles, then using the central limit theorem (CLT) on the
subensemble averages. Figure S4 shows example values and error estimates,
while also demonstrating the difficulties involved.

Figure S4:  The extraction of contact from simulation data – shown by example
of the $N=455852$ trap release (CE) simulation. Blue/cyan data: after
expansion at $t=689\mu$s, Red/orange: in situ initial condition at $t=0$. Each
data point corresponds to estimation of $c(k)$ given by (S7) using the density
of particles in the Bogoliubov field $\widehat{\Psi}_{B}({\mathbf{k}})$ in a
radial bin centred at $k=|{\mathbf{k}}|$. Error bars give the statistical
error from an ensemble of ${\mathcal{S}}=4000$ trajectories. Blue and red
points are the fitted part of the data. The vertical dashed line indicates
$k_{\rm inner}$, while the shaded region on the left shows the corresponding
scaled density $c_{\rm BEC}(k)=n_{\rm BEC}(k)k^{4}$ that would be calculated
from the condensate field $\phi({\mathbf{k}})$. The horizontal dashed lines
indicate the final estimates of $C_{\textrm{sim}}$ extracted from the data.

We can see that for a significant range of $k$ values the $c(k)$ estimate
gives fluctuations around a constant value. However, several difficulties in
extracting the overall trend are also evident. Low $k$ values are not
representative because the particles in the shaded region never emerge from
the condensate and are not measured in the experiment. Therefore we remove
data below a value $k_{\rm inner}$ from consideration. $k_{\rm inner}$ is
chosen such that it does not include any depletion that would significantly
overlap with the BEC in the final expanded cloud and be obscured by it.

This corresponds also to the energy at which the quasiparticle spectrum
becomes particle-like, since the mean field energy in the centre of the
trapped cloud is responsible for both effects.

High $k$ values, on the other hand, suffer from statistical error sufficiently
high as to make them useless, and so we also choose a maximum $k_{\rm outer}$
value, and only use $k\in[k_{\rm inner},k_{\rm outer}]$ to extract the contact
estimate. The extent of the useless high k region changes with time during the
simulation. To systematically adapt our fitting region to this, we choose
$k_{\rm outer}$ according to the calculated statistical error such that

$\Delta c(k,t)\leq\Delta_{\rm max}\quad\forall k\leq k_{\rm outer}(t),$ (S8)

and $\Delta_{\rm max}$ is chosen once for all times in a simulation. This
gives a time-dependent fitting region $k\in[k_{\rm inner},k_{\rm outer}(t)]$.
The final contact estimate $C_{\textrm{sim}}(t)$ is the mean of all $c(k,t)$
values in the fitting range. Error bars on $C_{\textrm{sim}}(t)$, $\Delta
C_{\textrm{sim}}(t)$, are obtained from a CLT estimate of the error in the
mean of the points used, after binning enough neighbouring points to encompass
the autocorrelation (i.e. so that neighbouring bins are uncorrelated). The
$\Delta_{\rm max}$ is chosen low enough so that the statistical error in
individual points $\Delta c(k,t)$ at the large $k$ end does not put off the
overall contact estimation $C_{\textrm{sim}}(t)\pm\Delta C_{\textrm{sim}}(t)$.
Operationally we choose the value of $\Delta_{\rm max}$ at which the error
introduced by adding successive high $k$ points introduces more uncertainty
than the reduction thanks to a larger ensemble.

$n_{o}N$ | $k_{\rm min}$ | $k_{\rm max}$ | $\phi_{c}$ | $N_{k_{\rm min},k_{\rm max}}$ | ratio $\mathcal{A}_{\rm sim}$
---|---|---|---|---|---
($10^{6}\mu$m-3) | ($\mu$m-1) | (deg) | final | final/in situ
Rapid release from trap (CE)
5.237(16) | 2.0 | 3.5 | 60 | 193(8) | 1.38(5)
12.09(6) | 2.25 | 3.5 | 60 | 362(10) | 1.43(4)
6.86(3) | 2.75 | 4.0 | 60 | 163(4) | 1.45(4)
19.90(8) | 3.0 | 4.0 |  40 |  423(8) |  1.37(3)
19.90(8) | 3.0 | 4.0 |  50 |  408(8) |  1.56(3)
19.90(8) | 3.0 | 4.0 | 60 | 368(7) | 1.80(4)
19.90(8) | 3.0 | 4.0 |  70 |  282(6) |  2.00(4)
19.90(8) | 3.0 | 4.0 |  80 |  157(5) |  2.19(6)
Rapid release, spherical (SE)
5.35(2) | 2.0 | 3.25 | 60 | 120(9) | 0.92(6)
12.23(6) | 2.0 | 3.0 | 60 | 256(11) | 1.05(4)
20.38(10) | 3.0 | 4.0 | 60 | 232(10) | 1.10(5)
Slow ramp down of trap (CS)
5.237(16) | 2.0 | 3.5 | 60 | 55(10) | 0.37(6)
12.09(6) | 2.0 | 3.5 | 60 | 103(21) | 0.31(5)
6.86(3) | 2.5 | 4.0 | 60 | 50(8) | 0.34(6)
19.90(0) | 2.5 | 4.0 | 60 | 181(14) | 0.54(4)
Table S2:  Tail strength data in simulations. Based on final times in the
simulations described in Table S3, referenced by the value of $n_{0}N$.
$N_{k_{\rm min},k_{\rm max}}$ is calculated as per Eqn (14) in the main text.
The ratio $\mathcal{A}_{\rm sim}$ of tail strength in the expanded cloud
compared to Tan theory predictions is obtained by dividing $N_{k_{\rm
min},k_{\rm max}}$ at the final time by its value in situ at $t=0$. $k_{\rm
max}$ is chosen to still contain all $4\pi$ steradians inside the square
lattice, while $k_{\rm min}$ to avoid overlap with the expanding condensate.

We also found that the relative factors between Eqn. (14) and the simulated
tails (i.e. $C_{\textrm{sim}}/\mathcal{C}$ and $\mathcal{A}_{\rm sim}$) depend
on choice of cutoff angle $\phi_{c}$. We find that apparent tail strengths
$C_{\textrm{sim}}$ and $\mathcal{A}_{\textrm{sim}}$ are larger for smaller
collection regions that are more tightly concentrated about the strong
trapping axis, whereas larger collection angles (that include areas closer to
the weak axis) produce lower apparent tail strengths. Data on
$\mathcal{A}_{\rm sim}$ are shown in Table S2 while $C_{\textrm{sim}}$ was
1.6(1), 1.9(2), and 2.2(3) times $\mathcal{C}$ for $\phi_{c}$ values of 60,
70, and 80 degrees, respectively.

Trap at $t=0$ | Peak density $n_{0}$ | $N$ | $n_{o}N$ | timescale | Time $t$ | $C_{\textrm{sim}}(t)$ | $C_{\textrm{sim}}(t)/C_{\textrm{sim}}(0)$ | $N_{B}(t)$ |  $k_{\rm inner}$–$k_{\rm outer}$ | $\Delta_{\rm max}$
---|---|---|---|---|---|---|---|---|---|---
$\omega$ (Hz) | at $t=0$ ($\mu$m-3) |  | ($10^{6}\mu$m-3) | ($\mu$s) | ($\mu$s) | (pm-1) |  |  | ($\mu$m-1) | ($\mu$m-1)
In situ initial state (CT,ST): |  |  |  |  |
$425\times 425\times 45$ | 16.82(5) | 311360 | 5.237(16) |  $t_{B}$=19.8 | 0 | 0.0260(13) | – | 1391(7) | 2.1–4.2 | 0.025
$425\times 425\times 45$ | 21.32(10) | 567180 | 12.09(6) |  $t_{B}$=15.8 | 0 | 0.060(5) | – | 2860(13) | 2.2–4.75 | 0.032
$902\times 895\times 71$ | 32.17(13) | 213293 | 6.86(3) |  $t_{B}$=8.5 | 0 | 0.034(4) | – | 1324(6) | 2.9–4.9 | 0.020
$902\times 895\times 71$ | 43.66(18) | 455852 | 19.90(8) |  $t_{B}$=9.3 | 0 | 0.100(3) | – | 3047(12) | 3.25–5.4 | 0.032
$201\times 201\times 201$ | 16.88(6) | 316766 | 5.35(2) |  $t_{B}$=19.8 | 0 | 0.025(3) | – | 1396(7) | 2.1–4.1 | 0.025
$201\times 201\times 201$ | 21.32(10) | 573650 | 12.23(6) |  $t_{B}$=14.2 | 0 | 0.060(9) | – | 2741(11) | 2.2–4.15 | 0.04
$393\times 393\times 393$ | 44.15(21) | 461514 | 20.38(10) |  $t_{B}$=7.7 | 0 | 0.101(5) | – | 3107(12) | 3.0–5.0 | 0.025
Rapid release from trap (CE,SE), $\tau_{\rm release}=37.5\mu$s: |  |  |  |  |
$425\times 425\times 45$ | 16.82(5) | 311360 | 5.237(16) | $\tau_{\rm release}$ | 1583 | 0.0386(22) | 1.49(11) | 910(20) | 2.1–2.8 | 0.020
$425\times 425\times 45$ | 21.32(10) | 567180 | 12.09(6) | $\tau_{\rm release}$ | 1346 | 0.090(7) | 1.51(17) | 1880(40) | 2.3–3.05 | 0.020
$902\times 895\times 71$ | 32.17(13) | 213293 | 6.86(3) | $\tau_{\rm release}$ | 689 | 0.057(5) | 1.69(14) | 950(20) | 2.8–3.7 | 0.03
$902\times 895\times 71$ | 43.66(18) | 455852 | 19.90(8) | $\tau_{\rm release}$ | 689 | 0.171(10) | 1.66(12) | 2300(40) | 3.1–4.1 | 0.06
$201\times 201\times 201$ | 16.82(5) | 316766 | 5.35(2) | $\tau_{\rm release}$ | 1900 | 0.025(3) | 1.02(15) | 780(30) | 1.9–2.65 | 0.025
$201\times 201\times 201$ | 21.32(10) | 573650 | 12.23(6) | $\tau_{\rm release}$ | 1900 | 0.066(7) | 1.11(20) | 1630(50) | 2.3–3.15 | 0.05
$393\times 393\times 393$ | 44.15(21) | 461514 | 20.38(10) | $\tau_{\rm release}$ | 973 | 0.110(5) | 1.09(8) | 1780(50) | 3.1–3.8 | 0.05
Slow ramp down of trap (CS): |  |  |  |  |
$425\times 425\times 45$ | 16.82(5) | 311360 | 5.237(16) |  $t_{\rm ramp}$=2375 | 2375 | 0.010(3) | 0.40(12) | 1750(140) | 1.9–3.45 | 0.030
$425\times 425\times 45$ | 21.32(10) | 567180 | 12.09(6) |  $t_{\rm ramp}$=2375 | 2375 | 0.020(3) | 0.35(7) | 2700(300) | 2.0–3.9 | 0.040
$902\times 895\times 71$ | 32.17(13) | 213293 | 6.86(3) |  $t_{\rm ramp}$=1216 | 1216 | 0.011(6) | 0.3(2) | 1520(150) | 2.5–4.2 | 0.040
$902\times 895\times 71$ | 43.66(18) | 455852 | 19.90(8) |  $t_{\rm ramp}$=1216 | 1216 | 0.044(5) | 0.45(6) | 3800(400) | 2.5–4.1 | 0.035
Instantaneous switch-off of the trap: |  |  |  |  |
$902\times 895\times 71$ | 43.66(18) | 455852 | 19.90(8) |  $t_{\rm ramp}$=0 | 689 | 0.0187(6) | 1.80(6) | 2380(40) | 3.1–4.05 | 0.05
Table S3:  Main simulation data and parameters for fitting of $k^{-4}$ tail
amplitude $C_{\rm sim}$. The time $t$ for which data are calculated is counted
relative to the start of the trap release. Abbreviations (CT,CE,…) as in
Figure 3a of the main text. The range $k_{\rm inner}$–$k_{\rm outer}$ here was
chosen as explained in the text and used to obtain the estimate and
uncertainty for $C_{\textrm{sim}}$. In all cases, ${\mathcal{S}}=4000$
trajectories averaged. $N_{B}$ is the number of Bogoliubov field particles as
per (S1).

## Toy two-mode model of escape

In a uniform gas in the Bogoliubov approximation, the $\widehat{a}_{k}$ mode
is coupled only to $\widehat{a}^{\dagger}_{-k}$ and the condensate. The
Bogoliubov-de Gennes equations for these modes can then be written [81]

$\frac{d}{dt}\left[\begin{array}[]{c}\widehat{a}_{k}(t)\\\
\widehat{a}^{\dagger}_{-k}(t)\end{array}\right]=-\frac{i}{\hbar}\left[\begin{array}[]{cc}\hbar^{2}k^{2}/2m+2gn(t)-\mu(t)&gn(t)\\\
-gn(t)&-\hbar^{2}k^{2}/2m-2gn(t)+\mu(t)\end{array}\right]\left[\begin{array}[]{c}\widehat{a}_{k}(t)\\\
\widehat{a}^{\dagger}_{-k}(t)\end{array}\right].$ (S9)

Notice the provision for a time-dependent background condensate density
$n(t)$. We stay in the real-particle basis rather than the Bogoliubov
quasiparticles $\widehat{b}_{k}$ which allows to avoid calculating time-
dependent changes of the coefficients $u_{k}$, $v_{k}$. The chemical potential
is $\mu(t)=gn(t)$ [81].

The equations (S9) can be used as input for equations of motion of the low
order moments $\rho(k)=\langle\widehat{a}^{\dagger}_{\pm k}\widehat{a}_{\pm
k}\rangle$, $A(k)=\langle\widehat{a}_{k}\widehat{a}_{-k}\rangle$, and
$A^{*}(k)=\langle\widehat{a}^{\dagger}_{k}\widehat{a}^{\dagger}_{-k}\rangle$,
for example
$d\rho(k)/dt=\langle(d\widehat{a}^{\dagger}_{k}/dt)\widehat{a}_{k}\rangle+\langle\widehat{a}^{\dagger}_{k}(d\widehat{a}_{k}/dt)\rangle$.
Assuming equal initial occupations $\rho(k,0)=\rho(-k,0)$, one obtains an
evolution equation for two coupled quantities $\rho(k,t)$ and
$A(k,t)=A_{r}(k)+iA_{i}(k)$. It is

$\frac{d}{dt}\left[\begin{array}[]{c}\rho(k,t)\\\ A_{r}(k,t)\\\
A_{i}(k,t)\end{array}\right]=\frac{1}{\hbar}\left[\begin{array}[]{ccc}0&0&-2gn(t)\\\
0&0&2(\hbar^{2}k^{2}/2m-gn(t))\\\
-2gn(t)&-2(\hbar^{2}k^{2}/2m-gn(t))&0\end{array}\right]\left[\begin{array}[]{c}\rho(k,t)\\\
A_{r}(k,t)\\\ A_{i}(k,t)\end{array}\right]+\left[\begin{array}[]{c}0\\\ 0\\\
-gn(t)\end{array}\right].$ (S10)

The initial conditions are $\rho(k,0)=v_{k}^{2}$, $A(k,0)=v_{k}u_{k}$,
corresponding to the Bogoliubov ground state. Taking $n(t)$ constant one
obtains the solution (16).

A more realistic model is obtained by estimating the $n(t)$ that particles in
the $k$ mode experience as the trap is released and the condensate expands.
Hydrodynamic expansion of a condensate with a Thomas-Fermi profile suddenly
released from the trap [94] gives a self-similar expansion with widths
approximated by $W_{j}(t)=W_{j}(0)\sqrt{1+(\omega_{j}t)^{2}}$. The initial
widths being the Thomas Fermi radii in the trap
$R_{j}=(1/\omega_{j})\sqrt{2gn_{0})/m}$, $j=\\{x,y,z\\}$. From this, and
assuming $\omega_{y}=\omega_{z}=\omega_{\perp}=\bar{\omega}\lambda^{1/3}$
(approximately true in our setup) the condensate density at any point in the
cloud decays as

$n_{c}(r,t)=\frac{n_{c}(r,0)}{(1+\omega_{\perp}^{2}t^{2})\sqrt{1+\omega_{x}^{2}t^{2}}},\qquad\text{where}\qquad
n_{c}(r,0)=\left\\{\begin{array}[]{@{\ }l@{\text{\quad if
}}l}n_{0}\left[1-(r/R_{\perp})^{2}\right]&r<R_{\perp}\\\ 0&r\geq
R_{\perp}\end{array}\right.$ (S11)

Escape is much easier, however, for atoms near the edge of the cloud than in
the centre, because they can avoid reabsorption by faster free-flight from the
high density region. We categorise the initial position of the studied section
of the gas with the quantity

$R_{0}=\frac{\sqrt{y^{2}+z^{2}}}{R_{\perp}}=\frac{\omega_{\perp}\sqrt{y^{2}+z^{2}}}{\sqrt{2gn_{0}/m}}.$
(S12)

The initial density is then

$n(0)=\left\\{\begin{array}[]{@{\ }l@{\text{\quad if
}}l}n_{0}(1-R_{0}^{2})&R_{0}\leq 1\\\ 0&R_{0}>1\end{array}\right.$ (S13)

The motion of the atoms is estimated by assuming they travel in the outward
direction at velocity $\hbar k/m$ relative to the hydrodynamic expansion
(which also accelerates them to some degree). Hence their position is
estimated as

$r(t)=R_{0}R_{\perp}\sqrt{1+\omega_{\perp}^{2}t^{2}}+\frac{\hbar kt}{m}.$
(S14)

The condensate at this position comes from hydrodynamic expansion of the
initial condensate at position

$r_{c0}(t)=\frac{r(t)}{\sqrt{1+\omega_{\perp}^{2}t^{2}}}=R_{0}R_{\perp}+\frac{\hbar
kt}{m\sqrt{1+\omega_{\perp}^{2}t^{2}}}.$ (S15)

The local density at position $r$ at time $t$ is then $n_{c}(r_{c0}(t),t)$
according to (S11).

Finally, to roughly include ramp time in the analysis, consider the density in
a trapped gas that adiabatically follows the trap evolution
$\omega_{j}(t)=\omega_{j}(0)e^{-t/\tau_{\rm release}}$. The overall density
for a piece of gas that stays in the same relative position in the cloud will
then follow

$n_{c}^{\rm adiabatic}(t)=n(0)e^{-3t/\tau_{\rm release}}.$ (S16)

To incorporate this effect into the model it is combined (somewhat ad hoc)
with the expression $n_{c}(r_{c0}(t),t)$ for instantaneous release to give the
final estimate of local density for the $k$ momentum escaping atoms used in
(S10) as

$n(t)=n(0)e^{-3t/\tau_{\rm release}}+(1-e^{-3t/\tau_{\rm release}})\
n_{c}(r_{c0}(t),t).$ (S17)

with the help of (S11) and (S15).