Pair formation in quenched unitary Bose gases

TL;DR

This paper investigates how a quenched unitary Bose gas develops many-body-induced pairs, revealing a process similar to Cooper pairing, with results aligning with experiments for fast sweeps but highlighting the need for higher-order effects in other regimes.

Contribution

It introduces a many-body model including three-body losses to describe pair formation in quenched Bose gases, demonstrating a phenomenon analogous to Cooper pairing.

Findings

Formation of extended pairs due to many-body correlations.

Quantitative agreement with experiments for fast sweeps.

Deviations at slower sweeps suggest higher-order effects are important.

Abstract

We study a degenerate Bose gas quenched to unitarity by solving a many-body model including three-body losses and correlations up to second order. As the gas evolves in this strongly-interacting regime, the buildup of correlations leads to the formation of extended pairs bound purely by many-body effects, analogous to the phenomenon of Cooper pairing in the BCS regime of the Fermi gas. Through fast sweeps away from unitarity, we detail how the correlation growth and formation of bound pairs emerge in the fraction of unbound atoms remaining post sweep, finding quantitative agreement with experiment. We comment on the possible role of higher-order effects in explaining the deviation of our theoretical results from experiment for slower sweeps and longer times spent in the unitary regime.

Figures

Equations

$$k_n{a}_{\mathrm{eff}}=1.58+3.44{\left(\frac{t_n}{t_{\mathrm{hold}}}\right)}^2,$$

$$\hat{H}=\sum _{\mathbf{k}}\frac{{\hslash }^2{k}^2}{2m}{\hat{a}}_{\mathbf{k}}^{†}{\hat{a}}_{\mathbf{k}}+\sum _{\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{q}}{V}_{\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{q}}{\hat{a}}_{\mathbf{k}+\mathbf{q}}^{†}{\hat{a}}_{{\mathbf{k}}^{\prime }-\mathbf{q}}^{†}{\hat{a}}_{{\mathbf{k}}^{\prime }}{\hat{a}}_{\mathbf{k}},$$

$$i\hslash {\dot{\psi }}_0$$

$$\hslash {\dot{\rho }}_{\mathbf{k}}$$

$$i\hslash {\dot{\kappa }}_{\mathbf{k}}$$

$$h_{\mathbf{k}}=\frac{{\hslash }^2{k}^2}{2m}+2g\left(|\zeta (\mathbf{k}){|}^2|{\psi }_0{|}^2+\sum _{\mathbf{q}\ne 0}|\zeta (\mathbf{k}-\mathbf{q}){|}^2{\rho }_{\mathbf{q}}\right)$$

$${\Delta }_{\mathbf{k}}=g\zeta (2\mathbf{k})\left({\zeta }^{\ast }(0){\psi }_0^2+\sum _{\mathbf{q}\ne 0}{\zeta }^{\ast }(2\mathbf{q}){\kappa }_{\mathbf{q}}\right)$$

$${\hat{b}}_0^{†}\equiv \sum _{\mathbf{k}}\frac{{\varphi }_{\ast }(k)}{\sqrt{2}}{\hat{a}}_{-\mathbf{k}}^{†}{\hat{a}}_{\mathbf{k}}^{†},$$

$$\begin{array}{rl}\frac{2N_{\mathrm{mol}}}{N_{\mathrm{in}}}& =V\sum _{\mathbf{k}}|{\varphi }_{\ast }(k){|}^2\left(|{\psi }_0{|}^4{\delta }_{\mathbf{k},0}+\frac{2}{V}{\rho }_{\mathbf{k}}^2\right)\\ & +V{\left|\sum _{\mathbf{k}}{\varphi }_{\ast }(k){\kappa }_{\mathbf{k}}^{\ast }\right|}^2\\ & +V\sum _{\mathbf{k}}2\text{mbox}Re\left[{\varphi }_{\ast }(0)[{\psi }_0^{†}{]}^2{\kappa }_{\mathbf{k}}[{\varphi }_{\ast }(k){]}^{\ast }\right],\end{array}$$

$${\varphi }_{\ast }(k)=\frac{\sqrt{\mathcal{N}{a}_{\ast }^3}}{1+(k{a}_{\ast }{)}^2},$$

$${\left(a_{\ast }-a_{\mathrm{bg}}\right)}^2\left(a_{\ast }-\bar{a}\right)=\frac{\hslash \Delta B{a}_{\mathrm{bg}}}{2mR},$$

$$E_{2B}^{(\nu )}{\varphi }_{\nu }^{\mathrm{R}}(\mathbf{k})$$

$${\kappa }_{\mathbf{k}}=\sum _{\nu }{c}_{\nu }{\varphi }_{\nu }^{\mathrm{R}}(\mathbf{k})\iff c_{\nu }=\sum _{\mathbf{k}}[{\varphi }_{\nu }^{\mathrm{L}}(\mathbf{k}){]}^{\ast }{\kappa }_{\mathbf{k}},$$

$$\begin{split}\frac{2N_{\mathrm{D}}}{N_{\mathrm{in}}}&=V\left|\sum_{\mathbf{k}}\phi_{\ast}(k)[\phi_{\mathrm{D}}^{\mathrm{R}}(\mathbf{k})]^{\ast}\sum_{\mathbf{q}}[\kappa_{\mathbf{q}}]^{\ast}\phi_{\mathrm{D}}^{\mathrm{L}}(\mathbf{q})\right|^{2}\\ &+V\sum_{\mathbf{k}}\,2\mbox{Re}\Big{[}\phi_{*}(0)[\psi_{0}^{\dagger}]^{2}\left(\sum_{\mathbf{q}}[\phi_{\mathrm{D}}^{\mathrm{L}}(\mathbf{q})]^{\ast}\kappa_{\mathbf{q}}\right)\\ &\times\phi_{\mathrm{D}}^{\mathrm{R}}(\mathbf{k})[\phi_{*}(k)]^{\ast}\Big{]}.\end{split}$$

$$\frac{\dot{N}(t)}{N(t)}$$

$$K_3^{\mathrm{eff}}(n(t))$$

