Two-Step Production of Resonant Bose-Einstein Condensates

TL;DR

This paper proposes a two-step method for producing resonant Bose-Einstein condensates that maximizes transfer probability by using an intermediate scattering length, addressing challenges like heating and losses.

Contribution

It introduces a novel two-step scheme with an optimal intermediate scattering length for efficient resonant BEC production.

Findings

Optimal intermediate scattering length is approximately 3.16N^{-2/3} times the harmonic oscillator length.

Maximum transfer probability is about 1.03N^{-1/6}.

Two-step process improves stability and efficiency over sudden quench methods.

Abstract

Producing a substantial and stable resonant Bose-Einstein condensate (BEC) has proven to be a challenging experimental task due to heating and three-body losses that may occur even before the gas comes to thermal equilibrium. In this paper, by considering only two-body correlations, we note that a sudden quench from small to large scattering lengths may not be the best way to prepare a resonant BEC. As an alternative, we propose a two-step scheme that involves an intermediate scattering length, between $0$ and $\infty$, which serves to maximize the transfer probability of $N$ bosons of mass $m$ in a harmonic trap with frequency $\omega$. We find that the intermediate scattering length should be $a\approx3.16N^{-2/3}\sqrt{\hbar/(m\omega)}$ to produce an optimum transition probability of $1.03N^{-1/6}$.

Equations

$$\Psi ({\mathbf{r}}_1,{\mathbf{r}}_2,\dots {\mathbf{r}}_{\mathbf{N}})=\prod _{\mathbf{i}=1}^{\mathbf{N}}{\varphi }_{\mathbf{a}}({\mathbf{r}}_{\mathbf{i}}).$$

$${\Psi }_{\mathrm{res}}({\mathbf{r}}_1,{\mathbf{r}}_2,\dots {\mathbf{r}}_{\mathbf{N}})=\prod _{\mathbf{i}=1}^{\mathbf{N}}{\psi }_{\infty }({\mathbf{r}}_{\mathbf{i}}).$$

$$P=|⟨\Psi |{\Psi }_{\mathrm{res}}⟩{|}^2={\left({\left|\int d^{3}r{\varphi }_a^{\ast }(\mathbf{r}){\psi }_{\infty }(\mathbf{r})\right|}^2\right)}^N$$

$${\rho }^2=\frac{1}{N}\sum _{i<j}^N{r}_{ij}^2.$$

$$H_{\mathrm{rel}}=-\frac{{\hslash }^2}{2m}\left[\frac{1}{{\rho }^{3N-4}}\frac{\partial }{\partial \rho }{\rho }^{3N-4}\frac{\partial }{\partial \rho }-\frac{{\Lambda }_{N-1}^2}{{\rho }^2}\right]+\frac{1}{2}m{\omega }^2{\rho }^2,$$

$$\Psi ={\rho }^{-(3N-4)/2}\sum _{\{\lambda \}}{F}_{\{\lambda \}}(\rho )Y_{\{\lambda \}}(\rho ;\Omega ),$$

$$\Psi ={\rho }^{-(3N-4)/2}{F}_{\nu }(\rho )Y_{\nu }(\rho ;\alpha ).$$

$$\left[-\frac{{\hslash }^2}{2m}\frac{d^2}{d{\rho }^2}+V^{\mathrm{diag}}(\rho )+V_{\nu }(\rho )\right]F_{\nu }(\rho )=E_{\mathrm{rel}}{F}_{\nu }(\rho ),$$

$$V^{\mathrm{diag}}(\rho )$$

$$V_{\nu }(\rho )$$

$$|a,n⟩={\rho }^{-(3N-4)/2}{F}_n^a(\rho ){\Phi }^a(\rho ;\Omega ).$$

$$V^{\infty }(\rho )=\frac{{\hslash }^2}{2m{\rho }^2}\left(\frac{9N^2}{4}+3c_0{N}^{8/3}\right)+\frac{1}{2}m{\omega }^2{\rho }^2,$$

$${\rho }^0$$

$${\rho }^{\infty }$$

$$V^0(\rho )$$

$$V^{\infty }(\rho )$$

$$E_{nK}=\hslash \omega \left(2n+K+\frac{3N-3}{2}\right),n=0,1,2,...$$

$$F^0(\rho )$$

$$F^{\infty }(\rho )$$

$${\left|⟨0,0,|\infty ,0⟩\right|}^2$$

