### Experimental measurements

Our experimental sequence began with BECs consisting of between \(2\times 10^5\) and \(5 \times 10^5\) \(^4\) \(\hbox {He}^*\) 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 _{\mathrm{release}}}\approx 38\,\upmu\)s. The condensates were allowed to expand for 2 ms 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\(^*~\)are challenging because the absence of a known Feshbach resonance precludes control over the contact \({\mathscr {C}}\propto ((a N_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 Eq. (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. Eq. (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.8 eV) internal energy^{50} of the metastable \(2^{3}S_1\) state, \(\hbox {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^{49,53,54,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 Fig. 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~\mu {{\text {m}}}\) and \(36~\mu {{\text {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}

$$\begin{aligned} \frac{n_T(k)}{{(2\pi )^3}} =\frac{N_T}{\zeta (3)} ~\left( \frac{\lambda _{dB}}{2\pi }\right) ^3 g_{3/2}\left( \exp \left( -\frac{k^2 \lambda _{dB}^2}{4\pi }\right) \right) , \end{aligned}$$

(12)

wherein the thermal de Broglie wavelength \(\lambda _{dB} = \sqrt{2\pi \hbar ^2/(m k_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^{\mathrm{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.

### 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 (Eq. (12)) with a power-law term as per

$$\begin{aligned} {n(k)} = n_T(k) + \frac{C_{\alpha }}{k^{\alpha }}, \end{aligned}$$

(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 Eq. (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\le \alpha \le 4.6\), but the data at hand cannot precisely determine the exponent \(\alpha\) (nor \(C_\alpha\)), as detailed in the supplement.

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_{{\text {min}}}, k_{{\text {max}}})\),

$$\begin{aligned} N_{k_{{\text {min}}},k_{{\text {max}}}} =\frac{{\mathscr {C}}}{2\pi ^2}\left( \frac{1}{k_{{\text {min}}}}-\frac{1}{k_{{\text {max}}}}\right) . \end{aligned}$$

(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_{{\text {min}}}\) and \(k_{{\text {max}}}\), Eq. (14) has the form

$$\begin{aligned} N_{k_{{\text {min}}},k_{{\text {max}}}} = \Lambda N_0n_0, \end{aligned}$$

(15)

(c.f. Eq. (10)). We can thus test Eq. (15) directly by measuring the number of counts detected in the interval \((k_{{\text {min}}},k_{{\text {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_0n_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_{{\text {min}}},k_{{\text {max}}}} = \Lambda _{{\text {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_0N_0\), and against the high energy tails being due to some technical noise. The correlation coefficient between the variables \(N_{k_{\mathrm{min}},k_{\mathrm{max}}}\) and \(N_0n_0\propto (N_0^7{\bar{\omega }}^6)^{1/5}\) is 0.9. We conclude that the product \(N_0n_0\) is a predictor of the high energy population, which is consistent with Eq. (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~,p = 0.54\)), as is the central density \(n_0\) alone (\(r^2 = 0.4~,p = 0.04\)). Accordingly, the particular nonlinear scaling of detected counts with the predictor \(N_0n_0\) is concordant with the tails’ originating in the quantum depletion, and inconsistent with any technical noise that we know of.

The gradient \(\Lambda _{{\text {fit}}}\) is of particular interest because it can be predicted using Eq. (14). Given a region of interest (ROI) over which we count atoms, one can calculate \(\Lambda _{{\text {pred}}} = 32 \varepsilon a^2(k_{{{\text {min}}}}^{-1}-k_{{{\text {max}}}}^{-1})/7\), where \(\varepsilon\) is the total detection efficiency. In our experiment, \(\varepsilon \approx 0.23(5)\%\) (see “Methods” section for details). In comparison with the predicted value \(\Lambda _\mathrm {pred} = 2.7(6) \times 10^{-7}\) (units of \(\upmu {{\text {m}}}^3\)/atom), we find that the empirical fit disagrees with the predicted slope by a factor of \({\mathscr {A}}_\mathrm {exp}=\Lambda _{{\text {fit}}}/\Lambda _{{\text {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_0N_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 \le \alpha \le 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_0N_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 Fig. 2b,c, as a counterpoint to power law fits over \(k_{\mathrm{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

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. 3a) 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. 3a), 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. Figure 3b shows the time evolution of the tail amplitude \(C_{\mathrm{sim}}\) extracted from a \(n(k)=C_{\mathrm{sim}}/k^4\) fit to the simulated density. In this configuration the steady-state value of the momentum tails was a factor of \(C_{\mathrm{sim}}/{\mathscr {C}}=\)1.64(9) above the predictions of Eq. (10). An analysis of the occupation of the tails according to (14), gives very similar factors \({\mathscr {A}}_{\mathrm{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\,\upmu \hbox {s}\) in the centre of the trapped cloud is comparable to the trap release time \(\tau _{\mathrm{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 Ref.^{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 (Eq. 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_{\mathrm{sim}}\) as well as the the tail strength \(N_{k_{\mathrm{min}},k_{\mathrm{max}}}\) from (14) decreased roughly as predicted—see the dashed line in Fig. 3b and Supplementary Table S2, strongly supporting the hypothesis that adiabaticity is needed for agreement with the results of Ref.^{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 Eq. (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 Ref.^{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~\mu {{\text {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^{44,62,63}.

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)\rightarrow \approx {\mathscr {C}}k/(k^2-2gn_0m/\hbar ^2)^{5/2}\). This leads, for example, to a doubling of the apparent contact \(C_{\mathrm{sim}}\) at \(k\approx 6/\upmu\)m for clouds with \(n_0=39\,\upmu \mathrm{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_{{\text {sim}}}\) and \({\mathscr {A}}_{\mathrm{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 {{\text {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_{\mathrm{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'<n_0\) at \(t=0\) and then stays constant, has the occupation of each *k* mode evolve as

$$\begin{aligned} \rho (k,t) = \rho (k,0) – \frac{gn[\varepsilon _0(k)^2-\varepsilon (k)^2]}{4\varepsilon (k)^2\varepsilon _0(k)}\,[1-\cos 2\varepsilon (k) t]. \end{aligned}$$

(16)

Here, \(\varepsilon _0(k)\) and \(\varepsilon (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. 4a. 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\varepsilon (k)\).

Figures 4b–e show the behaviour of a more careful toy model given by Eqs. (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 _{\mathrm{release}}\lesssim 100\,\upmu \hbox {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.