$$i\hslash {\dot{\psi }}_0$$

$$\hslash {\dot{\rho }}_{\mathbf{k}}$$

$$N_{\mathrm{free}}(t_{\mathrm{hold}},1/R)=N(t_{\mathrm{hold}})-2N_{\mathrm{mol}}(t_{\mathrm{hold}},1/R),$$

Full text

I Introduction

In ultracold quantum gases, precision control of magnetically tunable Feshbach resonances makes it possible to tune the effective interaction strength, characterized by the $s$-wave scattering length $a$ Chin et al. (2010). As $a$ becomes much larger than the interparticle spacing $n^{-1/3}$, where $n$ is the atomic density, the gas enters the unitary regime $(n|a|^{3}\gg 1)$. At unitarity ($|a|\to\infty$), interactions between atoms are as strong as allowed by quantum mechanics. Moreover, the macroscopic properties of unitary quantum gases appear insensitive to microscopic physics and therefore paradigmatic for other strongly correlated systems, including the inner crust of neutron stars and the quark-gluon plasma Schäfer and Teaney (2009); Braaten and Hammer (2003). The universality of the unitary Fermi gas has been both theoretically and experimentally well established over the past two decades Zwerger (2011). Under the universality hypothesis, the unitary Bose gas is also expected to behave similarly, with thermodynamic properties and relations that scale continuously solely with the “Fermi”scales constructed from powers of $n$, including the Fermi wave number $k_{n}=(6\pi^{2}n)^{1/3}$, energy $E_{n}=\hbar^{2}k_{n}^{2}/2m$, and time $t_{n}=\hbar/E_{n}$, where $m$ is the atomic mass Ho (2004).

Unlike their fermionic counterparts, at unitarity three bosons may form an infinite series of bound Efimov trimers Braaten and Hammer (2006) with characteristic finite size set by the three-body parameter $\kappa_{*}$  Wang et al. (2012); Naidon et al. (2014); Mestrom et al. (2019). Whereas Pauli repulsion suppresses three-body losses for fermions, the Efimov effect leads to a catastrophic $a^{4}$ scaling of three-body losses near unitarity, and therefore the unitary Bose gas is inherently unstable. In Refs. Makotyn et al. (2014); Klauss et al. (2017); Eigen et al. (2017, 2018), this barrier was overcome through a fast quench from the weakly interacting to the unitary regime, where the establishment of a steady state was observed before heating dominates. Time-resolved studies of the single-particle momentum distribution in Ref. Eigen et al. (2018) revealed that the theoretically predicted prethermal state Yin and Radzihovsky (2016); Sykes et al. (2014); Kain and Ling (2014) transitions to steady state prior to being overcome by heating. Although these findings, combined with studies of loss dynamics in Refs. Makotyn et al. (2014); Klauss et al. (2017); Eigen et al. (2017), are consistent with the universality hypothesis, a macroscopic population of Efimov trimers was observed in Ref. Klauss et al. (2017). Understanding the role of the Efimov effect Colussi et al. (2018a); D’Incao et al. (2018) and dynamics of higher-order correlations Kira (2015); Colussi et al. (2018b); Van Regemortel et al. (2018) in the quenched unitary Bose gas remains, however, an ongoing pursuit in the community.

The difficulties of probing the system at unitarity require that experiments return to the more stable and better-understood weakly interacting regime. During the course of the experiment, we have to distinguish different types of atomic pairs: (i) pairs of atoms with opposite momentum, analogous to Cooper pairs in Fermi gases, (ii) embedded dimers at unitarity whose size is determined by the mean interparticle separation, and (iii) weakly bound molecules away from unitarity, whose size is determined by the scattering length.

According to the experimental procedure of Refs. Makotyn et al. (2014); Klauss et al. (2017); Eigen et al. (2017, 2018), illustrated in Fig. 1, a Bose gas is initially quenched from the weakly interacting to the unitary regime, held there for a variable time $t_{\mathrm{hold}}$, and finally probed again in the weakly interacting regime. Here, the size of a molecule is much smaller than the mean interparticle separation, and the distinction between unbound and bound atoms is physically meaningful again Góral et al. (2004). In the unitary regime, unbound pairs progressively localize onto the scale of the interparticle spacing, purely due to many-body effects Colussi et al. (2018b). The nature of these embedded dimers is reflected by a universal time-dependent size $a_{\mathrm{eff}}$, fit to the universal form

[TABLE]

which indicates a transition from unbound ($a_{\mathrm{eff}}\to\infty$) to bound ($a_{\mathrm{eff}}\sim k_{n}^{-1}$) on Fermi timescales, as we will discuss in Sec. II.3. It is interesting to note the analogy of pair formation in the quenched unitary Bose gas to pair formation in the unitary Fermi gas Koetsier et al. (2009), which is at the center of the so-called BCS-BEC crossover. When entering this crossover from the Bardeen-Cooper-Schrieffer (BCS) side, fermionic pairs, loosely bound by the medium, smoothly evolve into tightly bound molecules that are stable even without the medium, when passing through to the Bose-Einstein condensation (BEC) side, while the effective atomic interaction changes from attractive to repulsive Zwerger (2011). For these experiments, a very successful technique was employed utilizing fast magnetic field sweeps to effectively project the fermionic pairs onto molecules throughout the whole crossover regime Regal et al. (2004); Hodby et al. (2005); Altman and Vishwanath (2005); Matyjaśkiewicz et al. (2008); Ketterle and Zwierlein (2008); Strinati et al. (2018).

In this work, we quench an initially pure Bose condensate to unitarity and track the resultant dynamics up to the level of two-body correlations, while including universal three-body losses phenomenologically. We then model the final step shown in Fig. 1 by a fast-sweep projection technique in the spirit of Ref. Altman and Vishwanath (2005), count the number of remaining unbound atoms, and compare quantitatively our results with the experimental findings of Ref. Eigen et al. (2017). Unlike in the experiment, in our model, we are able to distinguish between three-body losses and formation of molecules when determining the number of remaining unbound atoms. Through this ability, we estimate the universal three-body loss-rate coefficient by refitting the experimental data of Ref. Eigen et al. (2017). We also compare the predictions of our model for the number of unbound atoms with the results of that work, finding generally good agreement for fast ramp rates and for slower ramp rates at earlier times ($t_{\mathrm{hold}}\lesssim 0.5t_{n}$). As correlations grow and the condensate becomes increasingly depleted for longer times spent in the unitary regime, we highlight the dominant contribution of the embedded dimers in the number of unbound atoms detected after fast-sweep projection away from unitarity.