$$\le {\left|\int d\rho F^0(\rho )F^{\infty }(\rho )\right|}^2\approx \exp \left(-1.3N^{4/3}\right),$$

$$V^a(\rho )\underset{N\gg 3}{\approx }{V}^0(\rho )+\frac{{\hslash }^2}{m}{d}_0{N}^{7/2}\frac{a}{{\rho }^3},$$

$$V^a({\rho }_{1n})$$

$$V^a({\rho }_{2n})$$

$$a^{\ast }$$

$$E^{\ast }$$

$$|{\Psi }^a(t)⟩$$

$$=\sum _{n=0}^{\infty }|a,n⟩⟨a,n|0,0⟩e^{-i{E}_{n}t/\hslash },$$

$$P(t)=|⟨\infty ,0|{\Psi }^a(t)⟩{|}^2,$$

$$⟨\infty ,0|{\Psi }^a(t)⟩=\sum _{n=0}^{\infty }⟨\infty ,0|a,n⟩⟨a,n|0,0⟩e^{-i{E}_{n}t/\hslash }.$$

$$⟨a,n|0,0⟩$$

$$\approx F^0({\rho }_{1n})\sqrt{\frac{d{E}_n}{dn}}\sqrt{\frac{1}{|\partial V^a/\partial \rho {|}_{{\rho }_{1n}}}},$$

$$⟨\infty ,0|a,n⟩$$

$$\approx (-1{)}^n{F}^{\infty }({\rho }_{2n})\sqrt{\frac{d{E}_n}{dn}}\sqrt{\frac{1}{|\partial V^a/\partial \rho {|}_{{\rho }_{2n}}}},$$

$$⟨{\Phi }^a|{\Phi }^0{⟩}_{\Omega }$$

$$⟨{\Phi }^{\infty }|{\Phi }^a{⟩}_{\Omega }$$

Full text

I Introduction

Recent experimental efforts have sought to prepare a Bose-Einstein condensate (BEC) of ultracold atoms in a regime where the two-body scattering length $a$ is infinite Makotyn14_NP ; Klauss17_PRL ; Eigen17_PRL ; Fletcher13_PRL ; Fletcher17_Science . Such a situation is termed a “resonant” (or sometimes “unitary”) BEC. It represents an unusual situation, inasmuch as the perturbative parameter $na^{3}$ – where $n$ is the number density – is no longer small and the usual field-theoretic ideas struggle to be useful. This circumstance has hatched a variety of alternative theoretical descriptions, which are in general agreement about the nature of the gas, yet differ in details Song09_PRL ; Cowell02_PRL ; Lee10_PRA ; Borzov ; Diederix11_PRA ; FZhou ; Yin13_PRA ; Ding17_PRA ; vanHeugten ; Rossi_PRA89 ; Sykes14_PRA ; Smith14_PRL ; Sze18_PRA .

On the experimental side, producing the resonant BEC is problematic, since the rate of three-body recombination grows rapidly with scattering length. In the resonant limit, this rate ultimately saturates, but at a large value that ensures the heating and ultimate destruction of the gas within milliseconds. Under these circumstances, a semblance of the approach to equilibrium can be teased out Eigen17_PRL ; Fletcher17_Science ; Eigen18_Nature , while the loss can be understood as a few-body process incorporating local physics of the gas Sykes14_PRA ; DIncao18_PRL ; Bedaque ; Braaten06_PhysRep .

In order to perform an experiment of this kind at all, the resonant BEC must therefore be produced quickly. A typical experimental protocol starts with the gas at a small value of scattering length, then rapidly ramps the value of a magnetic field near a Fano-Feshbach resonance so that $a\rightarrow\infty$ within microseconds. This represents the essentially instantaneous projection of the many-body state at small $a$ onto a collection of many-body states at $a=\infty$.

While this rapid ramp is essential for defining time zero of the nonequilibrium dynamics, it it not necessarily the best protocol for generating a true resonant BEC. To see this, at least qualitatively, it is useful to regard the gas within a mean-field-like description. Consider a gas of $N$ identical bosons, each initially in some single-particle orbital $\phi_{a}({\bf r})$, corresponding to the small initial scattering length $a$ (the function $\phi_{a}$ could be the ground state solution to the Hartree-Fock equations for the Bose system, for example). The many-body wave function is then, to a good approximation,

[TABLE]

