Resonance triplet dynamics in the quenched unitary Bose gas

The quenched unitary Bose gas is a paradigmatic example of a strongly interacting out-of-equilibrium quantum system, whose dynamics become difficult to describe theoretically due to the growth of non-Gaussian quantum correlations. We develop a conserving many-body theory capable of capturing these effects, allowing us to model the post-quench dynamics in the previously inaccessible time regime where the gas departs from the universal prethermal stage. Our results show that this departure is driven by the growth of strong lossless three-body correlations, rather than atomic losses, thus framing the heating of the gas in this regime as a fully coherent phenomenon. We uncover the specific few-body scattering processes that affect this heating, and show that the expected connection between the two-body and three-body contacts and the tail of the momentum distribution is obscured following the prethermal stage, explaining the absence of this connection in experiments. Our general framework, which reframes the dynamics of unitary quantum systems in terms of explicit connections to microscopic physics, can be broadly applied to any quantum system containing strong few-body correlations.

Introduction.-It is a well known fact that the quantum many-body problem scales exponentially with the number of particles.Hence, despite the rapid rise in available computing power over the last few decades, exact descriptions typically remain out of computational reach [1,2].In many cases however, simplifying assumptions can be made about the underlying processes of the system, namely that they obey Gaussian statistics and can be described by long-lived quasiparticles and collective motions of the medium using just a few degrees of freedom [3].Such simplifications are at the heart of our understanding of a wide range of quantum phenomena, ranging from Fermi liquids and polarons to superfluids and conventional superconductors [4][5][6][7][8][9][10].
Due to the large degree of control over the interaction strength, ultracold atomic gases provide a versatile quantum simulator for probing non-Gaussian physics [39][40][41].Additionally, the underlying microscopic physics can be encoded exactly in many-body models of these systems, enabling quantitative comparisons between theory and experiment [42,43].At the frontier of such approaches is the description of ultracold atomic many-body systems featuring non-perturbative few-body effects [44][45][46][47][48][49].* Corresponding author: j.v.d.kraats@tue.nlHere, a series of experiments exploring the quench of a degenerate Bose gas to the unitary regime [50][51][52][53], where interactions are as strong as allowed by quantum mechanics, raise important fundamental questions regarding the interplay of integrability, ergodicity, and few-body correlations in strongly correlated systems far out of equilibrium.Immediately following the quench, the dynamics are Gaussian and integrable in nature, resulting in the formation of a universal prethermal stage [54][55][56][57].At later times, non-Gaussian correlations develop and integrability is broken, facilitating the transition of the system towards a global equilibrium.While this behavior is well characterized for weakly interacting systems [54], the analogous breaking of integrability for quenches to the strongly interacting regime poses a highly nontrivial theoretical problem due to the appearance of nonperturbative phenomena such as strong few-body scattering, bound states of the medium [57][58][59], and rapid threebody losses [60].Furthermore, the quench is associated with the formation of an infinite number of three-body bound Efimov states, whose role and impact remains a subject of active research [59,[61][62][63][64].
In this Letter, we utilize the method of cumulants to formulate a conserving theory of the quenched unitary Bose gas, capturing the non-Gaussian correlations forming in the intermediate time regime that follows the prethermal stage but preceeds the crossover into a thermal gas.Through direct comparisons with experiment, we characterize the prethermal departure as a consequence of the growth of strong few-body correlations, rather than incoherent particle losses, which broaden the momentum distribution and result in lossless heating of the system.Then, motivated by the absence of a powerlaw tail in the asymptotics of the momentum distributions observed in Refs.[50,53], we examine how this expected behavior set by the universal contact relations [65][66][67][68][69][70][71] is obscured as the system departs the prethermal arXiv:2309.04587v2[cond-mat.quant-gas]28 Mar 2024 regime.
Model.-We consider a system of N identical bosons of mass m in a cubic volume V , which occupy single-particle states with momentum k annihilated by operators âk .Two such atoms may couple to form a energeticallyclosed channel molecule with center of mass momentum 2k and internal energy ν > 0, annihilated by the molecular operator b2k .Crucially, we neglect any direct interaction between free atoms, such that all scattering processes in the energetically-open channel are mediated by this molecular state.The resulting Hamiltonian reads [72][73][74], with k ≡ |k| .The interaction coupling the atomic and molecular states is modelled by a separable potential with strength g and step-function form factor ζ(2q) ≡ θ(Λ − q), where θ is the Heaviside step function such that Λ represents a cut-off on the relative two-body momentum.By an analytic solution of the two-body problem, the model parameters can be linked to physical quantities via the renormalization relations [42,75,76], where a is the s-wave scattering length, R * is the characteristic length scale set by the atom-molecule transition rate [77], and the momentum cut-off is related to the van der Waals length as Λ ∼ 1/r vdW .At unitarity (a −1 → 0), the dressed molecular state becomes degenerate with the scattering threshold correspondent with a Feshbach resonance [40].We express all observables in the Fermi scales n /2m and t n = ℏ/E n , where n = N/V is the gas density.In a many-body context the relative importance of the molecular state is quantified by the resonance width R * k n [78].
For both the atomic and molecular fields we adopt the U(1) symmetry breaking picture of a Bose-Einstein condensate (BEC), where the singlets ⟨a k ⟩ = δ k0 √ V ψ a and ⟨b k ⟩ = δ k0 √ V ψ m are described by atomic and molecular condensate wave functions ψ a and ψ m respectively.To model the quench scenario, we assume that for t < 0 the gas is an ideal atomic BEC, such that |ψ a | 2 = n.At t = 0, the system is instantaneously quenched to unitarity such that na 3 ≫ 1.Then, the far out-of-equilibrium BEC state suffers quantum depletion and sequentially generates higher order correlations in the gas, which we track using a cumulant expansion [51,[79][80][81][82][83].Truncating the expansion at the level of two-body atomic correlations we obtain the model of Refs.[43,75], with vanishing background scattering length.This resonance doublet model includes the cumulants (denoted with subscript c), which describe the single-particle momentum distribution and pairing field respectively.Its integrable equations of motion are equivalent to the two-channel Hartree-Fock-Bogoliubov equations [3,43].In this work, we introduce an extension of this model referred to as the resonance triplet model.Here we include the additional molecular and mixed-channel cumulants, , and the non-Gaussian tripling field, In the short-range or vacuum limits, the three-body correlations κ am k and R a k,q obey a coupled-channel threebody Schrödinger equation and thus represent closed and open-channel components of the three-body wave function [76].As shown in Ref. [59], R a k,q becomes macroscopically occupied following the quench and furnishes an order parameter signalling the formation of a Bose-Einstein condensate of Efimovian triples, which motivates its inclusion in the theory.The cumulants χ k and n m k ensure conservation of energy, which is possible in the resonance triplet model due to the cubic form of the Hamiltonian (Eq.( 1)).This is a crucial difference between the present model and the single-channel triplet model of Ref. [57], where an inclusion of three-body correlations induces an unphysical leak of energy to the tail of the momentum distribution, thereby prohibiting that model from studying the intermediate time window considered in this Letter.
Prethermal departure.-Webegin by examining the dynamics of the momentum distribution in Fig. 1, comparing against experimental results to look for clues about the causes of the prethermal departure and the nature of the intermediate time regime.Since higher-order correlations develop sequentially following the quench, the doublet and triplet models are initially equivalent [80].In this early-time regime, atomic pair excitations are directly generated by condensed molecules, as depicted diagrammatically in Fig. 2(a).As shown in Ref. [75] Following the initial excitation growth, a momentum dependent plateau is reached, signifying the quasi-steady prethermal stage characterized by approximate equilibration of macroscopic observables and emblematic of integrable dynamics [54,57].Our key finding is that the FIG. 1. Broadening of the single-particle momentum distribution following the quench.In Figs . (a-c) we plot the dynamics of the excited state population n a k for a set of momenta.Resonance doublet(triplet) model results shown with dash-dotted(solid) lines, matched with the normalization of the experimental data of Ref. [53], shown with grey circles.In (d) we show in a similar format the dynamics of the average kinetic energy per particle ⟨ϵ⟩ k for the broad resonance case R * kn = 0.3, comparing with the experimental data of Ref. [52].Black dotted line shows the expected T ∼ t 2/13 scaling for loss induced heating [52].
rapid depletion of the condensate for broad resonances spurs the growth of R a and hence of non-Gaussian threebody correlations, leading to a momentum dependent departure from the prethermal stage visible in Fig. 1(a-c).The shape of the departure shows remarkable qualitative agreement with the experiment of Ref. [53], with just slightly slower growth of excitations due to the mismatch in resonance width.In contrast, the integrable resonance doublet model remains in the prethermal stage at all times, consistent with the single-channel doublet model studied in Ref. [57].
Physically, the broadening of the momentum distribution shown in Fig. 1 arises from the significant interaction energy injected into the system by the quench, which is gradually converted into kinetic energy via two and threebody scattering and thus results in lossless correlationinduced, rather than recombination induced, heating.To quantify this lossless heating, we follow experiment [52], and examine the averaged kinetic energy per particle [57], The restriction to a maximum momentum k accounts for the limited resolution in the experiment [52].From Fig. 1(d) it is clear that the Gaussian statistics of the resonance doublet model are unable to capture the experimental departure from the prethermal plateau that occurs by t ∼ t n .The resonance triplet model does capture this departure, following the experimental results until the crossover where the gas was experimentally found to pass from degenerate to nondegenerate regimes.This agreement further strengthens the interpretation of the dynamics at intermediate times as correlation (rather than loss) dominated.
Correlations in the intermediate time regime.-Havingprovided evidence for the lossless nature of the intermediate time regime, we now study the dominant processes in the system during this time.The difference between the resonance doublet and resonance triplet dynamics in Fig. 1 arises predominantly from a distinct non-Gaussian process in the resonance triplet model, shown diagrammatically in Fig. 2(b) [76].Together, the processes in Fig. 2 seed the post-quench growth of short-range atomic two-and three-body correlation functions, given as, Here d = ψ(0) ψ(0) and t = ψ(0) ψ(0) ψ(0) represent the local dimer and trimer field respectively, with the field operator ψ(0) = (1/ √ V ) k âk .Due to their short-range nature, the value of the correlation functions is expected to be fully specified by universal relations derived from few-body physics.Hence they can be directly related to the two-body and three-body contact densities C 2 and C 3 , which form a measure of the probability for finding two and three particles in close proximity [65][66][67]70].In our two-channel model, we derive the following relations FIG. 3. Dynamics of two-body and three-body contact densities C2 in (a) and C3 in (b), for different values of the resonance width R * kn.Results from the resonance doublet(triplet) model shown with dash-dotted(solid) lines.For comparison we show with black dashed lines the linear early time growth of C2 as derived in Ref. [84] for broad resonances and the long time constant value of C2 from Ref. [85].For C3 we show the range of quadratic early time growths derived in Ref. [63], which correspond to maximum (kn/κ * ∼ 1) and minimum (κ * ≪ kn) enhancement of C3 due to the Efimov effect, again in the broad resonance limit.
from effective field theory [76], Here H, J, H ′ are known log-periodic functions of the Efimovian binding wavenumber κ * and the cut-off Λ [70].
As shown in Fig. 3, the prethermal departure is indeed correlated with a significant increase in C 3 , indicating the introduction of strong non-Gaussian three-body correlations in this regime.Simultaneously, the influence of surrounding particles on clustered pairs leads to a significant decrease of the two-body contact [57].As shown in Ref. [59], the value of C 2 and C 3 following the quench is determined predominantly by the magnitude of macroscopic order parameters associated with condensation of pairs and triples.These findings are confirmed in our model by the development of a substantial fraction of condensed triples for t > t n [76].We note that the increase of C 3 for broad resonances is connected with an increased overlap of the size of the Efimov trimer, quantified by the binding wavenumber κ * , and the Fermi scale k n .Consistent with Figs.1(a-c), the value of C 3 is decreased for narrow resonances, where κ * ≪ k n represents the previously unexplored opposite limit to the universal regime κ * ≫ k n examined in Ref. [59].
Asymptotics in the intermediate time regime.-Wenow focus on the large-momentum tail of n a k , whose expected behavior due to universal relations derived from few-body physics was not observed in experiments [50,53].Specifically, the following asymptotic behavior FIG. 4. Dynamics of the tail of the single-particle momentum distribution for the broad resonance R * kn = 0.3.In the left(right) panels we compare the resonance doublet(triplet) models with the asymptotic prediction in Eq. ( 10), shown with a black dashed line.For the sake of comparison the doublet results are replotted in the right-hand panels.For t/tn = 1.0 and 2.0 we also compare with the experimental data of Ref. [53], shown with the green circles.
is predicted at thermal equilibrium, where F (k) is a log-periodic function specified by κ * [70].While Eq. ( 10) provides a possible route to extract the values of the contacts from n a k , a 1/k 4 tail has so far not been observed in quench experiments [50,53].Fits to Eq. (10) were made in Ref. [71], but found values for C 2 considerably larger than expected from theoretical calculations [57,85,86].Additionally it has been shown that in lower dimensional systems the scaling in Eq. ( 10) can be disrupted by non-local correlations [87].
Motivated by this disagreement between existing theory and experiment, we compare in Fig. 4 the singleparticle momentum distribution obtained from our models with the asymptotic prediction in Eq. (10).In the resonance doublet model, where all dynamics are Gaussian, C 3 vanishes trivially and one can show analytically that the expected 1/k 4 power law due to C 2 is always present [76].In the resonance triplet model, we observe that the non-Gaussian processes responsible for the departure from the prethermal stage for t ≳ t n damp the oscillatory behavior of n a k , consistent with the findings of Ref. [54] at weak interactions and in agreement with the experimental data of Ref. [53].At the same time, the significant increase of excitations for k > k n obscures the expected power laws for all the examined momenta k/k n ≤ 4. Hence, once the gas has exited the prethermal stage, the asymptotic expansion in Eq. (10) no longer captures the momentum distribution over the considered range.As experimental results are lacking beyond this range due to poor signal-to-noise ratio [50], our findings suggest that significant caution should be exercised when fitting Eq. ( 10) out of equilibrium, and explain why a power-law tail was not seen in Refs.[50,71].It is important to note however that our results do not invalidate Eq. ( 10), but rather push its possible applicability to larger momenta.For sufficiently large values of k/k n , our distributions do converge consistently with the value of C 2 obtained from Eq. ( 8), but also exhibit strong finite range features due to the density regimes considered [76].We have confirmed that the distributions in Fig. 4 are density independent.
Conclusion-In this Letter, we have shown how a conserving many-body model constructed from a selection of microscopically relevant Gaussian and non-Gaussian correlations, is able to elucidate the dynamics of quenched unitary Bose gases in the time window succeeding the integrable dynamics associated with the prethermal stage.In the future, this general framework can be used in studying the wide-array of other quantum systems that exhibit non-Gaussian physics, including for example strongly interacting ultracold mixtures and polarons [29,48,49,88,89], trions in semiconductors [90][91][92][93][94], nuclear matter [95][96][97][98], and Rydberg atom arrays [99].
In this section we provide further detail regarding the two-body interaction used in the Hamiltonian in Eq. ( 1) in the main text, see also Ref. [S1].We consider the stationary two-body scattering problem in the basis of states |k, −k⟩ describing two counter propagating plane waves with relative momentum k.In addition we define an energetically open and energetically closed channel, correspondent with distinct internal states of the particles, and fix the internal energy of the open channel to zero.We adopt the Feshbach formalism [S2], defining projection operators P and Q onto the open channel and closed channel subspaces respectively, with properties P P = P, Q Q = Q, P + Q = I and P Q = 0.A general two-body operator Ô may be decomposed as, We consider the atomic projection P |ψ⟩ = |ψ P ⟩, where Ĥ |ψ⟩ = E |ψ⟩.By inserting projection operators one can rewrite to obtain, where VP Q = VQP is the two-body interaction coupling the open and closed channel subspaces, and ĜQ (E) = (E − ĤQQ ) −1 the closed-channel Green's function.Equation (S2) motivates the introduction of an effective potential, Here we have made the important assumption that the resonant bound state |ϕ⟩ is the sole eigenstate of the closedchannel Hamiltonian, such that ĜQ (E) = (E − ν) −1 |ϕ⟩ ⟨ϕ|.This isolated resonance approximation is valid provided that the state |ϕ⟩ is well separated from other closed-channel states and its binding energy is large compared to the scattering energy.From Eq. S3 we formulate the open-channel component of the two-body transition matrix, which we evaluate for infinitesimally complex energies to skirt the poles on the real axis, and where  S4) can be formally solved to obtain, where we have defined g = √ 2β ⟨ϕ|ζ⟩ as the coupling strength.The matrix elements of the uncoupled two-body Green's function Ĝ2B 0 (z) = (z − Ĥ0 ) −1 are evaluated as, which holds to second order in k/Λ.Here we have introduced the scattering length a and resonance length R * as, which may be inverted to obtain the renormalization relations quoted in Eq. ( 2) in the main text.The above expression for R * follows the definition in Ref. [S3].Finally we can relate Λ to the range of the potential by considering the energy E D of the Feshbach dimer, signified by a pole in the transition matrix for a > 0. If we define a dimer wavenumber to second order in κ D /Λ.In systems with van der Waals interactions, the dimer binding energy near resonance is universally determined as , where ā = 0.995978 r vdW is the mean scattering length [S5].Substituting into Eq.(S9) one finds, where the last equality holds in the resonant limit a ≫ ā, R * , which is always true in our calculations.The universal zero-range limit is achieved when Λ far exceeds all other momentum scales, which has the important implication in our model that the ratio r vdW /R * should be small to avoid non-universal finite range effects.For all computations in the main text we set Λ such that nr 3 vdW = 5.8 • 10 −7 , which is sufficiently dilute for our results to be universal in the regime of resonance strengths we examine.

II. MANY-BODY MODEL
In this section we provide further detail regarding our many-body formalism and the cumulant equations of motion.

A. Molecular operators
If we define ĉk as the operator that annihilates a particle with momentum k in the closed spin-channel, then the composite molecular operator is defined as where ϕ(q) ≡ ⟨q|ϕ⟩.The operator bk obeys canonical commutation relations provided that the state |ϕ⟩ is strongly localized with respect to the density scales of the gas [S1].Since we assume that the state |ϕ⟩ is the only eigenstate of the closed channel Hamiltonian H QQ , Eq. S11 may be inverted as ĉ k 2 +q ĉ k 2 −q = √ 2ϕ(q) bk .Hence we can convert any expectation value containing an even-numbered product of single-particle closed-channel operators into a product of molecular operators, whilst odd-numbered products vanish.

B. Cumulant equations of motion
To track the post quench dynamics, we solve the Heisenberg equations of motion , with Ôp an arbitrary p-body operator and Ĥ given by Eq. ( 1 where α, β can represent either atomic or molecular operators.The resulting cumulant expansion transforms Eq. (S12) to a system of coupled ordinary differential equations for all the individual cumulants.Due to the cubic structure of Ĥ, a cumulant containing n operators couples to cumulants containing at most n + 1 operators.This is a major simplification compared to the more commonly used quartic single-channel Hamiltonian, where coupling occurs also to cumulants with n + 2 operators [S9].
The dynamics of the atomic and molecular condensate wave functions are given by the two-channel Gross-Pitaevskii equations, Upon including the atomic excitation density and atomic pairing field as defined in Eq. ( 3) in the main text, one obtains the two-channel Hartree-Fock-Bogoliubov equations, Here ε a k = ℏ 2 k 2 /(2m) gives the atomic kinetic energy.In the equation of motion for κ a k , we recognize the Bose enhancement factors (1 + n a k ), which alter the coupling to excited states with non-zero occupation.In the resonance triplet model, the additional cumulants introduced in Eqs. ( 4) and ( 5) in the main text obey the equations of motion, Here Ŝ is a three-body symmetrization operator defined as Ŝ = 1 + P+ + P− , where P+(−) define (anti)symmetric cyclic permutation operators of the three particle indices.We have also defined the molecular kinetic energy The molecular and anomalous excitation densities n m k and χ k are necessary inclusions in the model to ensure that the total system energy is always conserved.If these cumulants are excluded, the growth of the three-body cumulant κ am k leads to an unphysical increase in the total kinetic energy, similar to the violation observed in single-channel triplet models [S9].As we show in Sec.IV, the cumulants n m k and χ k are actually fully determined from short-range physics, meaning that they have predominant amplitude in the large momentum modes k ∼ Λ.This necessarily implies that the value of these cumulants will be strongly sensitive to the short-range detail of the two-body interaction, and thus non-universal in nature.For this reason we do not assign physical significance to the values of n m k and χ k , but treat them as necessary quantities required to regularize the short-range behavior of the theory.A detailed analysis of the dynamics of n m k and the influence of higher order correlations neglected in the resonance triplet model will be the subject of future work.
The molecular pairing field κ m k is four-body in nature, and thus negligibly small over the time regimes we consider.This is especially true for broad resonances, where its amplitude is further suppressed by the short molecular lifetime.It is included in the resonance triplet model for the sake of completeness at the level of second order correlations, and because its simulation is trivial in terms of the computational cost, yet allows us to confirm the above expectations numerically as well.We find that, for all analysis presented in the main text, κ m k can be safely excluded without altering the results.

III. THREE-BODY EQUATIONS IN SHORT-RANGE OR VACUUM LIMITS
In this section, we show that the equations of motion as introduced in Sec.II converge to the three-body Schrödinger equation in the short-range or vacuum limits.The analogous convergence to the two-body Schrödinger equation in the resonance doublet model has been shown in Ref. [S1], and may be compared to the similar convergence of the BCS gap equation [S10].Subsequently we discuss the numerical solution to the vacuum equations, and show that the binding energy of the Efimov trimer can be recognized as a non-universal scale in our many-body dynamics.
A. Derivation of coupled-channel three-body Schrödinger equation Consider a general three-body state |Ψ⟩ with energy E, such that Ĥ |Ψ⟩ = E |Ψ⟩.In the two-channel scenario, this state has four components, Here |Ψ aaa ⟩ denotes the component with all three particles in the open channel and |Ψ acc ⟩ , |Ψ cac ⟩ , |Ψ cca ⟩ are the atom-dimer components which constitute the closed three-body channel.The open and closed channel states are coupled through the two-body interaction V P Q as introduced in Sec.I, such that the three-body Schrödinger equation can be written as, Ĥaaa,aaa Ĥaaa,acc Ĥaaa,cac Ĥaaa,cca Ĥacc,aaa Ĥacc,acc 0 0 Ĥcac,aaa 0 Ĥcac,cac 0 Ĥcca,aaa 0 0 Ĥcca,cca We now adopt a three-particle plane wave basis |k, q, q ′ ⟩, and calculate the projection Ψ a (k, q) = ⟨k, q, −k − q|Ψ aaa ⟩.
From equation (S24), the closed-channel components of the three-body state can be directly related to |Ψ aaa ⟩ as, We now insert a full set of eigenfunctions of H acc,acc .Since the closed channel only contains the single molecular state |ϕ⟩, we then obtain,  , such that we can suppress the particle indices going forward.We reintroduce the coupling constant g = √ 2β ⟨ϕ|ζ⟩, and obtain the coupled system, Here we have used parity symmetry, i.e.Ψ a (−q, −k) = Ψ a (q, k).Equation (S28) constitutes a set of two coupled (integral) equations that solve the three-body problem.We note that this set is formally equivalent to the equations derived in Ref. [S11], although we obtain them by a different route.

B. Equations of motion in short-range or vacuum limits
We now consider the coupled evolution of the three-body cumulants R a k,q and κ am k,q in the limit of large momenta (short range).Then we find the equations of motion, The above set is fully equivalent to Eq. (S28) upon adopting the ansatz, Hence the short-range limit of our cumulant model converges to the three-body Schrödinger equation as expected.Note that the set (S29) can also be obtained by taking the vacuum limit where all populations vanish.At longer range or in the presence of a medium the ansatz in Eq. (S30) corresponds to an approximation where the density effects are treated as quasi-stationary [S12-S14].

C. Efimov trimer spectrum
In principal, one can now solve Eq. (S28) to obtain the binding energies E (n) * of the Efimov trimers at unitarity.However, since the two-body problem with our separable potential admits an analytic solution, it is advantageous to instead formulate the three-body equations in terms of the two-body transition matrix.First, we generalize Eq. (S5) off the energy shell [S15], which we only have to evaluate for z < 0.Here E c is the threshold energy of the closed channel, and v c is the closedchannel interaction strength such that VQQ = 1 √ V v c |ζ⟩ ⟨ζ|.Both E c and v c are fixed by choosing the asymptotic energetic separation between the channels at the bare resonance position ν = 0.The Efimov spectrum is independent of this energy, provided that it is chosen large enough such that the closed-channel is far away from a shape resonance.
Then, in order to obtain the trimer binding energy E * < 0, we follow the Faddeev approach and decompose the three-body state |Ψ⟩ as [S16-S18], (0) * ā = 0.226 for the ground state trimer obtained with a single-channel interaction, valid for R * ≪ ā.Dash-dotted lines show the theoretical prediction from Ref. [S11], which holds in the narrow resonance limit R * ≫ ā.In (b) we show the Efimovian oscillation of the three-body contact density C3, which increases in frequency as the Feshbach resonance broadens in accordance with panel (a).To extract the exact frequency, we repeat the procedure in Ref. [S9] and fit the growth of C3 to a function The obtained value of b3 then represents the characteristic oscillation frequency, which we plot in panel (c) and compare directly with the ground state Efimov frequency ω which is an integral version of the three-body Schrödinger equation, with the boundary condition that the three-body wavefunction vanishes at infinity.Here Ĝ3B 0 (E) is the uncoupled three-body Green's function, T (E) is the two-body transition matrix generalised to the three-body space, and P = P+ + P− [S17].Equation (S33) may be projected onto momentum states, after which it can be reduced to a one-dimensional integral equation, where κ * = m|E * |/ℏ 2 .After transferring to the continuum limit the angular part of the integration can be done analytically, after which κ * is obtained numerically by solving the integral equation using Gaussian quadrature.The result is shown for a range of R * /ā in Fig. S1(a).In the narrow resonance limit R * ≫ ā, we also plot the analytical prediction κ * R * = 2.6531 mod e π/s0 obtained in Ref. [S11].We additionally compare the result with the ground state three-body parameter κ * ā = 0.226 obtained from a single-channel interaction, correspondent with the broad resonance limit R * ≪ ā.Both limits match well with our numerical results.

IV. GROWTH OF THE EXCITED STATE FRACTION
From Eq. (S16) we recognize that the atomic excitation density is sourced by two distinct scattering processes.To make this explicit, we write ℏ ∂ ∂t n a k = s k + s k , where, s Evidently the tail of s k is splits into a two-body and a three-body contribution.The diagrammatic representations of the associated processes are shown in Fig. 2 of the main text.As we show in Sec.V, the two-body contribution can be directly related to the two-body contact density C 2 .are shown with vertical lines.As the density of the gas decreases the resonance moves to higher momenta, signifying its finite-range nature.At later times the width of the resonance also increases since it scales with the magnitude of κ am k .
We now consider the second source s k , which is neglected in the resonance doublet model.Its diagrammatic representation is shown in Fig. S2(b).From the equation of motion (S19) for χ k we derive that s k will vanish in the formal zero-range limit Λ → ∞.For finite values of Λ, s k induces an asymmetric "resonance" feature in the tail of the momentum distribution, see Fig. S3.As we also mentioned in Sec.II, this resonance compensates for the additional energy introduced into the system via κ am k , ensuring that the total energy is always conserved.Evidently however, the actual shape of the peak has unphysical features, most notably a negative excitation density for momenta k > k res .As n m k is sourced purely by s k , its value is similarly unphysical in the current model.The negativity can be attributed to the truncation of the cumulant hierarchy, in which the open channel analogue of χ, ⟨a † a † a⟩, is neglected.Inclusion of this cumulant however, induces hierarchical couplings to higher order correlations that significantly complicate the model.For the scope of the present work, it is sufficient that the non-universal effects of s (2) k are negligible for the momenta of interest k/k n < 4 which is indeed the case as shown in Fig. S3 and indicates that our results in Fig. 4 of the main text are universal.In addition we note that the population densities n kres form a practically negligible fraction of the total population (∼ 10 −5 n).

FIG. 2 .
FIG. 2. Two-body (a) and three-body (b) scattering processes that drive the depletion of the atomic condensate.Atomic(molecular) states shown with single(double) lines, and condensed states shown in dashed red.Each vertex represents atom-molecule conversion.

1 √V
S26) Here |k, ϕ 23 ⟩ defines a product state in which a molecular state |ϕ 23 ⟩ composed of particles 2 and 3 moves as a plane wave with center-of-mass momentum −k.Meanwhile, atom 1 moves similarly as a plane wave with momentum k.We now define a momentum-space closed channel wave function Φ 23 (k) = ⟨k, ϕ 23 |Ψ acc ⟩ / √ 2, and express the coupling potential in the separable form H acc,aaa = |ζ 23 ⟩ β ⟨ζ 23 |, consistent with the conversion of atoms 2 and 3 into a closed-channel molecule.Repeating the above procedure for the other closed channel components |Ψ cac ⟩ , |Ψ cca ⟩ and projecting Eq. (S24) on the open-channel subspace, we finally obtain,

)FIG. S1 .
FIG. S1.Illustrations of Efimov physics embedded in the resonance triplet model.In (a) we plot the three-body binding wavenumber κ (n) * = m|E (n) * | /ℏ 2 as a function of R * for the three lowest lying Efimov trimer states n = 0, 1, 2, obtained from a numerical solution of Eq. (S34).The black dashed line shows the van der Waals universal three-body parameter κ * | /ℏ extracted from panel (a).The excellent match shows that the three-body Schrödinger equation is correctly included in the many-body dynamics.

ℏ 4 k 4 * m − 4mg 2 ℏ 2 k 2
Fig.S2.For now we focus on the first source s (1) k .To examine the tail of the momentum distribution, we expand the formal solution of the equation of motion (S17) for κ a k in powers of 1/k.Then, s Im iψ m ∂ ∂t ψ Im(ψ m κ am * k ψ a ).(S36)

FIG. S2 .
FIG. S2.Diagrammatic representation of s (1) k in (a) and s (2)k in (b).Atomic and molecular momenta shown with single and double lines respectively.Atomic(molecular) states shown with single(double) lines, and condensed states shown in dashed red.Since we take the imaginary part, each source consists of the balance between forward and backward processes.

FIG. S3 .
FIG. S3.Single-particle momentum distribution n a k for different values of the density scale nr 3 vdW , at t = tn in (a) and t = 2 tn in (b).Computed with resonance width R * kn = 0.4.The universal doublet curve is shown in grey, and the positions kres of the resonances due to s (2) k