The organization of this work is as follows. In Sec. II, we outline our many-body model (Sec. II.1), adapt the technique of fast-sweep projection from Ref. Altman and Vishwanath (2005) (Sec. II.2) for Bose gases, and develop the theory of bound pairs in the unitary regime discussed in Ref. Colussi et al. (2018b) (Sec. II.3). In Sec. III, three-body losses are introduced phenomenologically into our many-body model, and in Sec. IV, we discuss the results of our model and compare them with the experimental findings of Ref. Eigen et al. (2017). We conclude in Sec. V and comment on prospects for future study.

II Model

II.1 Many-body equations

We model a uniform gas of identical spinless bosons interacting via pairwise interactions described by the single-channel many-body Hamiltonian

[TABLE]

where $V_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}=(g/2)\zeta(\mathbf{k}-\mathbf{k}^{\prime}+2\mathbf{q})\zeta^{\ast}(\mathbf{k}-\mathbf{k}^{\prime})$ is a non local separable potential with interaction strength $g$, step-function form factor $\zeta(\mathbf{k})=\theta(\Lambda-|\mathbf{k}|/2)$, and finite cutoff $\Lambda$, giving rise to a finite-range interaction both in momentum and position space. This model is suitable for describing open-channel dominated Feshbach resonances, which includes all degenerate unitary Bose gas experiments to date Makotyn et al. (2014); Fletcher et al. (2017); Klauss et al. (2017); Eigen et al. (2017, 2018). To fix the free parameters of the separable potential, we first set the strength of the potential, $g=U_{0}\Gamma$, where $U_{0}=4\pi\hbar^{2}a/m$ and $\Gamma=(1-2a\Lambda/\pi)^{-1}$, to reproduce the exact two-body $T$ matrix in the zero-energy limit Colussi et al. (2018b); Kokkelmans et al. (2002). To fix $\Lambda$, we follow Ref. Colussi et al. (2018b) and set $\Lambda=2/\pi\bar{a}$ to obtain finite-range corrections to the binding energy of the Feshbach molecule $E_{\mathrm{b}}\simeq-\hbar^{2}/m(a-\bar{a})^{2}$, valid only to first order in $1/\Lambda a$, and where $\bar{a}=0.955r_{\mathrm{vdW}}$ is the mean scattering length that depends on the van der Waals length $r_{\mathrm{vdW}}$, for a particular atomic species Chin et al. (2010). Consequently, at unitarity we obtain a finite interaction strength $g=-\pi^{3}\hbar^{2}\bar{a}/m$ for $a\to\infty$.

To model the condensate and excitations, we make the Bogoliubov approximation Dalfovo et al. (1999) and decompose the operator $\hat{a}_{\mathbf{k}}=\psi_{\mathbf{k}}+\delta\hat{a}_{\mathbf{k}}$ with $\braket{\delta\hat{a}_{\mathbf{k}}}=0$. We assume that only the atomic condensate is macroscopically occupied so that $\braket{\hat{a}_{\mathbf{k}}}=\psi_{0}\delta_{\mathbf{k},0}$ and consider only fluctuations of the excitations. These assumptions are valid provided the excited modes are not macroscopically occupied. Furthermore, we build our many-body theory from the cumulant expansion Kira (2015); Colussi et al. (2018b), which separates clusters of correlated particles within an interacting many-body system. Here, we consider only up to second-order clusters (correlations), described by the condensate wave function $\psi_{0}$ and the one-body $\rho_{\mathbf{k}}\equiv\braket{\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{\mathbf{k}}}$ and pairing $\kappa_{\mathbf{k}}\equiv\braket{\hat{a}_{-\mathbf{k}}\hat{a}_{\mathbf{k}}}$ density matrices for excitations Blaizot and Ripka (1986).

We derive the Hartree-Fock Bogoliubov (HFB) equations Blaizot and Ripka (1986) from the Heisenberg equation of motion for one or two $\hat{a}$-operator products and evaluate the expectation values in cumulant expansion, neglecting clusters of three or more particles. Summarily, if $\hat{\mathcal{O}}$ is a specific operator, using $i\hbar\braket{d\hat{\mathcal{O}}/dt}=\braket{[\hat{\mathcal{O}},\hat{H}]}$ one has

[TABLE]

where

[TABLE]

and

[TABLE]

are the Hartree-Fock Hamiltonian and the pairing field, respectively Blaizot and Ripka (1986). The HFB theory results in a mean-field description, typically suitable for the weakly interacting regime where $n|a|^{3}\ll 1$. However, here we formulate a finite-range HFB theory, which yields a finite mean-field energy $g$ at unitarity, as discussed above. We argue that this theory is applicable for strong interactions since $g$ remains small respect to $E_{n}$, or equivalently $nr_{\mathrm{vdW}}^{3}\ll 1$. Note that this condition is well satisfied for all experiments in the unitary regime to date ($nr_{\mathrm{vdW}}^{3}<10^{-5}$) Holland et al. (2001); Colussi et al. (2018b); Makotyn et al. (2014); Klauss et al. (2017); Eigen et al. (2017, 2018).