Similarly, on resonance each atom can be regarded as belonging to some different orbital wave function $\psi_{\infty}({\bf r})$. This could be obtained approximately, for example, by performing a Hartree-Fock calculation using a renormalized scattering length $a\propto n^{-1/3}$ Ding17_PRA ; Sze18_PRA ; vonStecher07_PRA ; FZhou ; Lee10_PRA ; Song09_PRL . Thus, at least up to a certain approximation, the desired resonant BEC is described by

[TABLE]

Then the probability that the initial state $\Psi$ produces the resonant BEC states $\Psi_{\rm res}$, assuming that the fast ramp results in a projection, is given by the square of their overlap

[TABLE]

Unless each of these overlap integrals is very close to one, the product of $N$ of them will be vanishingly small for typical experimental circumstances with $N>10^{3}$. For this reason, it appears that, while the sudden ramp to $a=\infty$ produces an interesting, nonequilibrium gas of strongly-interacting bosons, it is unlikely to generate the desired resonant BEC.

In this paper we present an alternative scheme for preparing a resonant BEC, which proceeds in two steps. In a first step, the scattering length is jumped quickly from a low initial value $a_{1}\approx 0$ to a modest intermediate value $a_{2}$. The sudden increase in scattering length causes the BEC to expand; when it reaches the size of the resonant BEC, the scattering length is suddenly jumped from $a_{2}$ to $a=\infty$. For a properly-chosen value of the intermediate scattering length $a_{2}$, we show that the fraction of atoms converted into a resonant BEC can be non-negligible Klauss17_Thesis .

To describe and carry out calculations of this scheme, we focus on an isotropic, harmonically trapped BEC and employ a coordinate-based representation of the BEC wave function. This representation presents the BEC as a wave packet subject to an effective potential energy surface (PES) Bohn98_PRA , and it has recently been shown to make a reasonable description of the BEC on resonance Ding17_PRA ; Sze18_PRA . It presents the dynamics as the time evolution of a wave packet obeying a linear Schrödinger equation. In these terms the two-step process is reminiscent of vibrational wave packet dynamics in molecular physics Garraway02 . It is also amenable to analytic approximations, which will yield simple estimates for the optimum value of the intermediate scattering length $a_{2}$, as well as the approximate yield of atoms in the resonant BEC at the end of the two steps.

II Potential Energy Surfaces, Quench from Non-Interacting to Resonant BEC’s

Here, we summarize the theory, which is detailed in Ref. Sze18_PRA .

We exploit a coordinate representation of BEC that is expressed in terms of potential energy surfaces (PES’s) analogous to Born-Oppenheimer (B.-O.) curves in molecular physics. We define a single, collective coordinate, the hyperradius $\rho$, which represents the size of the condensate. This hyperradius can be expressed as the root-mean-squared interparticle spacing for any configuration of atoms Sorensen02_PRA :

[TABLE]

All remaining coordinates, collectively denoted by $\Omega$, span a hypersphere of radius $\rho$ in a $(3N-4)$-dimensional configuration space. If the center of mass coordinates are separated, then the Hamiltonian describing the relative motion is given by Sorensen02_PRA

[TABLE]

where $m$ is the atomic mass and $\omega$ is the angular frequency of the isotropic harmonic trap. Thus the kinetic energy has a radial part and an angular part, the latter given in general by the grand angular momentum $\Lambda_{N-1}^{2}$. Two-body interparticle coordinates are encoded in the angular component. Realistic two-body potentials between atoms Das04_PRA ; Das07_PRA ; Chakrabarti08_PRA ; Lekala14_PRA , or boundary conditions with realistic scattering lengths Sorensen02_PRA ; Sorensen03_PRA ; Sorensen04_JPB ; Sogo05_JPB may be applied to solve the hyperangular component of $H_{\rm{rel}}$.

Under the B.-O. approximation, the hyperradius $\rho$ is treated as the slow coordinate. That is, at each value of $\rho$, the Schrödinger equation, $H_{\rm{rel}}\Psi(\rho,\Omega)=E_{\rm{rel}}\Psi(\rho,\Omega)$ with $E_{\rm{rel}}$ the energy of the relative motion, is solved in the coordinates $\Omega$ to yield a set of eigenenergies $V_{\nu}(\rho)$. A coupled set of differential equations is obtained if we expand the wave function $\Psi$ in adiabatic hyperangular basis,

[TABLE]