To simulate the first two steps of the experimental sequence illustrated in Fig. 1, Eqs. (3)-(5) are solved at fixed initial density $n_{\mathrm{in}}=N_{\mathrm{in}}/V$, where $N_{\mathrm{in}}\equiv N(t=0)$ with total atom number $N(t)$ in a volume $V$ 111The Fermi scales $k_{n}$, $E_{n}$, and $t_{n}$ are defined in terms of the initial density $n_{\mathrm{in}}$.. In particular, we consider experiments done in a box-trap, modelled as a uniform system Gaunt et al. (2013). We begin at $t=0$ from a pure condensate with $|\psi_{0}|^{2}=n_{\mathrm{in}}$. The scattering length is then ramped over $2$ $\mu$s to unitarity, where the system evolves for a varying amount of time, $t_{\mathrm{hold}}$. As the gas evolves at unitarity and in the absence of losses, the condensate fraction becomes depleted as correlated pair excitations are generated and counted by $\rho_{\mathbf{k}}$ as studied in Ref. Sykes et al. (2014). We expect that the increase of $\rho_{\mathbf{k}}$ beyond unity makes higher-order cumulants strongly driven and their inclusion in the model cannot be justified. Therefore, following Ref. Colussi et al. (2018b), we restrict our analysis to $t\leq 2t_{n}$ where $\rho_{\mathbf{k}}<1$ remains valid.

II.2 Fast-sweep projection away from unitarity

We finally model the third step of Fig. 1 with a projection of the many-body state at unitarity onto a molecular state at finite scattering length and count the number of molecules. Intuitively, in the limiting case of a sudden switch of the magnetic field, the number of molecules may be calculated, to good approximation, by simply projecting the state at unitarity onto molecules at the final magnetic field $B_{\mathrm{end}}$. For finite ramp rates $R=-dB/dt$, this approximation is not valid. In this case, the number of molecules may be calculated approximately by projection onto an effective molecular state $\phi_{\ast}$ with scattering length $a_{\ast}$ larger than the final scattering length $a_{\mathrm{end}}$ and intermediate to both the sudden and adiabatic cases, as detailed in Ref. Altman and Vishwanath (2005). This method provides an indirect measure of the buildup of correlations at unitarity. The conceptual problem of bound pairs in the unitary regime is revisited in Sec. II.3.

We construct a compound bosonic operator

[TABLE]

counting molecules away from unitarity with zero center of mass and relative momentum $\mathbf{k}$ of the constituent atoms, where $\phi_{\ast}(k)$ is a molecular wave function with a finite scattering length $a_{\ast}$ whose value will be specified shortly. By construction, the $\hat{b}$ operator satisfies $[\hat{b}_{0},\hat{b}_{0}]=[\hat{b}_{0}^{\dagger},\hat{b}_{0}^{\dagger}]=0$, and the canonical commutation relation $[\hat{b}_{0},\hat{b}_{0}^{\dagger}]=1+\sum_{\mathbf{k}}|\phi_{\ast}(k)|^{2}(\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}_{\mathbf{k}}+\hat{a}^{\dagger}_{-\mathbf{k}}\hat{a}_{-\mathbf{k}})$ is approximately well satisfied $\braket{[\hat{b}_{0},\hat{b}_{0}^{\dagger}]}\simeq 1$ away from unitarity, where the molecules are spatially much smaller than the interparticle spacing. We note also that the approach of counting composite bosons [Eq. (8)] has been also used extensively for counting fermionic pairs along the BEC-BCS crossover Nozières and Schmitt-Rink (1985); Altman and Vishwanath (2005); Perali et al. (2005).

We evaluate the expectation value $\braket{\hat{b}^{\dagger}_{0}\hat{b}_{0}}$, which can be expanded in terms of first- and second-order cumulants, consistently with the theory presented in Sec. II.1. Therefore, the molecular fraction is

[TABLE]

where $N_{\mathrm{in}}/2$ is the total possible number of molecules. At $t_{\mathrm{hold}}=0$ immediately following the completion of the quench, $|\psi_{0}|^{2}\approx n$ and $\rho_{\mathbf{k}}\approx\kappa_{\mathbf{k}}\approx 0$, and therefore only the first term on the right-hand side of Eq. (9) contributes. This contribution can be interpreted as the overlap of the molecular wave function and the atomic mean field Góral et al. (2004) and scales as $na_{\ast}^{3}$ proportional to the ratio of atomic and molecular volumes. This overlap must be insignificant so that $na_{\ast}^{3}<1$, and molecules can be separated from the many-body background. The remaining terms in Eq. (9) measure the overlap between molecular and pairing wave functions Avdeenkov and Bohn (2005) and reflect the development of correlations as the gas evolves in the unitary regime. We note that Eq. (9) is in agreement with the first-quantized multichannel description in position space found in Ref. Góral et al. (2004).

In the evaluation of Eq. (9), the molecular wave function has the universal form

[TABLE]

valid provided $a_{\ast}\gg r_{\mathrm{vdW}}$ Chin et al. (2010). The normalization constant $\mathcal{N}=4\pi^{2}/\{\arctan(\Lambda a_{\ast})-\Lambda a_{\ast}/[1+(\Lambda a_{\ast})^{2}]\}$ ensures that $\sum_{\mathbf{k}}^{\Lambda}|\phi_{\ast}(k)|^{2}=1$. As mentioned above, Eq. (10) is written in terms of $a_{\ast}$, whose value depends on $R$. $a_{\ast}$ represents the point at which the evolution of the system under the ramp changes from sudden to adiabatic, and the creation and dissociation of molecules is halted Altman and Vishwanath (2005); Matyjaśkiewicz et al. (2008). Quantitatively, this occurs when $E_{\mathrm{b}}/\hbar=E_{\mathrm{b}}^{-1}\dot{E_{\mathrm{b}}}$ is satisfied, where $E_{\mathrm{b}}=-\hbar^{2}/m(a-\bar{a})^{2}$ is the molecular binding energy including finite-range effects Chin et al. (2010). We obtain specific values of $a_{\ast}$ from the real solution of the third-order polynomial equation

[TABLE]

where $a_{\mathrm{bg}}$ is the background scattering length, and $\Delta B$ is the width of the Feshbach resonance Chin et al. (2010)222To derive Eq. (11), the two-channel expression $a=a_{\mathrm{bg}}-mG^{2}/(8\pi\hbar^{2}\nu)$ was used, where $G$ is the coupling between the two channels, $\nu=\Delta\mu(B-B_{0})$ is the detuning between collision energy and bound-state energy in the closed channel, and $\Delta\mu$ is the difference in magnetic moments between the two channels Duine and Stoof (2004); Kokkelmans and Holland (2002). This introduces an explicit dependence on $1/R$.. In the sudden limit, the initial state is projected onto the final scattering length $a_{\mathrm{end}}\equiv a(B_{\mathrm{end}})\gg a_{\ast}$ which is on the order of $a_{\mathrm{bg}}$ or $\bar{a}$ for $1/R\to 0$.

The dependence of $a_{\ast}$ on the ramp rate is shown by the solid red line in Fig. 2. Generally, larger values of $a_{\ast}^{-1}$ indicate a faster ramp and the many-body state at unitarity is projected onto more localized molecules. Consequently, $\phi_{\ast}(k)$ will be less pronounced at low momenta than for slower ramps, which can be seen in the insets of Fig. 2.

II.3 Embedded dimers at unitarity

To link the buildup of correlations at unitarity with the fast-sweep production of molecules, it is instructive to introduce a many-body length scale that can be compared with $a_{\ast}$. Here, we follow the approach outlined in Ref. Colussi et al. (2018b) and study embedded two-body bound states at unitarity.

To obtain the spectrum of these dimers embedded in the unitary Bose gas, the homogeneous part of Eq. (5), including only terms dependent on $\kappa_{\mathbf{k}}$, is solved as a two-body Schrödinger equation in the quasi-stationary limit Colussi et al. (2018b); Kira (2015). This approach is valid provided $\kappa_{\mathbf{k}}$ evolves faster than the density dynamics and scattering among clusters (inhomogeneous terms) are neglected, in which case one obtains an eigenvalue equation

[TABLE]

where $E_{2B}^{(\nu)}$ is a two-body eigenenergy and $\phi_{\nu}^{\mathrm{R}}(\mathbf{k})$ is a right-handed wave function  Kira (2015); Colussi et al. (2018b). The left-handed wave function $\phi_{\nu}^{\mathrm{L}}(\mathbf{k})$ is related via $\phi_{\nu}^{\mathrm{R}}(\mathbf{k})=(1+2\rho_{\mathbf{k}})\phi_{\nu}^{\mathrm{L}}(\mathbf{k})$, and they satisfy the usual orthogonality $\sum_{\mathbf{k}}[\phi_{\nu}^{\mathrm{L}}(\mathbf{k})]^{\ast}\phi_{\mu}^{\mathrm{R}}(\mathbf{k})=\delta_{\nu,\mu}$ and normalization $\sum_{\nu}[\phi_{\nu}^{\mathrm{L}}(\mathbf{k})]^{\ast}\phi_{\nu}^{\mathrm{R}}(\mathbf{q})=\delta_{\mathbf{k},\mathbf{q}}$ conditions.

It is illustrative to compare Eq. (12) with the Schrödinger equation for a Cooper pair in the BEC-BCS crossover, which depends instead on the Pauli-blocking factor $(1-2\rho_{\mathbf{k}})$  Nozières and Schmitt-Rink (1985); Strinati et al. (2018). Whereas the blocking factor in the BEC-BCS crossover theory forbids scattering at occupied intermediate states Fetter and Walecka (2003), the intermediate states for a Bose gas are Bose enhanced Proukakis et al. (1998). Both effects may lead to weakly bound pairs which are held together purely by many-body effects, whose presence was predicted in the finite-temperature phase diagram of the strongly interacting Bose gas Koetsier et al. (2009).

Following Ref. Colussi et al. (2018b), we track the gradual development in time of these embedded dimers, solutions of Eq. (12) with wave function $\phi_{\mathrm{D}}(k)$ and binding energy $E_{2B}^{\mathrm{D}}\equiv-\hbar^{2}/ma_{\mathrm{eff}}^{2}$. Equation (12) can be solved numerically, yielding $E_{\mathrm{2B}}^{D}$ as a function of time. The evolution of $E_{2B}^{\mathrm{D}}$ was fit in Ref. Colussi et al. (2018b), and we quote that result in Eq. (1). Initially, these dimers are basically unbound ($a_{\mathrm{eff}}\sim\infty$), but through the subsequent buildup of correlations and quantum depletion they are localized ($a_{\mathrm{eff}}\propto k_{n}^{-1}$) onto the Fermi scale and behave universally.

Comparing $a_{\ast}$ with $a_{\mathrm{eff}}$ provides a convenient way of characterizing the underlying physics of the fast-sweep projection. These scales are shown in Fig. 2, where the development of $a_{\mathrm{eff}}^{-1}$ as the gas evolves in the unitary regime is represented by the progression of horizontal lines in the shaded region. As discussed in Sec. II.2, the fast-sweep projection must be such that $k_{n}a_{\ast}\ll 1$ and therefore outside of the shaded region. These length scales may also be used to understand how the buildup of correlations influences the number of remaining unbound atoms after the fast-sweep projection. The evolution of $\phi_{\mathrm{D}}(k)$ with $t_{\mathrm{hold}}$ is shown along with $\phi_{\ast}(k)$ in Fig. 2 for two different ramp rates. The gradual localization of $\phi_{\mathrm{D}}(k)$ onto the Fermi scale leads to increasing overlap with $\phi_{\ast}(k)$. This behavior is more pronounced for slower ramps and for longer $t_{\mathrm{hold}}$. Therefore, we intuitively expect that embedded dimers make an increasing contribution to the overlap term in Eq. (9) and therefore the number of molecules produced by the fast-sweep projection.

To determine the role of the embedded dimers at unitarity, we decompose $\kappa_{\mathbf{k}}$ in the basis of $\phi_{\nu}^{\mathrm{R}}(\mathbf{k})$ as

[TABLE]

where the coefficient $c_{\nu}$ quantifies the relative weight of the component $\nu$ within the total $\kappa_{\mathbf{k}}$. We define the embedded dimer contribution $N_{\mathrm{D}}$ in Eq. (9), by evaluating only the component $\nu=D$ of Eq. (13)

[TABLE]