for some set of radial expansion functions $F_{\{\lambda\}}$, and $Y_{\{\lambda\}}$ are the eigenstates of $\Lambda_{N-1}^{2}$. However, by applying the B.-O. approximation we assume that the hyperradial kinetic energy operator does not affect $Y_{\{\lambda\}}$. Hence the differential couplings can be neglected. For simplicity, we also reduce the collective set of quantum numbers $\{\lambda\}$ to a single quantum number $\nu$, describing excitations in a single hyperangle $\alpha$, where $\sin\alpha=r_{12}/\left(\sqrt{2}\rho\right)$, which incorporates two-body correlations. Thus, we write the wave function as

[TABLE]

Using a single adiabatic function, the Schrödinger equation becomes a single ordinary differential equation in $\rho$:

[TABLE]

where

[TABLE]

$V^{\rm diag}$ is the diagonal potential whose ground state supports the non-interacting condensate wave function. $V_{\nu}$ represents the interaction potential since the state $|\nu\rangle$ is defined by the boundary conditions based on the scattering length which describes the two-body interaction. Exact calculation of the matrix element $\langle\nu|\Lambda_{N-1}^{2}|\nu\rangle$, which involves integration over the entire hypersphere, is not trivial. In this paper, we use the results of Ref. Sze18_PRA , where some convenient approximations have been made to obtain a meaningful outcome even for $a=\infty$. Also, we apply only this single matrix element that is representative of the interaction of the lowest hyperangular state of the condensate with any given $a$. Thus, for any scattering length $a$, we find a B.-O. potential $V^{a}(\rho)=V^{\rm diag}(\rho)+V_{\nu}(\rho)$ with associated hyperangular wave function $\Phi^{a}(\rho;\Omega)$. Vibrational states in the PES $V^{a}(\rho)$ constitute the radial wave functions $F_{n}^{a}(\rho)$, each vibration $n$ describing a breathing mode excited above the ground state condensate with $n=0$. The states relevant to our model are, therefore, defined by the scattering length $a$ and the number of breathing quanta,

[TABLE]

Figure  1 shows the B.-O. PES’s for the non-interacting ($a=0$) and resonant ($a=\infty$) cases for a gas of $N=100$ atoms. With $a=0$, $V^{0}(\rho)=V^{\rm diag}(\rho)$. This PES, the lowest curve on the left, is exact. The topmost curve on the right is an approximate surface for the resonant limit. Considering only two-body correlations explicitly, this surface is constructed based on Ref. Sze18_PRA . For the large $N$ limit, this is given by

[TABLE]

where $c_{0}\approx 2.122$ is a constant determined by the root of some transcendental equation. For realistic values of $N>10^{2}$, the centrifugal term with $9N^{2}/4$ can be safely neglected. The ground states of these PES’s represent the non-interacting and resonant BEC’s. Near their minima Sze18_PRA ,

[TABLE]

where $a_{ho}=\sqrt{\hbar/(m\omega)}$, we approximate these potentials as harmonic oscillators:

[TABLE]

In both cases, the excitation frequency of the radial breathing modes considered is exactly twice the trap frequency, $\omega_{b}=2\omega$. For non-interacting bosons, the energies are well-known and are given by Smirnov77

[TABLE]

and $K=0,1,2,...$ is the quantum number associated with the hyperangular component. For the resonant gas, the $2\omega$ frequency was anticipated by symmetry considerations in Refs. Castin_Physique5 ; Werner_PRA74 . Without considering three-body or higher order correlations, these references also emphasize that the B.-O. approximation is exact in the $a=\infty$ limit. Corrections beyond the B.-O. approximation arise because the adiabatic wave functions $\Phi$ change from one value of $\rho$ to the next. But this change is only effective if $\rho$ changes significantly on the scale of $a$, i.e., the corrections are of order $\rho/a$ and vanish in the infinite scattering length limit. Therefore, if the atoms could be prepared in the state $F^{\infty}\Phi^{\infty}$ that we describe, this state would be stable against non-adiabatic transitions to whatever other states there are that could lead to heating, loss, etc. This stability is likely reduced if we were to include explicit three-body correlations in the wave function.

From the harmonic oscillator nature of the potential curves in Eqs. (15) and  (16), the expected ground state hyperradial wave functions should be Gaussians centered at the minima and with root-mean-squared width of $a_{ho}/\sqrt{2}$:

[TABLE]

The unnormalized Gaussian functions $F^{0}$ and $F^{\infty}$ for $N=100$ are illustrated as Gaussian-shaped humps at the bottom of the $a=0$ and $a=\infty$ PES’s, respectively, in Fig. 1. From this picture, we see that the centers are far away from each other such that quenching the gas suddenly from $a=0$ to $a=\infty$ will yield a low transfer probability. That is, the probability of the atoms landing in the resonant BEC state $F^{\infty}$, upon a direct quench, is