Figure 3 shows the ratio between embedded dimers and total number of molecules as a function of $(k_{n}a_{\mathrm{eff}})^{-1}$ after fast-sweep projections for three ramp rates of experimental interest. We find that by $t_{\mathrm{hold}}\sim 2t_{n}$, when $(k_{n}a_{\mathrm{eff}})^{-1}\sim 0.4$, embedded dimers make up $\approx 60\%$ of the detected molecules. Therefore, the fast-sweep projection increasingly converts embedded dimers into weakly bound molecules away from resonance, as the gas spends more time at unitarity, agreeing with the intuitive overlap picture shown in Fig. 2. We note that the behavior shown in Fig. 3 is reminiscent of the monotonic conversion of fermions pairs into molecules along the BEC-BCS crossover as a function of the scattering length Salasnich et al. (2005); Strinati et al. (2018).

III Modeling three-body losses

The development of strong correlations at unitarity is also accompanied by strong losses Makotyn et al. (2014); Fletcher et al. (2017); Klauss et al. (2017); Eigen et al. (2017, 2018); D’Incao et al. (2018); Sykes et al. (2014). In Ref. Eigen et al. (2018), by focusing on the early-time dynamics of the tail of the single-particle momentum distribution for $k/k_{n}\gtrsim 0.8$, it was possible to experimentally distinguish between the formation of a steady-state and long-time heating. However, this separation was not possible experimentally in Refs. Makotyn et al. (2014); Eigen et al. (2017) for observables depending on the full range of momentum. In the present work, we model the findings of Ref. Eigen et al. (2017) and study the number of unbound atoms detected following the completion of the experimental sequence illustrated in Fig. 1. In particular, the number of unbound atoms decreases in time because of two main phenomena, which are difficult to distinguish experimentally: molecular formation and three-body losses. Therefore, the inclusion of losses is required to make a quantitative comparison.

We assume a universal form for three-body losses scaling as $n^{2/3}$ in degenerate Bose gases, ignoring possible log-periodicities due to the Efimov effect D’Incao et al. (2018); Eismann et al. (2016)

[TABLE]

This gives an effective three-body loss coefficient

[TABLE]

that satisfies the standard relation $\dot{N}(t)/N(t)=-K_{3}^{\mathrm{eff}}N(t)^{2}/V^{2}$ for a uniform system Chin et al. (2010). We treat the constant $A$ as a free parameter that is varied in Sec. IV in order to fit experimental data of Ref. Eigen et al. (2017). The form of Eq. (15) was found experimentally in Refs. Eigen et al. (2017); Klauss et al. (2017) and is theoretically motivated by the universal substitution $a^{4}\to a_{\mathrm{eff}}^{4}$ in the scaling law of $K_{3}^{\mathrm{eff}}$ for shallow dimers as was suggested in Refs. Colussi et al. (2018a); Stoof and van Heugten (2014). For clarity, the universal scaling $n^{2/3}$ was found to be valid only for $t_{\mathrm{hold}}\lesssim 4t_{n}$, for longer time the loss rate scales as $n^{26/9}$, consistent with the results for a thermal gas at unitarity Petrov and Werner (2015); Rem et al. (2013) and beyond the limit of validity of our model, as discussed in Sec. II.1.

To incorporate three-body losses into the HFB equations, we consider the time derivative of the atomic density $dn(t)/dt=d(N(t)/V)/dt=d(|\psi_{0}(t)|^{2})/dt+\sum_{\mathbf{k}\neq 0}d\rho_{\mathbf{k}}(t)/dt$ and in order to satisfy Eq. (15) we only have to modify Eqs. (3) and (4) with additional terms

[TABLE]

where $\cdots$ represents the lossless terms of the HFB equations. We note that a similar phenomenological approach has been used at the level of the Gross-Pitaevksii equation in Ref. Ancilotto et al. (2015) and also to describe the Bosenova in Refs. Altin et al. (2011); Snyder et al. (2012). These approaches, however, did not include density dependence in the three-body loss coefficient. We also note that it should be possible to go beyond this phenomenological approach through a proper inclusion of third-order correlations into an extension of the many-body model outlined in Sec. II. These matters remain the subject of future study.

IV Results

In this section, we compare the results of our model to the experimental data in Ref. Eigen et al. (2017). We approximate the box cylindrical trap used in that work as a homogeneous gas Gaunt et al. (2013) and numerically solve the HFB equations including losses [Eqs. (5), (17), and  (18)] for the 39K Feshbach resonance at $B=402$ G with $a_{\mathrm{bg}}=-29a_{0}$, $\Delta B=-52$ G, and $\bar{a}=61.7a_{0}$ Chin et al. (2010). To mimic the experimental setup, we fix the initial density $n_{\mathrm{in}}$ and simulate up to $t_{\mathrm{hold}}=2t_{n}$, which is the range of validity of our model as discussed in Sec. II. We then calculate the total number of unbound atoms after the fast sweep away from unitarity from Eq. (9) for the ramp rates used experimentally. We calculate the number of free (unbound) atoms as

[TABLE]

where the ramp-rate dependence is indicated explicitly.

Before discussing the results, we comment on the validity of our approach. For the 39K Feshbach resonance at $B=402$ G, we find that for the ramp rates and initial densities considered $2.0\leq(k_{n}a_{\ast})^{-1}\leq 6.7$, and therefore the fast-sweep projection method outlined in Sec. II.2 can be applied. Although not analyzed in this work, we estimate that this method can also be applied to model the fast-sweep projection studied in Ref. Klauss et al. (2017) with 85Rb 333For the 85Rb Feshbach resonance at $B=155$ G studied in Ref. Klauss et al. (2017) with $a_{\mathrm{bg}}=-443a_{0}$, $\Delta B=10.7$ G, and $\bar{a}=79.1a_{0}$ Chin et al. (2010), we find that for initial densities between $0.2$ and $5.8\times 10^{12}\,\mbox{cm}^{-3}$ the range is $1.6\leq(k_{n}a_{\ast})^{-1}\leq 4.7$ for $1/R=12.5\mu s/G$.. For smaller $1/R$ and hence smaller $a_{\ast}$, we follow in the spirit Ref. Matyjaśkiewicz et al. (2008) and check the expression of $E_{\mathrm{b}}$ used to calculate Eq. (11) against a coupled-channel calculation tho , finding discrepancies of less than $5\%$.

Our results for $N_{\mathrm{free}}$ are compared against the experimental findings of Ref. Eigen et al. (2017) as a function of $t_{\mathrm{hold}}$ for initial density $n_{\mathrm{in}}=2.7\times 10^{12}$ cm*-3* and ramp rates $0.3$ and $6$ $\mu$s/G, as shown in Fig. 4. At $t_{\mathrm{hold}}=0$, the small gap between the theoretical results for the two different ramp rates is due solely to the first term on the right-hand side of Eq. (9) which scales as $n_{\mathrm{in}}a_{\ast}^{3}$ and therefore varies with the ramp rate [see Sec. II.2]. At later times, pair correlations begin to develop, and the overlap between embedded dimers at unitarity $[\phi_{\mathrm{D}}(k)]$ with molecules away from resonance $[\phi_{\ast}(k)]$ increases, as illustrated in Fig. 2. Consequently, the decrease of $N_{\mathrm{free}}$ shown in Fig. 5 is due jointly to molecular formation and three-body losses.

The constant $A$ was estimated in Ref. Eigen et al. (2017) for 39K as $A=0.28(3)$ by assuming that the $0.3$ $\mu$s/G ramp projects the gas at unitarity only onto unbound atoms [solid green line in Fig. 4]. In our model, we separate the contributions of molecular formation and loss, and it is therefore possible to test this assumption and provide an independent estimation of $A$ using the approach outlined in Sec. III. We therefore adjust $A$ in the HFB equations including losses [Eqs. (5), (17), and  (18)], and refit the $0.3$ $\mu$s/G experimental data as shown in Fig. 4. For this specific ramp, we find a molecular fraction $\approx 10\%$, which is compatible with the experimental estimate in Ref. Eigen et al. (2017). By comparing three values $A=\{0.28,0.20,0.18\}$ to the $0.3$ $\mu$s/G experimental data, we find that $A=0.20$ provides the best fit of the experimental results over the full range of $t_{\mathrm{hold}}$ considered in this work. For the slower $6$ $\mu$s/G ramp, we find that $A=0.20$ gives good agreement at early times until roughly $t_{\mathrm{hold}}\gtrsim 0.5t_{n}$. We discuss possible sources of this discrepancy at longer $t_{\mathrm{hold}}$ at the conclusion of this section.

Our results for $N_{\mathrm{free}}$ over a range of $1/R$ are compared against the experimental findings in Ref. Eigen et al. (2017) as shown in Fig. 5. The results shown in Fig. 5 are at fixed $t_{\mathrm{hold}}=1.9t_{n}$, nearing the limit of validity of our model [see Sec. II.1]. The intuitive picture, discussed in Secs. II.2, II.3 and illustrated in Fig. 2 provides a way to understand our results particularly at this later time where the bound pairs at unitarity play a dominant role (see Fig. 3). For smaller ramp rates, the largest values of $N_{\mathrm{free}}$ shown in Fig. 5 result from the fast-sweep projection occurring further away from unitarity where the overlap between embedded dimers $[\phi_{\mathrm{D}}(k)]$ with molecules $[\phi_{\ast}(k)]$ becomes minimal. We find good agreement with experiment only for the fastest ramps considered using the refitted value $A=0.20$. In Ref. Eigen et al. (2017) the ramp-rate dependence of $N_{\mathrm{free}}$ is fit to a Landau-Zener exponential $N_{\mathrm{free}}=\alpha+\beta\ \mbox{exp}(-\gamma/R)$ Hodby et al. (2005), where they found $1/\gamma=2.2(3)$ $\mu$s/G. From fitting the $A=0.20$ theoretical data in Fig. 5 (dot-dashed purple curve), we find $1/\gamma=4.1$ $\mu$s/G. The possible sources of discrepancy for slower ramps will be discussed at the end of this section.

In addition, we analyze $N_{\mathrm{free}}$ over a range of initial densities, $n_{\mathrm{in}}$, and compare against the experimental results in Ref. Eigen et al. (2017). Taking the refitted value $A=0.20$, we follow experiment and vary $n_{\mathrm{in}}$ between $1.3\times 10^{12}$ and $4.0\times 10^{12}$ cm*-3*, measuring the difference $\Delta N(1/R)\equiv N_{\mathrm{free}}(t_{\mathrm{hold}},0.3\,\mu\mathrm{s/G})-N_{\mathrm{free}}(t_{\mathrm{hold}},1/R)$ as shown in Figs. 6(a), and 6(b). We note that for fixed $A$ this is equivalent in our model to the difference $\Delta N(1/R)=2(N_{\mathrm{mol}}(t_{\mathrm{hold}},1/R)-N_{\mathrm{mol}}(t_{\mathrm{hold}},0.3\,\mu\mathrm{s/G}))$. At $t_{\mathrm{hold}}=0$, $\Delta N$ is nonzero due to the first term of Eq. (9) scaling as $n_{\mathrm{in}}a_{\ast}^{3}$ that was also discussed earlier in connection with Fig. 4. At later times, the gradual separation of the $\Delta N$ curves shown in Fig. 6(a) can be understood by comparing the density-dependent and independent length scales $a_{\mathrm{eff}}$ and $a_{\ast}$, respectively. The many-body length scale $a_{\mathrm{eff}}\propto n_{\mathrm{in}}^{-1/3}$ is sensitive to changes in the initial density, whereas $a_{\ast}$ remains fixed by the ramp rate $1/R$. Consequently, the overlap between $\phi_{\mathrm{D}}(k)$ and $\phi_{\ast}(k)$ increases with $n_{\mathrm{in}}$, which results in the separation of the theoretical $\Delta N$ curves in Fig. 6(a), where $1/R=6\,\mu$s/G. In Fig. 6(b), we also compare our results for $\Delta N$ at fixed $n_{\mathrm{in}}$ for ramp rates $3$ and $6$ $\mu$s/G, in order to differentiate between $1/R$ and $n_{\mathrm{in}}$ dependencies. As before, we attribute the separation of the theoretical $\Delta N$ curves to the time dependence of the overlap between $\phi_{\mathrm{D}}(k)$ and $\phi_{\ast}(k)$ and the dominance of the bound pairs at unitarity at later times [see Fig. 3]. This separation is reflected also in the experimental data shown in Figs. 6(a) and 6(b). In general, our predictions in Figs. 6(a) and (b) match the experimental data well until we begin to underestimate $\Delta N$ compared to experiment at times $t_{\mathrm{hold}}\gtrsim 0.5t_{n}$.