[TABLE]

which is negligible for large $N$.

III The Two-Step Scheme

III.1 Franck-Condon Factors

The tiny overlap between $F^{0}$ and $F^{\infty}$ suggests that direct projection from $a=0$ to $a=\infty$ will not yield a good amount of resonant BEC. We then seek an intermediate state with finite, nonzero value of $a$. Such a PES, $V^{a}(\rho)$, is shown as the intermediate curve in Fig. 2. A good candidate for $V^{a}(\rho)$ is one that supports a set of vibrational excitations $n$ so that $|a,n\rangle$ has a good overlap with both the $|0,0\rangle$ and $|\infty,0\rangle$ states as shown in Fig. 2. Real BEC experiments have $N>10^{4}$ atoms. In Fig. 2, we use $N=100$ as an illustrative example. For larger $N$, $\rho^{0}$ and $\rho^{\infty}$ grow farther apart. One then needs to use higher vibrational states (with larger number of nodes $n$) to optimize the overlaps $\langle 0,0|a,n\rangle$ and $\langle a,n|\infty,0\rangle$. These squared overlaps $|\langle 0,0|a,n\rangle|^{2}$ and $|\langle a,n|\infty,0\rangle|^{2}$ are called Franck-Condon (FC) factors. Numerical calculations of these FC factors $|\langle 0,0|a,n\rangle|^{2}$ and $|\langle a,n|\infty,0\rangle|^{2}$, and $|\langle 0,0|a,n\rangle\langle a,n|\infty,0\rangle|^{2}$ are reflected as color-map plots in Fig. 3 for $N=100$; the x-axis is the scattering length, y-axis the vibrational state $n$, and the color indicates the transition probability. In general, for the first step from the non-interacting to the intermediate, the optimum transition occurs when $a$ is small and for low $n$ states, decreasing quickly with increasing $a$ and $n$ as shown in Fig. 3(a). For the second step from intermediate to final, the transition is optimum when $a$ and $n$ are larger, and diminishes slowly with decreasing $a$ and increasing $n$ as in Fig. 3(b). These two steps cannot be individually at their maxima under the same conditions. However, the best overall yield occurs when $a$ is still small relative to the oscillator length and for higher vibrational states. This is true for any large values of $N$. Further, the two-step transition probabilities seem to decrease as a function of $N$. See the transition probability for $N=1000$ in Fig. 4.

III.2 The Optimum Intermediate State

Since the intermediate state will have a small value of $a$, we can use a perturbative approximate expression for $V^{a}$. In the limits of perturbative $a\ll a_{ho}$ and large $N$, this is given by

[TABLE]

where $d_{0}=(3/4)\sqrt{3/\pi}\approx 0.733$. This potential can be well utilized by considering the classical inner and outer turning points, $\rho_{1n}$ and $\rho_{2n}$, of $V^{a}$ at some particular energy $E_{n}$ of $F_{n}$. At high vibrational states $n$, $\rho_{1n}$ and $\rho_{2n}$ can be approximated through

[TABLE]

where $V^{a}$ is dominated by the interaction term at small $\rho$, and by the trapping potential at large $\rho$. An effective two-step scheme is illustrated in Fig. 2. It is achieved when the inner turning point $\rho_{1n}$ of $V^{a}$ is near $\rho^{0}$ and the outer turning point $\rho_{2n}$ is near $\rho^{\infty}$. Thus, with $\rho_{1n}\approx\rho^{0}$, $\rho_{2n}\approx\rho^{\infty}$ and Eqs. (23) and  (22), the state which would give the maximum Franck-Condon overlap is one whose scattering length and energy are

[TABLE]

Using these approximations for $N=100$, the results are $a^{*}=0.145a_{ho}$ and $E^{*}=571.2\hbar\omega$ which are close to the exact calculations of $a^{*}=0.0859a_{ho}$ and $E^{*}=598.9\hbar\omega$, the latter set of values can be visually estimated through Figs. 3 and 2. Expressions  (24) and (25) become better estimates for larger $N$. For $N=1000$, the predicted results are $a^{*}=0.0316a_{ho}$ and $E^{*}=1.26(10^{4})\hbar\omega$, and the numerical computations give $a^{*}=0.0332a_{ho}$ and $E^{*}=1.28(10^{4})\hbar\omega$.

While this static picture provides overall orientation, it does not describe the dynamics involved. Roughly, upon the initial projection from $a=0$ to the intermediate value $a^{*}$, a wave packet is formed at $\rho_{1n}$. In approximately one half of the trap period, this wave packet propagates to $\rho_{2n}$, giving the condensate its maximum radial extent and preparing it for projection onto the resonant BEC state.

III.3 Wave Packet Dynamics

To describe the time dynamics, we express the initial state after the first step as a wave packet expanded in the basis of the vibrational states of the intermediate potential

[TABLE]

where at time $t=0$, $\Psi^{a}$ is at the ground state of the non-interacting potential with total energy $E\approx 3N\hbar\omega/2$. The probability the projection of the wave packet onto the desired resonant BEC ground state is given by

[TABLE]

where

[TABLE]

After extracting the most appropriate choice for the intermediate $a$, we compute this transition probability at different times with the unitary BEC model found in Ref.Sze18_PRA for $N=100$ and $N=1000$. Figures 5(a) and (b) show that the first maximum transition occurring at around $t_{m}\approx\pi/(2\omega)$, with $48\%$ transfer probability for $N=100$ and $36\%$ for $N=1000$. It takes about half a period, $T/2$, for the BEC to expand to resonance starting from the left side of the $V^{a}$; the breathing mode frequency is close to $2\omega$, thus the dwell time is $t_{m}\approx T/2=\pi/\omega_{b}=\pi/(2\omega)$.

Figure  6 shows how the size of the BEC with $N=100$ atoms, expressed in terms of the mean hyperradius $\langle\rho\rangle$, is changing over time. It starts with $\rho=\rho^{0}$, the size of the non-interacting gas, and reaches $\rho=\rho^{\infty}$, the size of the resonant BEC, at $t\approx t_{m}$. The peaks of $P(t)$ and $\langle\rho\rangle$ decrease slowly over time as the wave packet gradually dephases. It is, therefore, worthwhile to instigate the second projection, to resonance, at time $t=T/2$.

IV Large $N$ Limit

In calculating the $P(t)$ numerically, we notice that $P(t_{m})$ decreases with $N$. Determining how $P(t_{m})$ scales with $N$ is extremely useful. Here, we outline a method to get a good estimate for this scaling. The details are found in the Appendix, and the final result turns out to be simple.

Using the results from Appendix A and B, the overlap integrals in Eq. (28) are approximated to be

[TABLE]

where $dn/dE_{n}$ is the density of vibrational states in the intermediate potential. For the hyperangular parts of the wave function we approximate

[TABLE]

since $N$ is large and $a/\rho_{1n}$ is small.

Next, we convert the discrete sum in Eq. (28) into a continuum integral over the energy and evaluate it at $t=t_{m}$ around which the maximum transfer occurs. See Appendix C for details. The resulting transition amplitude is

[TABLE]

where $c_{0}$ and $d_{0}$ are defined in Eqs. (12) and  (21). Plots of $P(t_{m})$, calculated in this way, for different $N$ are shown in Fig. 7 as a function of the intermediate scattering length $a$. We see that the estimated maximum transfer for $N=10^{3}$ is $\sim 33\%$, which is close to what the exact calculation gives. The inset in Fig. 7 shows the sensitivity of the transition probability to the intermediate $a$ for $N=10^{5}$. The intermediate $a$ should at least be within $0.4\%$ from the optimum to get at least half of the maximum transfer. By maximizing Eq. (33) with respect to $a$, or by using Eqs. (24), (13), and (14), the optimum scattering length is found to be

[TABLE]

And the maximum transfer is

[TABLE]

To put this into context, for 85Rb in a trap with frequency $\omega=2\pi\times 10\,\text{Hz}$, the oscillator length is $a_{ho}=6.51\times 10^{4}\,a_{0}$. Starting with $N=10^{5}$ non-interacting atoms in the trap, the two-step process would be optimized for a scattering length of $a^{*}\approx 95.4\,a_{0}$.

V Conclusions and Prospects

We have presented a protocol designed to implant a nontrivial fraction of the trapped atoms into a resonant BEC. It remains to be understood what the consequences of this preparation step will be. It is not clear, for example, what further reorganization of the atoms might be necessary for the gas to resemble an equilibrium resonant BEC. It is equally unclear at present how three-body losses would differ in the resonant BEC than in a gas of equivalent density. A useful initial experiment might be to prepare the resonant BEC as proposed here, and compare its dynamics to that of a gas of equal initial density as the resonant BEC, but jumped suddenly to resonance.