We now address the deviation between our theoretical predictions presented in this section and the experimental results of Ref. Eigen et al. (2017) for the $3$ and $6$ $\mu$s/G ramps over longer timescales $t_{\mathrm{hold}}\gtrsim 0.5t_{n}$. In Ref. Eigen et al. (2017), it was experimentally observed that a degenerate Bose gas quenched to the unitary regime undergoes a universal crossover to the thermal regime by $t_{\mathrm{hold}}/t_{n}\approx 4.0$. In the thermal regime, the three-body loss rate $\dot{N}/N$ scales as $n^{26/9}$ Rem et al. (2013). However, the quantitative agreement between theory and experiment for the loss-dominated $0.3$ $\mu$s/G ramp shown in Fig. 4 is consistent with the $2/3$ power law in the degenerate regime [see Eq. (15)].

In Ref. Eigen et al. (2018), it was experimentally observed that momentum modes with $k/k_{n}\gtrsim 0.8$ reach a prethermal steady state and plateau by $t_{\mathrm{hold}}\sim t_{n}$ before long-time heating dominates. In our model the momentum modes described by $\rho_{\mathbf{k}}$ in the HFB equations (see Sec. II A) do not plateau as function of $t_{\mathrm{hold}}$ but oscillate in time, as in Refs. Gao et al. ; Muñoz de las Heras et al. (2019); Sykes et al. (2014), where the dynamics at unitarity is described through a time-dependent coherent-state pairing wave function ansatz equivalent to the HFB model Rançon and Levin (2014); Corson and Bohn (2015). This would be most apparent for the slowest $6$ $\mu$s/G ramp [see Fig. 6(a)] where $2.5\leq(k_{n}a_{\ast})^{-1}\leq 3.2$, and therefore it is possible that the physics behind the plateau are responsible for the deviation between theory and experiment.

Finally, from the experimental findings in Ref. Klauss et al. (2017), a macroscopic population of Efimov trimers, corresponding to $8\%$ of the initial state, was found after performing a fast-sweep projection away from unitarity. To estimate the potential relevance of Efimov trimers, we follow Refs. Colussi et al. (2018a); D’Incao et al. (2018) and compare the Fermi scale with the size of the nearby first-excited trimer $R_{3b}^{(1)}=(1+s_{0}^{2})^{1/2}e^{\pi/s_{0}}/(3/2)^{1/2}\kappa_{*}$, where $s_{0}\approx 1.00624$ and $\kappa_{*}=0.226/r_{\mathrm{vdW}}$ is the universal three-body parameter Braaten and Hammer (2006); Wang et al. (2012); Naidon et al. (2014). For the density range considered in Fig. 6, we estimate that $1.7\leq k_{n}R^{(1)}_{\mathrm{3b}}\leq 2.5$. Based on the qualitative findings in Ref. D’Incao et al. (2018), the first-excited Efimov trimer population is expected to grow more slowly than the dimer contribution to the molecular fraction, and this may be partially responsible for the deviation at later times 444Here, we reference specifically Figs. S3(b) and S3(c) in Ref. D’Incao et al. (2018).. However, in that work a breakdown of the Landau-Zener behavior was found for increasing $t_{\mathrm{hold}}$, which qualitatively disagrees with the experimental and theoretical results shown in Fig. 5 displaying this behavior. We leave, however, the possibility of resolving this deviation by either including into our many-body model three-body correlations or equilibrating collisions Van Regemortel et al. (2018) as inspiration for future work.

V Conclusion

In this work, we present a dynamical model of the degenerate Bose gas quenched to unitarity, which we compare against recent experimental results Eigen et al. (2017) for the number of unbound atoms remaining after a fast-sweep ramp away from the unitary regime. We adopt the method of Ref. Altman and Vishwanath (2005) from the study of Cooper pairs in the BEC-BCS crossover and project the many-body state in the unitary regime onto molecular states away from unitarity. As the Bose gas evolves in the unitary regime, the buildup of correlations and quantum depletion leads to the formation of pairs bound purely by many-body effects as studied in Ref. Colussi et al. (2018b). The size of these embedded dimers sets a new length scale given by the effective scattering length, and we draw the analogy with Cooper pairing in BCS theory Strinati et al. (2018). We find that this length scale and the development of the bound pairs at unitarity provide an intuitive way to frame both the theoretical results of our model and the experimental results of Ref. Eigen et al. (2017) for the number of unbound atoms remaining after a fast-sweep projection. In order to make a quantitative comparison with the experiment, we include three-body losses phenomenologically in our many-body model by assuming an effective universal three-body loss-rate coefficient and by refitting the experimental estimate of this parameter.

We find good quantitative agreement with experimental data from Ref. Eigen et al. (2017) for the fastest ramp considered in that work over the full range of times where our model remains valid. However, for slower ramps we begin to deviate quantitatively from the experimental findings at later times. We argue that this deviation may be due to the presence of Efimov trimers or from the equilibrating effect of collisions both of which are not described in our model. This motivates further development of our theoretical model to include higher-order correlations, which remains a subject of ongoing study.

VI Acknowledgments

The authors thank Christoph Eigen, Zoran Hadzibabic, and Robert P. Smith for inspiring discussion and for providing experimental data. We also acknowledge Thomas Secker, Paul Mestrom, and Denise Braun for useful discussion. This work is supported by Netherlands Organization for Scientific Research (NWO) under Grant No. 680-47-623.