This experiment would unfortunately be clouded by another issue. Consider, for example, that starting from a non-interacting BEC of $N=10^{4}$ atoms, our protocol is expected to transfer only one fifth of them to the resonant BEC. What becomes of the rest? They are presumably projected onto other quantum mechanical states of the system, each of which has its own dynamics and three-body loss rates. To address this, it is necessary to formulate a reliable theory of excited states, in our case in the hyperangular degrees of freedom. This pursuit is currently underway.

Acknowledgements.

M. W. C. S. and J. L. B. were supported by the JILA NSF Physics Frontier Center, grant number PHY-1734006.

Appendix A Franck-Condon Factors Using the Reflection Formula

Here, we evaluate overlap integrals

[TABLE]

Leading contribution to the Franck-Condon factors comes from the overlap of wave functions at the classical turning points, where the wave functions $F_{n}^{a}$ are sharply peaked. In between the turning points, the wave functions are highly-oscillating. Yet, we can consider that the projections of $F_{n}^{a}$ to $F^{0}$ and $F^{a}$ are still localized to the turning points since the latter wave functions are also localized (or close to zero where $F_{n}^{a}$ is wildly oscillating). The idea that the Franck-Condon factors can be estimated from properties of the potential near the turning points goes back to the early days of quantum mechanics Condon_PR28 ; Winans_28 . It is widely used in theories of optical and Raman transitions in molecules, and recently to photoassociation of cold atoms as well Suominen96 ; Weiner99 ; Bohn99 ; Julienne96 ; Boisseau00 . Out of these types of molecular spectroscopy studies, the reflection formula was developed Julienne96 ; Jablonski45 , which we will adapt.

We first express $F_{n}^{a}$ in terms of the energy-normalized wavefunction $F_{E}$ through

[TABLE]

which leads to $F_{n}^{a}=\sqrt{dE_{n}/dn}F_{E}$. Casting $F_{E}$ into phase-amplitude form, after Milne Milne30 ,

[TABLE]

where the amplitude $\zeta$ and phase $\beta$ satisfy

[TABLE]

with the wave vector

[TABLE]

The rapid oscillations of $F_{E}$ in  (39) will have negligible effect on the integrals in Eqs. (36) and (37), where $F_{n}^{a}$ is expressed in terms of $F_{E}$, except when $\rho$ is near a turning point which is also a point of stationary phase. Away from a turning point, it is sufficient to use the WKB approximations for the amplitude and phase:

[TABLE]

Near a turning point $\rho_{t}$, we expand the Milne phase to second order

[TABLE]

Now, with

[TABLE]

the integrand $F_{n}^{a}F^{0}$ is sharply localized around $\rho_{1n}$, the classical inner turning point. Thus,

[TABLE]

To evaluate the last integral, we use the formula

[TABLE]

Finally, we arrive at

[TABLE]

The other overlap factor (37) can be approximated in a similar fashion; it is given by

[TABLE]

where the $(-1)^{n}$ accounts for the sign of the rightmost amplitude around the outer turning point $\rho_{2n}$ of the vibrational state if we set the leftmost amplitude around $\rho_{1n}$ always positive as expressed in Eq. (53).

Appendix B Overlap between LOCV Hyperangular Wave Functions

To give a complete picture of the overlap between wave functions, the angular overlaps $\langle\Phi^{a}|\Phi^{0}\rangle_{\Omega}$ and $\langle\Phi^{\infty}|\Phi^{a}\rangle_{\Omega}$ should also be considered. Real calculation involves $3N-4$ dimensional integrals since this is the size of the hyperangular space. However, here, we only consider the one hyperangle, $\alpha$, that describes the two-body interactions, and the large $N$ case.

We start with a symmetrized Jastrow-type basis,

[TABLE]

where $\alpha_{ij}$ is parametrically related to the coordinate distance between to particles, $r_{ij}$ through $r_{ij}=\sqrt{2}\rho\sin\alpha_{ij}$; the function $\phi_{\nu}$ satisfies the Bethe-Peierls boundary condition which describes what happens when two particles are close to each other. The other boundary condition is set by treating $|\phi_{\nu}|^{2}$ as a pair correlation function such that if two atoms are more than distance $r_{d}=\sqrt{2}\rho\sin\alpha_{d}$ apart, then they become uncorrelated or $|\phi(\alpha_{ij}\geq\alpha_{d})|^{2}=1$. Therefore, within a region bounded by $\alpha_{d}$, there is on the average only one other atom (out of $N-1$) which can be seen by a fixed atom, or

[TABLE]

where $d\Omega=d\Omega_{N-1}=4\pi d\Omega_{\alpha}d\Omega_{N-2}$, and $d\Omega_{\alpha}=\sin^{2}\alpha\cos^{3N-7}\alpha\,d\alpha$. If $\alpha_{d}=\pi/2$, then the right side of  (56) should be one. The full form of the pair correlation function $g_{2}$ can be written as

[TABLE]

which is hard to evaluate. To lowest order, however, it is approximated to be $g_{2}(\alpha)=|\phi_{\nu}(\alpha)|^{2}$. This whole procedure outlined above describes a lowest order constraint variational (LOCV) method in hyperspherical coordinates; details can be found in Ref.Sze18_PRA . Given $\rho$ and the scattering length $a$, one can then find $\alpha_{d}$ and $\phi_{v}$. The angle $\alpha_{d}$ becomes extremely small as $N$ increases. Hence $\phi(\rho;\alpha_{ij})$ is one in large region of $\alpha_{ij}$ - this is an approximation that leads to $g_{2}(\alpha)=|\phi_{\nu}(\alpha)|^{2}$.

In the following derivations, we will also treat all the pair wave functions $\phi(\rho,\alpha_{i^{\prime}j^{\prime}})$ equivalent to unity, except one pair namely, $\phi(\rho,\alpha_{12})=\phi(\rho,\alpha)$. So,

[TABLE]

where the $\mathcal{N}$’s are some normalization constants so that $\langle\Phi^{0}|\Phi^{0}\rangle_{\Omega}=1$, $\langle\Phi^{a}|\Phi^{a}\rangle_{\Omega}=1$, and $\langle\Phi^{\infty}|\Phi^{\infty}\rangle_{\Omega}=1$, and Sze18_PRA

[TABLE]

The wave functions $\phi^{a}$ and $\phi^{\infty}$ identically approach unity for $\alpha>\alpha_{a}$ and $\alpha>\alpha_{\infty}$, which are given by

[TABLE]

Note that $\alpha_{a}$ and $\alpha_{\infty}$ are extremely small for large $N$ so that the integrals in Eqs. (58) and (59) are over large part of the $\alpha$-space where $\phi^{a}$ and $\phi^{\infty}$ are unity. The constants $A$ and $B$ are determined from the continuity boundary condition at $\alpha_{a}$ and $\alpha_{\infty}$:

[TABLE]

We then find

[TABLE]

Finally, after a series of algebraic steps and careful bookkeeping of $N$-scaling of the relevant parameters, we find

[TABLE]

which are our approximations for $\langle\Phi^{a}|\Phi^{0}\rangle_{\Omega}$ and $\langle\Phi^{\infty}|\Phi^{a}\rangle_{\Omega}$, respectively. For large $N$, these quantities are both essentially equal to one.

Appendix C The Transition Amplitude

We evaluate the transition amplitude at $t=t_{m}\approx\pi/(2\omega)$ at large $N$. In terms of the Franck-Condon factors derived in Appendix A, we write the transition amplitude as

[TABLE]

with $\omega_{n}\approx(2+\Delta_{n})n\omega$, where $\Delta_{n}<1$ ($\Delta_{n}\ll 1$ for small $a$). Thus,

[TABLE]

Also, using Eqs. (22) and  (23),

[TABLE]

Converting the discrete sum into an integral over energy, $\sum_{n}\rightarrow\int dE$, and using the form of $F^{0}$ and $F^{\infty}$ in Eqs. (18) and (19), and noting that the resulting integrand is strongly peaked at $E^{*}\approx\sqrt{3c_{0}}N^{4/3}\hbar\omega/2\approx 1.26N^{4/3}\hbar\omega$ (see Eqs. (25) and (14)), we get

[TABLE]

with $\rho_{1}\approx(\frac{\hbar^{2}}{m}d_{0}N^{7/2}a)^{1/3}E^{-1/3}$ and $\rho_{2}\approx\sqrt{2E/(m\omega^{2})}$ from Eqs. (22) and  (23). Now, $F^{\infty}(\rho_{2}(E))$ is a peaky function of $E$. We can then use the saddle point approximation to solve the integral in Eq. (80):

[TABLE]

Finally, expressing $\rho^{0}$ and $\rho^{\infty}$ in terms of $N$,

[TABLE]