Number of closed-channel molecules in the BEC-BCS crossover

# Number of closed-channel molecules in the BEC-BCS crossover

F. Werner Present address: Department of Physics, University of Massachusetts, Amherst, MA 01003, USALaboratoire Kastler Brossel, Ecole Normale Supérieure, UPMC and CNRS, 24 rue Lhomond, 75231 Paris Cedex 05, France    L. Tarruell Present address: Institute for Quantum Electronics, ETH Zürich, 8093 Zürich, Switzerland Laboratoire Kastler Brossel, Ecole Normale Supérieure, UPMC and CNRS, 24 rue Lhomond, 75231 Paris Cedex 05, France    Y. Castin Laboratoire Kastler Brossel, Ecole Normale Supérieure, UPMC and CNRS, 24 rue Lhomond, 75231 Paris Cedex 05, France
###### Abstract

Using a two-channel model, we show that the number of closed-channel molecules in a two-component Fermi gas close to a Feshbach resonance is directly related to the derivative of the energy of the gas with respect to the inverse scattering length. We extract this quantity from the fixed-node Monte Carlo equation of state and we compare to the number of closed-channel molecules measured in the Rice experiment with lithium [Partridge et al., Phys. Rev. Lett. 95, 020404 (2005)]. We also discuss the effect of a difference between the trapping potentials seen by a closed-channel molecule and by an open-channel pair of atoms in terms of an effective position-dependent scattering length.

###### pacs:
03.75.SsDegenerate Fermi gases - 67.90.+z Other topics in quantum fluids and solids

## 1 Introduction

It is now experimentally possible to prepare two-component atomic Fermi gases at low temperature with a fully controlable interaction strength, in a regime where the interaction range is negligible and the atomic interactions are characterized by a single parameter, the s-wave scattering length a between two opposite “spin” atoms. The value of a can be adjusted at will thanks to a magnetically induced Feshbach resonance. The weakly interacting limits k_{F}a\to 0^{+} and k_{F}a\to 0^{-}, where k_{F} is the Fermi wavevector of the gas, correspond respectively to the BEC limit (a Bose-Einstein condensate of dimers, observed in BEC_Inns (); BEC_JILA (); BEC_MIT (); SalomonCrossover (); Hulet ()) and the BCS limit (a “condensate” of Cooper pairs approximately described by the BCS theory). Furthermore, the experiments can access the so-called crossover regime between BEC and BCS, where the gas is strongly interacting (k_{F}|a|\gtrsim 1), which includes the celebrated unitary limit k_{F}|a|=\infty, where the gas acquires fully universal properties Thomas1 (); Jin (); Salomon1 (); GrimmCrossover (); Grimm_gap (); ThomasVirielExp (); GrimmModes (); Ketterle_vortex (); JinPotentialEnergy (); thomas_entropie (); HuletPolarized (); Ketterle_unba ().

In this context, early theoretical studies Timmermans_phys_rep (); Holland0 (); Holland () put forward a many-body Hamiltonian which accurately models the microscopic two-body physics of the Feshbach resonance. The interaction potential is described by two channels, an open channel and a closed channel. The atoms exist in the form of fermionic particles in the open channel and in the form of bosonic short range molecules in the closed channel. Two atoms may be converted into a short range closed-channel molecule and vice versa due to the interchannel coupling. The corresponding Hamiltonian is

 \displaystyle H_{2}=\int d^{3}r\left\{\sum_{\sigma=\uparrow,\downarrow}\psi_{% \sigma}^{\dagger}(\mathbf{r})\left[-\frac{\hbar^{2}}{2m}\Delta+U(\mathbf{r})% \right]\psi_{\sigma}(\mathbf{r})\right.\\ \displaystyle\left.+\psi_{b}^{\dagger}(\mathbf{r})\left[E_{b}(B)-\frac{\hbar^{% 2}}{4m}\Delta+U_{b}(\mathbf{r})\right]\psi_{b}(\mathbf{r})\right\}\\ \displaystyle+\Lambda\int d^{3}r_{1}d^{3}r_{2}\,\chi(\mathbf{r}_{1}-\mathbf{r}% _{2})\left\{\psi_{b}^{\dagger}[(\mathbf{r}_{1}+\mathbf{r}_{2})/2]\psi_{% \downarrow}(\mathbf{r}_{1})\psi_{\uparrow}(\mathbf{r}_{2})+{\rm h.c.}\right\}% \\ \displaystyle+g_{0}\int d^{3}r_{1}d^{3}r_{2}d^{3}r_{3}d^{3}r_{4}\,\chi(\mathbf% {r}_{1}-\mathbf{r}_{2})\chi(\mathbf{r}_{3}-\mathbf{r}_{4})\\ \displaystyle\delta\left(\frac{\mathbf{r}_{1}+\mathbf{r}_{2}}{2}-\frac{\mathbf% {r}_{3}+\mathbf{r}_{4}}{2}\right)\psi_{\uparrow}^{\dagger}(\mathbf{r}_{1})\psi% _{\downarrow}^{\dagger}(\mathbf{r}_{2})\psi_{\downarrow}(\mathbf{r}_{3})\psi_{% \uparrow}(\mathbf{r}_{4}) (1)

where the atomic fields \psi_{\sigma}(\mathbf{r}) obey the usual fermionic anticommutation relations and the field \psi_{b}(\mathbf{r}) describing the closed-channel molecules obeys the usual bosonic commutation relations. \Lambda gives the coupling between the closed and open channels, with a short range cut-off function \chi(\mathbf{r}), of range b, and such that \int d^{3}r\,\chi(\mathbf{r})=1. A typical value for b is in the nanometer range 111 b can be estimated by the van der Waals length, that is the length that one forms with \hbar, m and the C_{6} coefficient of the van der Waals interaction Koehler_review (). . The coupling constant g_{0} represents the background atomic interaction in the open channel, which is conveniently modeled by a separable potential with the same cut-off function \chi. The energy E_{b} is the energy of a closed-channel molecule, in the absence of the coupling \Lambda, counted with respect to the dissociation limit of the open channel. This energy is adjusted with a magnetic field B, using the fact that the different magnetic moments in the open and closed channels lead to a differential Zeeman shift. The resulting effective magnetic moment is

 \mu_{b}\equiv\frac{dE_{b}}{dB}. (2)

U(\mathbf{r}) and U_{b}(\mathbf{r}) are the trapping potentials experienced by the atoms and the closed-channel molecules respectively.

On the contrary, recent Monte Carlo simulations of the many-body problem in the BEC-BCS crossover simply use single channel model Hamiltonians DeanLee (); bulgacQMC (); zhenyaPRL (); zhenyaNJP (); Juillet (); BulgacCrossover (); zhenyas_crossover (). In most unbiased Quantum Monte Carlo simulations, a lattice model is used:

 \displaystyle H_{1}=\sum_{\mathbf{r}}b^{3}\left\{\sum_{\sigma=\uparrow,% \downarrow}\psi_{\sigma}^{\dagger}(\mathbf{r})\left[-\frac{\hbar^{2}}{2m}% \Delta+U(\mathbf{r})\right]\psi_{\sigma}(\mathbf{r})\right.\\ \displaystyle\left.+g_{0}\psi_{\uparrow}^{\dagger}(\mathbf{r})\psi_{\downarrow% }^{\dagger}(\mathbf{r})\psi_{\downarrow}(\mathbf{r})\psi_{\uparrow}(\mathbf{r}% )\right\}. (3)

Here b is the lattice spacing (which coincides with the “range” of the discrete delta interaction potential), \Delta is a discrete representation of the Laplacian, and g_{0} is a coupling constant adjusted to have the desired scattering length: As was shown e.g. in zhenyaPRL (); zhenyaNJP (); YvanHoucheslowD (); boite (),

 \frac{1}{g}-\frac{1}{g_{0}}=\int_{[-\pi/b,\pi/b]^{3}}\frac{d^{3}k}{(2\pi)^{3}}% \frac{1}{2\epsilon_{\mathbf{k}}} (4)

where g=4\pi\hbar^{2}a/m is the effective coupling constant and \epsilon_{\mathbf{k}} the energy of a single particle with wavevector \mathbf{k} on the lattice. In fixed-node Monte Carlo calculations, one rather uses an interaction potential in continuous real space, e.g. a square well potential Pandha (); Giorgini (); Reddy (); Lobo (); Lobo_unb (); Giorgini_unb ().

There is now a consensus about the fact that, in the zero-range limit, there is a convergence of predictions for the two models H_{1} and H_{2}, in particular for the equation of state of the gas. The zero-range limit means, for single channel models, that k_{F}b\ll 1, b\ll\lambda, and b\ll|a|, where k_{F} is the Fermi wavevector and \lambda is the thermal de Broglie wavelength.222 The assumption b\ll|a| is not necessary on the BCS side if one restricts to macroscopic observables such as the total energy, but it is necessary for microscopic observables such as the ones considered in section 3.1. For two-channel models, one has to impose the same conditions, not only for the interaction range b but also for the effective range r_{e} (one indeed has |r_{e}|\gg b on narrow Feshbach resonances PetrovBosons (); YvanVarenna ()); we shall see in appendix A that one also needs to impose a condition on the background scattering length, and that all conditions are satisfied in the experiment Hulet ().333 A simple formal definition of the zero-range limit for the homogeneous gas is to take the limit of vanishing density k_{F}\to 0 while tuning the magnetic field B in such a way that 1/(k_{F}a) remains constant, and keeping constant all parameters of the Hamiltonian other than B ; one then also has to take the limit T\to 0 so that T/T_{F} remains constant.

One can however consider observables which exist in reality and are included in the two-channel model, but which are absent in a single channel model. The most natural example is the number of closed-channel molecules, which was recently measured at Rice using laser molecular excitation techniques Hulet ().

In this paper, we show that the equilibrium number of closed-channel molecules N_{b} is directly related to the equation of state of the gas, more precisely to the derivative of the gas energy (or free energy at non-zero temperature) with respect to 1/a. Paradoxically, we can thus use the equation of state calculated in Giorgini () within a single channel model, in order to predict N_{b} in the crossover and compare to Rice measurements, see section 2. Since the derivative of the energy with respect to 1/a is related to other observables involving atomic properties only Tan_nkREVTEX (); Tan_dEREVTEX (), it is also possible in principle to access N_{b} by a pure atomic measurement, see section 3. We conclude in section 4.

## 2 Prediction for the number of closed-channel molecules

### 2.1 General result

We first suppose that the system is in an arbitrary eigenstate |\psi\rangle of the two-channel Hamiltonian H_{2}, with eigenergy E. From Hellmann-Feynman theorem we have

 \frac{dE}{dB}=\langle\psi|\frac{dH_{2}}{dB}|\psi\rangle. (5)

The only magnetic field dependent part of the Hamiltonian H_{2} is the bare closed-channel molecule energy E_{b}(B) so that

 \frac{dE}{dB}=\mu_{b}N_{b} (6)

where

 N_{b}=\int d^{3}r\langle\psi_{b}^{\dagger}(\mathbf{r})\psi_{b}(\mathbf{r})\rangle (7)

is the mean number of closed-channel molecules.

We then eliminate the magnetic field by using as a parameter the scattering length a experienced by the trapped atoms rather than B. In the general case where the trapping potential U_{b}(\mathbf{r}) for the closed-channel molecule differs from twice the trapping potential U(\mathbf{r}) experienced by an open-channel atom, the scattering length of two atoms around point \mathbf{r} depends on \mathbf{r}, see appendix B. To lift the ambiguity we thus take for a the scattering length in the trap center taken as the origin of coordinates, \mathbf{r}=\mathbf{0}. We then rewrite (6) as

 N_{b}=\frac{dE}{d(1/a)}\,\frac{d(1/a)}{dB}\,\frac{1}{\mu_{b}}. (8)

Assuming that \mu_{b} is magnetic field independent, a as a function of B can be calculated explicitly for H_{2} by solving the zero-energy two-body scattering problem in free space (see appendix A), and by including the energy shifts due to the trapping potentials U_{b}(\mathbf{0}) and U(\mathbf{0}), as explained in the appendix B. The function takes the form

 a=a_{\rm fs}(B-\delta B_{0}) (9)

where we have introduced the function giving the free space scattering length as a function of B,

 a_{\rm fs}(B)=a_{\rm bg}\left(1-\frac{\Delta B}{B-B_{0}}\right). (10)

Here a_{\rm bg} is the background scattering length, B_{0} is the location of the Feshbach resonance in the absence of trapping potential, \delta B_{0}=[2U(\mathbf{0})-U_{b}(\mathbf{0})]/\mu_{b} the shift of the resonance location due to the trapping potential, and \Delta B the resonance width.

Equation (8) is directly applicable at zero temperature. At non-zero temperature, one has to take a thermal average of (8), keeping in mind that the mean value of a derivative is not necessarily equal to the derivative of the mean value. In the canonical ensemble one can check the exact relation Bogoliubov_livre ()

 \left\langle\frac{dE}{d(1/a)}\right\rangle=\left(\frac{dF}{d(1/a)}\right)_{T}=% \left(\frac{d\langle E\rangle}{d(1/a)}\right)_{S} (11)

where \langle\ldots\rangle is the thermal average, the derivative of the free energy F is taken for a fixed temperature T, and the derivative of \langle E\rangle is taken for a fixed entropy S. We thus have in the canonical ensemble

 \langle N_{b}\rangle=\left(\frac{d\langle E\rangle}{d(1/a)}\right)_{S}\,\frac{% d(1/a)}{dB}\,\frac{1}{\mu_{b}}. (12)

### 2.2 Analytical results in a trap in limiting cases

We now restrict to a spin balanced gas where the number of particles is equal to N/2 in each spin component. We also assume that the gas is at zero temperature (or at temperatures much smaller than the Fermi temperature T_{F}), in a harmonic trap

 U(\mathbf{r})=U({\bf 0})+\frac{1}{2}m\sum_{\alpha=x,y,z}\omega_{\alpha}^{2}r_{% \alpha}^{2}, (13)

in the macroscopic regime where k_{B}T_{F} is much larger than the oscillation quanta \hbar\omega_{\alpha}. The usual definition of the Fermi temperature T_{F} and of the Fermi wavevector k_{F} is then

 k_{B}T_{F}=\frac{\hbar^{2}k_{F}^{2}}{2m}=(3N)^{1/3}\hbar\bar{\omega} (14)

where \bar{\omega} is the geometric mean of the three oscillation frequencies \omega_{\alpha}. In this macroscopic regime, we rely on the local density approximation to calculate dE/d(1/a) for the trapped gas: This is expected to be exact in the large N limit, and to already hold for the Rice experiment Hulet () since k_{B}T_{F}/(\hbar\omega_{\alpha})\gtrsim 9.

The key ingredient of the local density approximation is the equation of state of the interacting homogeneous gas. At zero temperature, and in the zero-range limit detailed in the introduction, the equation of state is expected to be a universal function of the gas density and of the scattering length, independent of the microscopic details of the interaction. Furthermore, the spatial dependence of the scattering length in presence of the trapping potential, due to U_{b}(\mathbf{r})\neq 2\,U(\mathbf{r}), is expected to be negligible in the zero-range limit, as argued in the appendix B. Under these assumptions, the local-density approximation for the trapped gas leads to the universal form

 -\frac{dE}{d(1/a)}=\frac{\hbar^{2}k_{F}}{m}\,N\,\mathcal{F}\left(\frac{1}{k_{F% }a}\right), (15)

where \mathcal{F} is a dimensionless function, expressed in the appendix C in terms of the dimensionless universal function f giving the energy per particle \epsilon^{\rm hom} of the homogeneous interacting gas:

 \epsilon^{\rm hom}=\epsilon_{0}^{\rm hom}f\left(\frac{1}{k_{F}^{\rm hom}a}% \right). (16)

Here \epsilon_{0}^{\rm hom} and k_{F}^{\rm hom} are the energy per particle and the Fermi wavevector of the ideal Fermi gas with the same total density n as the interacting gas:

 \displaystyle k_{F}^{\rm hom} \displaystyle= \displaystyle(3\pi^{2}n)^{1/3} (17) \displaystyle\epsilon_{0}^{\rm hom} \displaystyle= \displaystyle\frac{3}{10}\frac{\hbar^{2}\left(k_{F}^{\rm hom}\right)^{2}}{m}. (18)

In the BEC regime 0<k_{F}a\ll 1 we have LeyronasLHY (); HuangYang (); LeeYang ()

 \displaystyle f\left(\frac{1}{k_{F}^{\rm hom}a}\right) \displaystyle= \displaystyle-\frac{5}{3(k_{F}^{\rm hom}a)^{2}} \displaystyle+ \displaystyle\frac{5}{18\pi}k_{F}^{\rm hom}a_{d}\left[1+\frac{128}{15\sqrt{6% \pi^{3}}}(k_{F}^{\rm hom}a_{d})^{3/2}\right]+\ldots

where the first term comes from the dimer binding energy \hbar^{2}/(ma^{2}), the second term is the mean field interaction energy among the dimers, proportional to the dimer-dimer scattering length a_{d}, and the last term is the (bosonic) Lee-Huang-Yang correction. The solution of the 4-fermion problem gives PetrovPRL (); Leyronas ()

 a_{d}\simeq 0.60a. (20)

The local-density approximation then leads to

 \mathcal{F}\left(\frac{1}{k_{F}a}\right)=\frac{1}{k_{F}a}+\frac{5^{2/5}}{2^{12% /5}\cdot 7}\left(\frac{a_{d}}{a}\right)^{2/5}(k_{F}a)^{7/5}+\ldots (21)

where we have omitted the Lee-Huang-Yang correction.

In the BCS limit 0<-k_{F}a\ll 1, we have HuangYang (); LeeYang (); Abrikosov (); Galitski ()

 \displaystyle f\left(\frac{1}{k_{F}^{\rm hom}a}\right) \displaystyle= \displaystyle 1+\frac{10}{9\pi}k_{F}^{\rm hom}a (22) \displaystyle+ \displaystyle\frac{4(11-2\ln 2)}{21\pi^{2}}(k_{F}^{\rm hom}a)^{2}+\ldots

where the first term is the ideal gas result, the second term is the Hartree-Fock mean field term, and the last term is the (fermionic) Lee-Huang-Yang correction. This leads to 444 For a small positive scattering length 0<k_{F}a\ll 1, Eqs.(22,23) also apply to the atomic gas state; since the interaction potentials considered here have a two-body bound state for a>0, this atomic gas state is only metastable LudoYvanToyModel (), the ground state being the BEC of dimers considered previously [cf. Eq.(21)].

 \displaystyle\mathcal{F}\left(\frac{1}{k_{F}a}\right)=\frac{512}{945\pi^{2}}(k% _{F}a)^{2}\Big{[}1+\\ \displaystyle\left(\frac{256}{35\pi^{2}}-\frac{63+189\ln 2}{1024}\right)k_{F}a% +\ldots\Big{]}. (23)

Around the unitary limit 1/(k_{F}a)\to 0 we use the expansion

 f\left(\frac{1}{k_{F}^{\rm hom}a}\right)=\xi-\frac{\zeta}{k_{F}^{\rm hom}a}+\ldots, (24)

which gives for the trapped gas:

 \mathcal{F}(0)=\frac{2^{7}}{5^{2}\cdot 7\cdot\pi}\frac{\zeta}{\xi^{1/4}}\simeq 0% .28 (25)

where we took the estimates

 \displaystyle\xi \displaystyle\simeq \displaystyle 0.41 (26) \displaystyle\zeta \displaystyle\simeq \displaystyle 0.95\ . (27)

The estimate for \xi is obtained by averaging the predictions of Giorgini () (\xi=0.42(1)), of Juillet () (\xi=0.449(9)), of Reddy () (\xi=0.42(1)) and of BulgacCrossover () (\xi=0.37(5)). The estimate for \zeta is obtained from the calculation of the short range pair correlation of Lobo (), see subsection 3.1, and is close to the value \zeta\simeq 1.0 that may be extracted directly from the values of f in Giorgini () by approximating the derivative by a two-point formula. 555The values (26,27) of \xi and \zeta can be inserted into the expressions BulgacModes (); CombescotLeyronasComment () for the derivatives of the frequencies of the hydrodynamic collective modes with respect to 1/(k_{F}a) taken at unitarity. The resulting slope is compatible with the experimental data for the radial mode GrimmModes (); ThomasModePRA (), while the agreement with the data for the axial mode GrimmModeVieux () is only marginal. The data from GrimmModes () agree very well for all positive values of a with the theoretical result AstrakharchikMode () obtained from the fixed-node Monte Carlo equation of state Giorgini (), while in ThomasModePRA (); GrimmModeVieux () finite-temperature effects may play a role GrimmModes (); AstrakharchikMode (); Urban ().

### 2.3 Comparison to Rice experiment

At Rice Hulet () the quantity Z=N_{b}/(N/2) was measured for a lithium gas, by resonantly exciting the closed-channel molecules with a laser that transfers them to another, short-lived molecular state; the molecule depletion reflects into a reduction of the atom number that can be measured.

We now compare the Rice results to our prediction (8). We first insert (9,10) into (8) which gives

 N_{b}=Nk_{F}R_{*}\mathcal{F}\left(\frac{1}{k_{F}a}\right)\left(1-\frac{a_{\rm bg% }}{a}\right)^{2} (28)

where the function \mathcal{F} is defined in (15). The length R_{*} depends on the interchannel coupling and is related to the width of the Feshbach resonance as:

 R_{*}=\frac{\hbar^{2}}{ma_{\rm bg}\mu_{b}\Delta B}. (29)

For the B_{0}\simeq 834G Feshbach resonance in {}^{6}{\rm Li}, we have \mu_{b}\simeq 2\mu_{B} where \mu_{B} is the Bohr magneton Koehler_review (), \Delta B=-300G, a_{\rm bg}=-1405a_{0} where a_{0}=0.0529177 nm is the Bohr radius InnsbruckLi (); this gives R_{*}=0.0269 nm. From the Rice data for N_{b}/N, using (28), we thus get values of \mathcal{F}, i.e. of dE/d(1/a), that we compare in Fig.1 to theoretical predictions, given in the local-density approximation either by the expansions (21,23) in the weakly interacting regime, or by (25) at the unitary limit, or by an interpolation of the Monte Carlo results of Giorgini (); Lobo () in the strongly interacting regime (see appendix D for details) 666 The error bars in Fig. 1 include the systematic theoretical uncertainty resulting from the leading order correction to (10) given in InnsbruckLi (), but are largely dominated by the experimental error bars from Hulet () because we did not include the three experimental points from Hulet () which are deep in the BEC regime [1/(k_{F}a)>5]; in this regime, two-body theories which go beyond our two-channel model give good agreement with the Rice measurements Hulet (); Koehler_review (); StoofZ (). . The agreement is satisfactory, the experimental points being scattered around the theoretical prediction. In particular, a possible shift \delta B_{0} of the resonance location due to the fact that U_{b}(\mathbf{r})\neq 2U(\mathbf{r}), see (9), an effect not explicitly included in Hulet () for the calculation of 1/(k_{F}a) as a function of B, and which would manifest itself in horizontal shifts of the experimental data with respect to theory, does not show up in Fig.1. For the considered broad resonance of lithium, this is not surprising, see appendix B.

Our approach even allows to predict the time evolution of the particle number, in presence of the depleting laser, under the assumption that the molecular depletion rate is so weak that the equilibrium relation (28) holds at all times. Neglecting the slow off-resonant free-bound photoassociation process Hulet (), we have

 \frac{dN}{dt}=-2N_{b}(t)\frac{\Omega^{2}}{\gamma} (30)

where \Omega^{2}/\gamma is the effective decay rate of the closed-channel molecule, the Rabi frequency \Omega and the spontaneous emission rate \gamma being defined in Hulet ().777Eq. (30) can be obtained from the following master equation for the density matrix \hat{\rho} of the gas: d\hat{\rho}/dt=(\Omega^{2}/\gamma)\int d^{3}r\left[\psi_{b}(\mathbf{r})\hat{% \rho}\psi_{b}^{\dagger}(\mathbf{r})-\frac{1}{2}\{\psi_{b}^{\dagger}(\mathbf{r}% )\psi_{b}(\mathbf{r}),\hat{\rho}\}\right]+[H_{2},\hat{\rho}]/(i\hbar). N in (30) stands for N_{\rm f}+2N_{b}, where N_{\rm f} is the number of fermions; thus N\simeq N_{\rm f} in the regime considered in this paper.

Simple expressions can be obtained in limiting cases: In the BEC limit, the dominant term of the expansion (21) gives

 N(t)=N(0)\,e^{-\Gamma_{0}t}; (31)

in the unitary limit, (25) gives

 N(t)=\frac{N(0)}{\left(1+\Gamma_{0}t/6\right)^{6}}; (32)

and in the BCS limit the dominant term of (23) gives

 N(t)=\frac{N(0)}{\left(1+\Gamma_{0}t/2\right)^{2}}. (33)

Here, \Gamma_{0}=-[dN(0)/dt]/N(0) is the inital atom loss rate, which can be expressed in terms of the initial atom number N(0) using Eqs. (28,30). We see that in general N(t) is not an exponential function of time. For the last point on the BCS side in Hulet (), the (unpublished) data for N(t) are better fitted by Eq.(33) than by an exponential; however the value of \Gamma_{0} resulting from this new fitting procedure does not differ significantly from the published value HuletPrivate ().

## 3 Related observables

In this section, contrarily to subsections 2.2 and 2.3, we return to the general case and we do not assume that the temperature is \ll T_{F} and that the numbers of atoms in each spin component are equal. For concreteness we take the lattice model (3), even if the reasonings may be generalized to a continuous space model and to a two-channel model.

### 3.1 Short range pair correlations

As was shown in Tan_nkREVTEX (); Tan_dEREVTEX () and recently rederived in Braaten (), in the zero interaction range limit, the derivative of the gas energy with respect to the inverse scattering length can be expressed in terms of the short range behavior of the pair distribution function of opposite spin particles:

 -\left(\frac{d\langle E\rangle}{d(1/a)}\right)_{S}=\frac{4\pi\hbar^{2}}{m}\int d% ^{3}\mathbf{R}\,\underset{r\to 0}{\tilde{\lim}}r^{2}g_{\uparrow\downarrow}(% \mathbf{R}+\mathbf{r}/2,\mathbf{R}-\mathbf{r}/2) (34)

where the pair distribution function is given by

 g_{\uparrow\downarrow}(\mathbf{r}_{1},\mathbf{r}_{2})=\langle\psi_{\uparrow}^{% \dagger}(\mathbf{r}_{1})\psi_{\downarrow}^{\dagger}(\mathbf{r}_{2})\psi_{% \downarrow}(\mathbf{r}_{2})\psi_{\uparrow}(\mathbf{r}_{1})\rangle. (35)

Since we are dealing here with a finite range interaction model, we have taken a zero-range limit in (34), introducing the notation

 \underset{r\to 0}{\tilde{\lim}}=\lim_{r\to 0}\lim_{b\to 0}. (36)

Note that the order of the limits is important. The limit in (34) is reached for a distance r much smaller than |a|, the mean distance between particles, and the thermal de Broglie wavelength \lambda, but still much larger than the interaction range b.

We thus see that a measurement of the pair distribution function of the atoms gives access, using (8), to the number of closed-channel molecules. Experimentally, the pair distribution function was measured via the photoassociation rate in a 1D Bose gas Gangardt (); Weiss (), or simply by dropping the gas on a detector with high enough spatio-temporal resolution Shimizu (); Aspect_bosons (); Aspect_fermions (). The possibility of measuring the pair distribution function in the BEC-BCS crossover was studied in Lobo (); LoboScaling (). In the unitary limit, by inserting into Eq.(34) the value of \tilde{\lim}\,r^{2}g_{\uparrow\downarrow} calculated in Lobo () with the fixed-node Monte Carlo technique, we get the value (27) for the parameter \zeta defined in (24).

The relation (34) is the three-dimensional version of the relation obtained in Lieb () for a one-dimensional Bose gas with contact interaction, using the Hellmann-Feynman theorem. We now show that the Hellmann-Feynman technique also provides a simple derivation of (34). Taking the derivative of the gas eigenergy E with respect to the effective coupling constant g=4\pi\hbar^{2}a/m, for a fixed lattice spacing b, the Hellmann-Feynman theorem gives

 \displaystyle\left(\frac{d\langle E\rangle}{dg}\right)_{S} \displaystyle= \displaystyle\frac{dg_{0}}{dg}\left\langle\sum_{\mathbf{r}}b^{3}\psi_{\uparrow% }^{\dagger}\psi_{\downarrow}^{\dagger}\psi_{\downarrow}\psi_{\uparrow}\right\rangle (37) \displaystyle= \displaystyle\frac{dg_{0}}{dg}\sum_{\mathbf{r}}b^{3}g_{\uparrow\downarrow}(% \mathbf{r},\mathbf{r}), (38)

where g_{0} is a function of g as given in (4), and where we used the relation \langle dE/dg\rangle=(d\langle E\rangle/dg)_{S} [cp. (11)]. In a gas, with an interaction range much smaller than the mean distance \sim 1/k_{F} between particles and than the thermal de Broglie wavelength \lambda, the pair distribution function is dominated by two-body physics if |\mathbf{r}_{1}-\mathbf{r}_{2}| is much smaller than 1/k_{F} and \lambda (and also much smaller than a in the BEC regime), a property contained in the Yvon Yvon () and Waldmann-Snider Snider () ansatz of kinetic theory, which also appeared in the context of quantum gases Lobo (); Huang_livre (); Holzmann (); Leggett_IHP ():

 g_{\uparrow\downarrow}(\mathbf{r}_{1},\mathbf{r}_{2})\simeq C\left(\frac{% \mathbf{r}_{1}+\mathbf{r}_{2}}{2}\right)\,|\phi(\mathbf{r}_{1}-\mathbf{r}_{2})% |^{2} (39)

where \phi is the two-body zero-energy free space scattering state for two atoms, here in the spin singlet state. With the lattice model we find 888 An elegant derivation is to enclose two atoms in a fictitious cubic box with periodic boundary conditions of volume L^{3}, with L\gg|a|,b. To leading order in a/L, the wavefunction is \phi_{\rm box}(\mathbf{r}_{1}-\mathbf{r}_{2})=\phi(\mathbf{r}_{1}-\mathbf{r}_{% 2})/L^{3/2} with an energy given by the “mean-field” shift E_{\rm box}=g/L^{3}. The Hellmann-Feynman theorem (37) then gives the result (40). An alternative derivation is to directly calculate \phi in Fourier space [see (47)] which gives \phi(\mathbf{r})=1-g_{0}\phi(\mathbf{0})\int_{[-\pi/b,\pi/b]^{3}}[d^{3}\mathbf% {k}/(2\pi)^{3}]\exp(i\mathbf{k}\cdot\mathbf{r})/(2\epsilon_{\mathbf{k}}); this, combined with (4), gives (40).

 |\phi(\mathbf{0})|^{2}=\frac{dg}{dg_{0}}, (40)

so that Eq.(38) simplifies to

 \left(\frac{d\langle E\rangle}{dg}\right)_{S}\simeq\sum_{\mathbf{r}}b^{3}C(% \mathbf{r}). (41)

On the other hand, at interparticle distances much larger than b, the zero-energy scattering wavefunction \phi behaves as 1-a/|\mathbf{r}_{1}-\mathbf{r}_{2}| so that

 \underset{r\to 0}{\tilde{\lim}}r^{2}g_{\uparrow\downarrow}(\mathbf{R}+\mathbf{% r}/2,\mathbf{R}-\mathbf{r}/2)=a^{2}C(\mathbf{R}). (42)

### 3.2 Tail of the momentum distribution

As was shown in Tan_dEREVTEX () and recently rederived in Braaten (), in the zero-range limit, the momentum distribution of the gas has a large momentum tail scaling as 1/k^{4}, with a coefficient common to the two spin states and proportional to -dE/d(1/a):

 \underset{k\to+\infty}{\tilde{\lim}}k^{4}n_{\sigma}(k)=-\frac{4\pi m}{\hbar^{2% }}\left(\frac{d\langle E\rangle}{d(1/a)}\right)_{\!S} (43)

where

 \underset{k\to+\infty}{\tilde{\lim}}=\underset{k\to+\infty}{\lim}\ \underset{b% \to 0}{\lim}, (44)

\sigma=\uparrow or \downarrow and the momentum distribution is normalized as

 \int\frac{d^{3}k}{(2\pi)^{3}}n_{\sigma}(\mathbf{k})=\langle N_{\sigma}\rangle, (45)

where \langle N_{\sigma}\rangle is the mean number of particle in the spin state \sigma. When combined with (8), this relation reveals that the number of closed-channel molecules may be deduced from a measurement of the momentum distribution of the atoms.

The momentum distribution in a two-component Fermi gas was actually already measured JILAnk (); TarruellVarenna (); JILA_T_finie (), simply by switching off the interaction and measuring the cloud density after a ballistic expansion. Two practical problems however arise, the effect of the non-zero switch-off time of the interaction JILAnk (); Holland_nk () and the signal to noise ratio in the far tails of the distribution; we have therefore not been able to extract the value of dE/d(1/a) from the published data for n_{\sigma}(k). The value of dE/d(1/a) has not been extracted either from the fixed-node Monte Carlo calculations of n_{\sigma}(k) performed in giorgini_nk ().

The relation (43) is a three-dimensional version of the relation derived for a one-dimensional Bose gas in Olshanii (). We now present a very simple rederivation of (43).

In a gas with short range interactions, i.e. b\ll 1/k_{F} and b\ll\lambda (and b\ll a on the BEC side), the tail of the momentum distribution is dominated by binary collisions in the spin singlet channel between a particle with a large momentum \mathbf{k} and a particle with a momentum \simeq-\mathbf{k}, so that at large k

 n_{\sigma}(\mathbf{k})\simeq\mathcal{B}|\tilde{\phi}(\mathbf{k})|^{2} (46)

where \mathcal{B} is momentum-independent and \tilde{\phi}(\mathbf{k}) is the Fourier transform \sum_{\mathbf{r}}b^{3}e^{-i\mathbf{k}\cdot\mathbf{r}}\phi(\mathbf{r}) of the zero-energy two-body free space scattering state \phi(\mathbf{r}). For the contact interaction in the lattice model the two-particle Schrödinger’s equation takes the form boite ()

 2\epsilon_{\mathbf{k}}\tilde{\phi}(\mathbf{k})+g_{0}\phi(\mathbf{r}=\mathbf{0}% )=0, (47)

where \epsilon_{\mathbf{k}} is the free wave dispersion relation in the lattice model, so that \tilde{\phi}(\mathbf{k}) drops as 1/\epsilon_{\mathbf{k}} for large k, and the associated kinetic energy diverges for b\to 0. Eq.(46) implies that the kinetic energy of the gas is proportional to the kinetic energy of the zero-energy two-body scattering state in the zero-range limit:

 E_{\rm kin}[{\rm gas}]\simeq\mathcal{B}\cdot E_{\rm kin}[\phi]. (48)

Since the total energies have a finite limit for b\to 0, we can replace the kinetic energies by the opposite of the interaction energies in the above formula:

 E_{\rm int}[{\rm gas}]\simeq\mathcal{B}\cdot E_{\rm int}[\phi]. (49)

The interaction energy of the gas is g_{0}\left\langle\sum_{\mathbf{r}}b^{3}\psi_{\uparrow}^{\dagger}\psi_{% \downarrow}^{\dagger}\psi_{\downarrow}\psi_{\uparrow}\right\rangle, equal to (d\langle E\rangle/dg)_{S}\,g_{0}(dg/dg_{0}) according to (37). The interaction energy of \phi is simply g_{0}|\phi(\mathbf{0})|^{2}, equal to g_{0}(dg/dg_{0}) according to (40). The common factor g_{0}(dg/dg_{0}) simplifies and

 \left(\frac{d\langle E\rangle}{dg}\right)_{S}=\mathcal{B}. (50)

It remains to relate the coefficient \mathcal{B} to \tilde{\lim}_{k\to+\infty}k^{4}n_{\sigma}(k), using (46). In the zero-range limit, the scattering state \phi(\mathbf{r}) tends to 1-a/r, so that, by Fourier transform, \tilde{\phi}(\mathbf{k}) tends to (2\pi)^{3}\delta(\mathbf{k})-4\pi a/k^{2}. We then get (43).

## 4 Conclusion

We have calculated the number N_{b} of closed-channel molecules in the whole BEC-BCS crossover for a two-component Fermi gas. Our result is in satisfactory agreement with the experimental results from Rice Hulet (). The expression (8) for N_{b} is proportional to the quantity -dE/d(1/a) where E is the energy of the gas; this quantity is universal in the zero interaction range limit. At unitarity we find N_{b}/N=k_{F}R_{*}\mathcal{F}(0) where \mathcal{F}(0) is a universal constant given in (25) and the length R_{*} is related to the resonance parameters by (29); in the zero range and zero effective range limit where k_{F}b\ll 1 and k_{F}|r_{e}|\ll 1 we have k_{F}R_{*}\ll 1 [see Eq.(65)], we thus confirm that N_{b}\ll N i.e. the gas mainly populates the open channel. In the BCS limit, -dE/d(1/a) and thus N_{b} are determined by the Hartree-Fock mean field energy, and do not tend exponentially to zero with 1/(k_{F}|a|), contrarily to the prediction of usual BCS theory Hulet (). This is also related to the fact that the BCS theory does not predict the correct short-distance pair correlations Iacopo ().

The quantity -dE/d(1/a) is related to the short-distance behavior of the pair correlation function and to the large-k tail of the momentum distribution, as discovered by Tan Tan_nkREVTEX (); Tan_dEREVTEX () and rederived in Braaten () and in this paper. The quantity -dE/d(1/a) is also proportional to the average radiofrequency shift Zwerger (); Baym (). The general idea that various short-range two-body quantities such as the number of closed-channel molecules, the interaction energy of the gas and the average radiofrequency shift are the product of the value of the considered observable in the zero-energy two-body scattering state and a single universal many-body quantity was given by Leggett Leggett_IHP (). Our work shows that the Rice experiment Hulet () is the first direct measurement of this fundamental quantity in the BEC-BCS crossover.

It would be interesting to study the effect of a finite temperature. Experimentally, this could be achieved by measuring the number of closed-channel molecules together with the density profile of the trapped gas, the released energy Thomas1 (); SalomonCrossover (), or the entropy JinPotentialEnergy (); thomas_entropie (). Theoretically, Eq.(12) can be used when more Monte Carlo results at finite temperature will become available. A useful tool in this context is the virial theorem: at unitarity it directly gives the total energy in terms of the density profile ThomasVirielExp (); thomas_entropie (); WernerSym (); ThomasVirielTheorie () while for a finite a it also involves [d\langle E\rangle/d(1/a)]_{S} TanViriel (); WernerViriel ().

###### Acknowledgements.
We thank A. Bulgac, F. Chevy, M. Colomé-Tatché, R. Hulet, M. Jona-Lasinio, F. Laloë, P. Naidon, S. Nascimbène, G. Partridge, L. Pricoupenko and C. Salomon for very useful discussions and comments, and G. Partridge et al. for providing us with the data of Hulet () in numerical form. Our group is a member of IFRAF.

Note: We became aware of the related work BraatenLong () while completing this paper, and of Zhang_Leggett () while revising it.

## Appendix A Two-body analysis of the two-channel model in free space

### A.1 Scattering amplitude

In this appendix we briefly discuss the two-body scattering properties for our model Hamiltonian H_{2} in free space, that is in the absence of trapping potential (U_{b}(\mathbf{r})\equiv 0,U(\mathbf{r})\equiv 0). The scattering problem may be considered in the center of mass frame so that we take the two-body scattering state

 |\psi\rangle=\beta b_{\mathbf{0}}^{\dagger}|0\rangle+\int\frac{d^{3}k}{(2\pi)^% {3}}A(\mathbf{k})a_{\mathbf{k}\uparrow}^{\dagger}a_{\mathbf{k}\downarrow}^{% \dagger}|0\rangle, (51)

with the unknown amplitude \beta in the closed channel and the unknown momentum dependent amplitude A(\mathbf{k}) in the open channel. The annihilation operators for the bosonic molecules, b_{\mathbf{k}}, and for the fermionic atoms a_{\mathbf{k}\sigma}, \sigma=\uparrow or \downarrow, are the Fourier transforms of the corresponding field operators \psi_{b}(\mathbf{r}) and \psi_{\sigma}(\mathbf{r}), where the Fourier transform of a function f(\mathbf{r}) is \tilde{f}(\mathbf{k})=\int d^{3}r\exp(-i\mathbf{k}\cdot\mathbf{r})f(\mathbf{r}). Inserting the state (51) in Schrödinger’s equation (H_{2}-E)|\psi\rangle=0 and projecting respectively on the molecular subspace and the atomic subspace leads to

 \displaystyle(E_{b}-E)\beta+\Lambda\gamma=0 (52) \displaystyle\left(\frac{\hbar^{2}k^{2}}{m}-E\right)A(\mathbf{k})+\tilde{\chi}% (\mathbf{k})(\beta\Lambda+g_{0}\gamma)=0 (53)

where we used the fact that the function \chi(\mathbf{r}) is real and even, and where we introduced the amplitude

 \gamma=\int\frac{d^{3}k}{(2\pi)^{3}}\tilde{\chi}(\mathbf{k})A(\mathbf{k}). (54)

We now specialize to a state corresponding to the scattering of a spin \uparrow atom with incoming wavevector \mathbf{k}_{0} onto a spin \downarrow atom with incoming wavevector -\mathbf{k}_{0}, so that E=\hbar^{2}k_{0}^{2}/m is the energy of the incoming state. Eliminating \beta in terms of \gamma with (52) and using (53) leads to

 A(\mathbf{k})=(2\pi)^{3}\delta(\mathbf{k}-\mathbf{k_{0}})+\frac{\gamma\tilde{% \chi}(\mathbf{k})[g_{0}+\Lambda^{2}/(E-E_{b})]}{E+i0^{+}-\hbar^{2}k^{2}/m} (55)

where we have singled out the incoming state, and the scattered component is guaranteed to correspond to a pure outgoing wave in real space by the introduction of +i0^{+} in the denominator. Injecting (55) in (54) gives a closed linear equation for \gamma. From the known large r expression for the kinetic energy Green’s function we extract the scattering amplitude

 f_{k_{0}}=\frac{-m\tilde{\chi}(\mathbf{k}_{0})^{2}/(4\pi\hbar^{2})}{[g_{0}+% \frac{\Lambda^{2}}{\frac{\hbar^{2}k_{0}^{2}}{m}-E_{b}}]^{-1}-\int\frac{d^{3}k}% {(2\pi)^{3}}\frac{\tilde{\chi}^{2}(\mathbf{k})}{\frac{\hbar^{2}k_{0}^{2}}{m}+i% 0^{+}-\frac{\hbar^{2}k^{2}}{m}}}. (56)

Since we are considering a s-wave Feshbach resonance, the function \chi(\mathbf{r}) is rotationally invariant, a property that we have used in (56).

We now perform explicit calculations of various quantities characterizing the low energy scattering, such as the scattering length and the effective range. To this end, it is convenient as in YvanVarenna () to perform the specific choice

 \tilde{\chi}(\mathbf{k})=\exp(-k^{2}b^{2}/2) (57)

and to extend (56) to negative energies, setting k_{0}=iq,q>0 and E=-\hbar^{2}q^{2}/m. One then obtains an explicit expression

 -\frac{1}{f_{iq}}=e^{-q^{2}b^{2}}\left[\frac{4\pi\hbar^{2}/m}{g_{0}-\frac{% \Lambda^{2}}{E_{b}+\hbar^{2}q^{2}/m}}+\frac{1}{b\sqrt{\pi}}\right]+q\,\mathrm{% erf}\,(qb)-q (58)

where the last term -q is imposed by the optical theorem and \mathrm{erf} is the error function.

The inverse of the scattering length is obtained by taking q\to 0 in (58):

 \frac{1}{a}=\frac{4\pi\hbar^{2}/m}{g_{0}-\Lambda^{2}/E_{b}}+\frac{1}{b\sqrt{% \pi}}. (59)

The large E_{b} limit of this expression is the inverse of the background scattering length:

 \frac{1}{a_{\rm bg}}=\frac{4\pi\hbar^{2}/m}{g_{0}}+\frac{1}{b\sqrt{\pi}}. (60)

The location of the Feshbach resonance corresponds to an energy E_{b}^{0} of the closed-channel molecule such that 1/a=0, so that

 E_{b}^{0}=\frac{\Lambda^{2}}{g_{0}+4\pi^{3/2}\hbar^{2}b/m}. (61)

Supposing that the energy E_{b} of the closed-channel molecule varies linearly with the magnetic field B,

 E_{b}(B)=E_{b}^{0}+\mu_{b}(B-B_{0}), (62)

and introducing the (possibly negative) resonance width \Delta B such that

 \mu_{b}\Delta B=\frac{\Lambda^{2}}{g_{0}}-E_{b}^{0}, (63)

we get from (59) exactly the B-field dependence (10) for the scattering length a in free space.

We now consider the effective range r_{e} defined by the low-q expansion -1/f_{iq}=a^{-1}-q+\frac{1}{2}q^{2}r_{e}+O(q^{4}). Its general expression for our model is

 r_{e}=\frac{-8\pi\hbar^{4}\Lambda^{2}/m^{2}}{(g_{0}E_{b}-\Lambda^{2})^{2}}+% \frac{4b}{\sqrt{\pi}}-\frac{2b^{2}}{a}. (64)

Of particular interest is the value of the effective range right on resonance,

 r_{e}^{\rm res}=-2R_{*}+\frac{4b}{\sqrt{\pi}}, (65)

where the length R_{*} is always non-negative,

 R_{*}=\left(\frac{\Lambda}{2\pi bE_{b}^{0}}\right)^{2}. (66)

This definition of R_{*} agrees with Eq. (29).

A very convenient rewriting of the scattering amplitude f_{iq} is obtained by taking as independent parameters, in addition to the van der Waals length b, the two lengths a_{\rm bg},R_{*} rather than the bare parameters g_{0},\Lambda of the Hamiltonian. The quantity E_{b}-\Lambda^{2}/g_{0} appearing in (58) is conveniently rewritten using the quantity

 \displaystyle Q^{2} \displaystyle\equiv \displaystyle\frac{m}{\hbar^{2}g_{0}}\left(g_{0}E_{b}-\Lambda^{2}\right) (67) \displaystyle= \displaystyle\frac{-1}{a_{\rm bg}R_{*}}\frac{1}{1-a_{\rm bg}/a} (68)

which can be positive or negative. The ratio g_{0}/\Lambda is also easily eliminated with the identity

 \frac{g_{0}^{2}}{4\pi\Lambda^{2}}=R_{*}a_{\rm bg}^{2}. (69)

Our final expression for the scattering amplitude is thus

 -\frac{1}{f_{iq}}=\frac{e^{-q^{2}b^{2}}}{a}\left[1-\left(1-a/a_{\rm bg}\right)% \frac{q^{2}}{q^{2}+Q^{2}}\right]+q\,\mathrm{erf}\,(qb)-q (70)

i.e.

 -\frac{1}{f_{k}}=\frac{e^{k^{2}b^{2}}}{a}\left[1-\left(1-a/a_{\rm bg}\right)% \frac{k^{2}}{k^{2}-Q^{2}}\right]-ik\,\mathrm{erf}\,(-ikb)+ik. (71)

We deduce the effective range in terms of a,a_{\rm bg},R_{*} and b:

 r_{e}=-2R_{*}(1-a_{\rm bg}/a)^{2}+\frac{4b}{\sqrt{\pi}}-\frac{2b^{2}}{a}. (72)

### A.2 Conditions for reaching the zero-range limit

In a degenerate gas, we expect that the zero-range limit is reached and the equation of state E(1/a) becomes model independent provided the conditions

 \displaystyle kb \displaystyle\ll \displaystyle 1 (73) \displaystyle-\frac{1}{f_{k}} \displaystyle\simeq \displaystyle\frac{1}{a}+ik (74)

hold for k\sim k_{F} and for k\sim 1/a on the BEC side YvanVarenna ().

We now assume that (73) holds, and we derive validity conditions for (74) in the case of our two-channel model. Following the usual procedure, we use the second order expansion

 -\frac{1}{f_{k}}=\frac{1}{a}+ik-\frac{1}{2}r_{e}k^{2}+\dots, (75)

which reduces to (74) if the effective range is small enough to have

 \frac{1}{2}|r_{e}|k^{2}\ll\left|\frac{1}{a}+ik\right|. (76)

Taking into account the condition (73) as well as the expression (72) of r_{e} for the two-channel model, and using |i+1/(ka)|\geq 1 and |1+ika|\geq 1, we find that the condition (76) becomes equivalent to

 k^{2}R_{*}\left(1-\frac{a_{\rm bg}}{a}\right)^{2}\left|\frac{1}{a}+ik\right|^{% -1}\ll 1. (77)

A second condition is obtained by imposing that the higher order terms \dots in (75) are indeed negligible compared to k^{2}|r_{e}|/2. For this to hold we need not only (73), but also k^{2}\ll|Q^{2}| i.e.

 k^{2}|a_{\rm bg}|R_{*}\left|1-\frac{a_{\rm bg}}{a}\right|\ll 1. (78)

We conclude that to reach the universal zero-range limit, it is sufficient to fulfill the three conditions (73,77,78). All three conditions are satisfied for the experimental data on {}^{6}{\rm Li} from Hulet () reproduced in Fig.1, where 1/k_{F}>250 nm and |a|>200 nm: using the values |a_{\rm bg}|=74 nm and R_{*}=0.027 nm (see section 2.3) and taking the estimate b\simeq 2.1 nm obtained by imposing that the effective range on resonance (65) should be equal to the value 4.7 nm obtained in the multi-channel calculation of Strinati (), we find that the left hand sides of Eqs.(73,77,78) are respectively smaller than 10^{-2}, 3\cdot 10^{-4} and 7\cdot 10^{-5}.

To be complete, let us point out that the conditions (77,78) are not stricly necessary since it is possible that the \dots in (75) are not negligible compared to k^{2}|r_{e}|/2 and (74) nevertheless holds. We have chosen these conditions because they simplify the discussion, in particular they automatically ensure that the analytic continuation f_{iq} has a pole at q\simeq 1/a so that the dimer energy is \simeq-\hbar^{2}/(ma^{2}). For this last property, in the regime b\ll a, one directly gets from (70) the truly necessary condition

 \left|(1-a/a_{\rm bg})^{-1}+(1-a_{\rm bg}/a)^{-2}a/R_{*}\right|\gg 1. (79)

Close to the Feshbach resonance, a\gg|a_{\rm bg}|, this reduces to the simpler condition

 a\gg R_{*}. (80)

## Appendix B Effect of a difference between the trapping potentials seen by a closed-channel molecule and by an open-channel pair of atoms

In this appendix, we consider the physical consequences of the fact that, in general, the trapping potential U_{b}(\mathbf{r}) seen by a closed-channel molecule is not simply the trapping potential 2\,U(\mathbf{r}) seen by two separated atoms in the open channel:

 U_{b}(\mathbf{r})\neq 2\,U(\mathbf{r}), (81)

and we derive conditions for the gas to behave approximately as if we had U_{b}(\mathbf{r})=2\,U(\mathbf{r}). Note that a significant deviation in (81) is not a priori excluded for the laser traps used e.g. in Hulet (); InnsbruckLi () for lithium, with a laser wavelength \lambda_{L}\simeq 1\,\mum: Contrarily to the atomic lightshift, the lightshift of the closed-channel molecule may be sensitive to the laser detuning from transitions to bound states of the 1^{1}\Sigma_{u}^{+}(A) Born-Oppenheimer interaction potential between a ground-state atom and an excited atom, considering the fact that the minimum of this excited-ground interaction potential is separated from the dissociation limit of the ground-ground interaction potential by only \simeq 0.6\,\mu{\rm m}^{-1}<1/\lambda_{L} in spectroscopic units  ChemPhys ().

### B.1 Position dependent scattering length and effective magnetic field

At the two-body level, (81) leads to an effective position dependent scattering length for two opposite spin atoms. To define such a local scattering length around point \mathbf{r}_{0}, we consider the “tangential” free space problem where the trapping potentials are replaced by the constant values U_{b}(\mathbf{r}_{0}) and U(\mathbf{r}_{0}) and we call a(\mathbf{r}_{0}) the resulting scattering length. We simply have to revisit the calculations of appendix A: In (52) one has to replace E_{b} by E_{b}+U_{b}(\mathbf{r}_{0}) to include the trapping energy shift of the closed-channel molecule; similarly in (53) one has to replace \frac{\hbar^{2}k^{2}}{m} by \frac{\hbar^{2}k^{2}}{m}+2U(\mathbf{r}_{0}) to include the trapping energy shift of the two atoms; finally, for atoms with incoming wavevectors \pm\mathbf{k}_{0}, the scattering state eigenenergy is now E=\frac{\hbar^{2}k_{0}^{2}}{m}+2U(\mathbf{r}_{0}). The net effect on the scattering amplitude f_{k_{0}} in (56) is to apply the substitution

 E_{b}\longrightarrow E_{b}+U_{b}(\mathbf{r}_{0})-2U(\mathbf{r}_{0}), (82)

so that a(\mathbf{r}_{0}) is simply obtained by applying (82) to (59).

A powerful interpretation of this substitution procedure may be given for a magnetic field independent \mu_{b}: We can rewrite the term E_{b}(B)+U_{b}(\mathbf{r}) appearing in (1) as E_{b}[B_{\rm eff}(\mathbf{r})]+2\,U(\mathbf{r}) provided we set

 B_{\rm eff}({\bf r})=B+\frac{U_{b}(\mathbf{r})-2\,U(\mathbf{r})}{\mu_{b}}. (83)

Thus we can replace U_{b}(\mathbf{r}) by 2\,U(\mathbf{r}) in (1) provided we also replace the homogeneous magnetic field B in (1) by the inhomogeneous effective magnetic field B_{\rm eff}({\bf r}). One then immediately obtains the position dependence of the scattering length,

 a(\mathbf{r})=a_{\rm fs}[B_{\rm eff}(\mathbf{r})], (84)

the function a_{\rm fs}(B) being given in (10).

This effective magnetic field picture is also very efficient at the many-body level. In the grand canonical point of view, one adds the term -\mu\hat{N} to the Hamiltonian H_{2}, where \mu is the chemical potential, and the particle number operator \hat{N} is defined as \hat{N}=\hat{N}_{f}+2\hat{N}_{b}, \hat{N}_{f} giving the number of fermions and \hat{N}_{b} the number of bosons, since this is the quantity which commutes with H_{2}. In the local density approximation, in the general case U_{b}\neq 2U, we thus see that the gas behaves locally like a free space, homogeneous gas of chemical potential \mu-U(\mathbf{r}), in presence of the locally homogeneous magnetic field B_{\rm eff}({\bf r}).

### B.2 Shift of the resonance location

We now consider the case of an optical trap such as the one at Rice Hulet (): U(\mathbf{r}) and U_{b}(\mathbf{r}) are the lightshift potentials experienced by an atom and by a closed-channel molecule. Being both proportional to the laser intensity, they are mutually proportional:

 U_{b}(\mathbf{r})=2(1+\eta)U(\mathbf{r}) (85)

where the dimensionless parameter \eta is a convenient measure of the relative deviation between U_{b} and 2U. At large \mathbf{r} the laser intensity vanishes, so do U_{b} and U.

In practice, the gas is trapped around the intensity maximum of the laser field, located in \mathbf{r}=\mathbf{0}, and the extension of the cloud along the axial (resp. radial) direction is small compared to the Rayleigh length (resp. the waist) of the laser beam; thus, to leading order, U(\mathbf{r})\simeq U(\mathbf{0}), we can temporarily forget about the spatial variation of B_{\rm eff} and take B_{\rm eff}\simeq B-\delta B_{0}, where

 \delta B_{0}=\frac{-U_{b}(\mathbf{0})+2\,U(\mathbf{0})}{\mu_{b}}=-\frac{2\eta U% (\mathbf{0})}{\mu_{b}}. (86)

Since the free space scattering length a_{\rm fs}(B) only depends on the difference B-B_{0}, replacing B by B_{\rm eff} is equivalent to shifting the position B_{0} of the Feshbach resonance by the quantity \delta B_{0}. Introducing the length R_{*} defined in (29), the ratio of this shift to the width \Delta B of the Feshbach resonance may be written as

 \frac{\delta B_{0}}{\Delta B}=-\frac{2mU(\mathbf{0})}{\hbar^{2}}a_{\rm bg}R_{*% }\eta. (87)

In the Rice experiment, the trap depth -U(\mathbf{0}) is at most k_{B}\times 2.71~{}\mu{\rm K} HuletPrivateData (), so that we have

 \frac{|\delta B_{0}|}{|\Delta B|}<\frac{|\eta|}{7500} (88)

where we have used the values of a_{\rm bg} and R_{*} given below (29) and the mass of lithium. Except if \eta takes truly large values for the laser wavelength \lambda_{L}=1064 nm used in Hulet (), we expect only a small shift of the resonance location due to the laser trap, relative to the resonance width.

In Hulet () the value of the scattering length as a function of B was calculated essentially from the free space formula a_{\rm fs}(B) ignoring the shift \delta B_{0} due to the trap, and taking for the resonance location B_{0} the value measured in InnsbruckLi () in a laser trap with different laser wavelength \lambda_{L}=1030 nm and intensity. On the contrary, our definition of a as a function of B includes the shift \delta B_{0}, a=a_{\rm fs}(B-\delta B_{0}), as already defined in (9). This may introduce a systematic discrepancy between our theory curves and the experimental data in Fig.1, in the form of shifts of the abscissas of the experimental points. We thus define the horizontal shift

 \displaystyle\delta\left(\frac{1}{k_{F}a}\right) \displaystyle\equiv \displaystyle\frac{1}{k_{F}a_{\rm fs}(B-\delta B_{0})}-\frac{1}{k_{F}a_{\rm fs% }(B)} (89) \displaystyle= \displaystyle\frac{\delta B_{0}}{\Delta B}\,\frac{1}{k_{F}a_{\rm bg}}\frac{(1-% a_{\rm bg}/a)^{2}}{1-(1-a_{\rm bg}/a)\delta B_{0}/\Delta B}. (90)

To leading order in \eta (that is for |\delta B_{0}/\Delta B|\ll 1), and using the fact that the data from Hulet () reproduced in Fig. 1 is sufficiently close to resonance to have |1-a_{\rm bg}/a|<1.4, this simplifies to:

 \delta\left(\frac{1}{k_{F}a}\right)\simeq\frac{\delta B_{0}}{\Delta B}\,\frac{% 1}{k_{F}a_{\rm bg}}(1-a_{\rm bg}/a)^{2}. (91)

Using the parameters of the Rice experiment HuletPrivateData () and Eq.(87) this gives

 \left|\delta\left(\frac{1}{k_{F}a}\right)\right|<\frac{|\eta|}{1000}. (92)

Even without knowing \eta, we can thus resonably expect that this shift is negligible on the scale of Fig. 1, i.e. that

 \left|\delta\left(\frac{1}{k_{F}a}\right)\right|\ll 1. (93)

This is also consistent with the fact that no systematic horizontal shift is apparent in Fig.1 between the experimental data and the theory curves.

### B.3 Effect of the position dependence of B_{\rm eff}(\mathbf{r})

We now discuss the effect on the many-body properties of the gas of the inhomogeneity of B_{\rm eff}(\mathbf{r}), that is of a spatial dependence of the scattering length a(\mathbf{r}), see (84), now taking into account the \mathbf{r} dependence of U(\mathbf{r}). We rewrite (83) as

 B_{\rm eff}(\mathbf{r})=B+\frac{2\eta U(\mathbf{r})}{\mu_{b}}=B-\delta B_{0}+% \frac{2\eta[U(\mathbf{r})-U(\mathbf{0})]}{\mu_{b}}. (94)

When inserted in (84) this leads to the following position dependence for the scattering length:

 \frac{1}{a(\mathbf{r})-a_{\rm bg}}=\frac{1}{a-a_{\rm bg}}-\eta\frac{2mR_{*}}{% \hbar^{2}}[U(\mathbf{r})-U(\mathbf{0})], (95)

where we recall that a is the scattering length in the trap center. Right on resonance, |a|=+\infty, R_{*} is one of the contributions to the interaction effective range, see (65), so that the zero-range limit implies k_{F}|R_{*}|\ll 1. Since the typical value of U(\mathbf{r})-U(\mathbf{0}) over the cloud size is of the order of \hbar^{2}k_{F}^{2}/2m, we see that the small parameter k_{F}R_{*} appears in (95) so that we may expect that the position dependence of a(\mathbf{r}) is negligible in the zero-range limit (for a fixed value of \eta).

Let us investigate this expectation more quantitatively and beyond the unitary limit, however still in a limiting case with simplifying assumptions. First, we assume that we are sufficiently close to resonance to have |a(\mathbf{r})|\gg|a_{\rm bg}| over the atomic cloud size. This implies that the excursion of B_{\rm eff}(\mathbf{r}) over the cloud size is small as compared to |\Delta B|; since the excursion of U(\mathbf{r}) over the laser trap size is at most |U(\mathbf{0})|, we note that this condition is automatically satisfied for |\delta B_{0}|\ll|\Delta B|. Second, we shall assume that the gas may be treated in the zero temperature local density approximation. Third, we use the quadratic approximation (13) for U(\mathbf{r}), further supposing for simplicity that the harmonic trap is isotropic, a case to which one can always reduce within the local density approximation by a rescaling of the coordinates, see appendix C. Equation (95) then reduces to

 \frac{1}{a(r)}\simeq\frac{1}{a}-\eta\frac{mR_{*}}{\hbar^{2}}m\omega^{2}r^{2}, (96)

and the gas density profile n(\mathbf{r}) is given by

 \mu-U(\mathbf{0})=\frac{1}{2}m\omega^{2}r^{2}+\mu^{\rm hom}[n(\mathbf{r}),% \frac{1}{a(\mathbf{r})}] (97)

where \mu^{\rm hom}[n,1/a] is the chemical potential of a free space homogeneous gas of density n and scattering length a.

Let us derive simple conditions to have a weak effect of the scattering length inhomogeneity on the density profile. Some of these conditions will be expressed in terms of k_{F} defined by (14).

In the BCS regime, a<0 and k_{F}|a|<1, we have \mu-U(\mathbf{0})\simeq k_{B}T_{F}=\hbar^{2}k_{F}^{2}/2m. Over the Thomas-Fermi radius R of the cloud, such that \frac{1}{2}m\omega^{2}R^{2}\simeq k_{B}T_{F}, we see that the relative change of 1/a(r) is negligible if

 |1-a/a(R)|\simeq(k_{F}|a|)(|\eta|k_{F}R_{*})\ll 1. (98)

This is satisfied in the zero-range limit, for fixed \eta, according to the condition (77) written for k=k_{F}.

In the unitary limit k_{F}|a|\gg 1, the previous reasoning fails, because requiring that the relative change of a(r) is small is too stringent. Taking for simplicity the limiting case |a|=+\infty, it is sufficient to check that the gas remains locally unitary except in a small region near the edge of the cloud. First, to zeroth order in \eta, we calculate the local Fermi wavevector k_{F}^{\rm hom}(\mathbf{r}) using the unitary equation of state \mu^{\rm hom}[n,1/a=0]=\xi\hbar^{2}(k_{F}^{\rm hom})^{2}/2m, where k_{F}^{\rm hom} is defined in (17), and using the value of the chemical potential \mu-U(\mathbf{0})=\sqrt{\xi}k_{B}T_{F}. We obtain

 k_{F}^{\rm hom}(r)=\frac{k_{F}}{\xi^{1/4}}\left(1-\frac{r^{2}}{R^{2}}\right)^{% 1/2} (99)

where R is the unperturbed Thomas-Fermi radius such that m\omega^{2}R^{2}/2=\sqrt{\xi}k_{B}T_{F}. Then, we include the \eta position dependent term in 1/a(r) to calculate at which distance l from the edge of the cloud the gas starts leaving the unitary regime, that is

 k_{F}^{\rm hom}(R-l)|a(R-l)|=1. (100)

Assuming a priori l\ll R we obtain from (99)

 \frac{l}{R}\simeq\frac{\xi^{3/2}}{2}(\eta k_{F}R_{*})^{2}. (101)

We conclude a posteriori that l/R\ll 1 and that the spatial dependence of a(r) is negligible when

 |\eta|k_{F}R_{*}\ll 1. (102)

As in the BCS case, this is equivalent for a fixed |\eta| to the zero range condition (77) written here with k=k_{F} and |a|=+\infty.

Finally we turn to the BEC regime a>0 and k_{F}|a|<1. We then use the approximate equation of state

 \mu^{\rm hom}[n(r),1/a(r)]\simeq-\frac{\hbar^{2}}{2ma^{2}(r)}+\frac{\pi\hbar^{% 2}}{2m}\alpha a(r)n(r) (103)

where \alpha\simeq 0.60 is the proportionality factor in (20) between the dimer-dimer scattering length a_{d} and a. In particular, this requires a\gg R_{*}, see (80). Replacing 1/a(r) by (96) in the above expression, we find that the leading term in 1/a^{2}(r) gives rise to an external potential that sums up with the trapping potential to give

 \displaystyle U_{\rm eff}(r) \displaystyle\equiv \displaystyle U(r)-\frac{\hbar^{2}}{2ma^{2}(r)} (104) \displaystyle= \displaystyle U_{\rm eff}(0)+\frac{1}{2}\omega_{\rm eff}^{2}r^{2}-\frac{m(\eta R% _{*})^{2}}{2\hbar^{2}}(m\omega^{2}r^{2})^{2}. (105)

We see that the atomic oscillation frequency is renormalized into

 \omega_{\rm eff}=\left(1+\frac{2\eta R_{*}}{a}\right)^{1/2}\omega. (106)

To have a small effect of the scattering length inhomogeneity on the trapped gas, we first require \omega_{\rm eff}\simeq\omega so that

 |\eta|R_{*}/a\ll 1, (107)

a condition with no counterpart in the BCS regime and more stringent than (102) since k_{F}a<1 here. Again, we find that we recover, for a fixed |\eta|, the zero range condition (77) written here for k=1/a. Second we require that the r^{4} term in U_{\rm eff}(r) remains small as compared to the r^{2} one over the Thomas-Fermi radius R of the condensate of dimers, calculated here to zeroth order in \eta; we reach the second condition

 \frac{m(\eta R_{*})^{2}}{\hbar^{2}}m\omega^{2}R^{2}\simeq(\eta k_{F}R_{*})^{2}% \left(\frac{5}{64}k_{F}a_{d}\right)^{2/5}\ll 1. (108)

Since k_{F}a<1 here, this second condition is automatically ensured by the first condition (107). Third, we require that the relative variation of 1/a(r) is small over the Thomas-Fermi radius R of the condensate, to ensure that the spatially inhomogeneity of a(r) in the mean field term [\pi\hbar^{2}/(2m)]\alpha a(r)n(r) has a small impact on the density profile:

 |1-a/a(R)|\simeq(k_{F}a)(|\eta|k_{F}R_{*})\left(\frac{5}{64}k_{F}a_{d}\right)^% {2/5}\ll 1, (109)

which is also implied by (107).

To summarize, discussing perturbatively the effect of the position dependence of B_{\rm eff}(\mathbf{r}) on the gas density profile in the local density approximation, we find close to the Feshbach resonance (|a|\gg|a_{\rm bg}|, and also a\gg R_{*} in the BEC regime) that the effect is small under the condition (98) in the BCS regime, under the condition (102) in the unitary regime, and under the condition (107) in the BEC regime: In all three cases, this is the zero range condition (77), written respectively for k=k_{F},k=k_{F} and k=1/a, and multiplied by the parameter |\eta|. If |\eta| is not too large, neglecting the r dependence of the scattering length is thus automatically justified in the zero-range limit.

Interestingly, in the unitary limit, it is possible to relate \eta k_{F}R_{*} to the quantity \delta(1/k_{F}a) introduced in the previous subsection:

 |\eta|k_{F}R_{*}\xi^{1/2}\simeq\left|\delta\left(\frac{1}{k_{F}a}\right)\frac{% \mu-U(\mathbf{0})}{U(\mathbf{0})}\right|. (110)

Experimentally, the gas is confined at the bottom of the laser trap so that [\mu-U(\mathbf{0})]/|U(\mathbf{0})|\ll 1 (its value near unitarity at Rice is 0.15 HuletPrivateData ()). Together with the assumption (93) this implies that the right-hand-side of (110) is \lll 1, and that the position dependence of a(\mathbf{r}) has a negligible effect on the gas density.

A quantitative extension of the present qualitative discussion is to evaluate the density change \delta n(r) of the gas to first order in \eta. To this end, one simply expands (97) to first order in \eta around the solution n_{0}(r) with \eta=0, writing to first order n(r)=n_{0}(r)+\delta n(r) and \mu=\mu_{0}+\delta\mu:

 \delta n(\mathbf{r})=\frac{\delta\mu-\delta\left(\frac{1}{a}\right)(r)\partial% _{1/a}\mu^{\rm hom}[n_{0}(r),1/a]}{\partial_{n}\mu^{\rm hom}[n_{0}(r),1/a]} (111)

where we have introduced the variation of 1/a(r) around \eta=0 to first order in \eta, not supposing |a(r)|\gg|a_{\rm bg}| in (95):

 \delta\left(\frac{1}{a}\right)(r)=-\left(1-\frac{a_{\rm bg}}{a}\right)^{2}\eta% \frac{mR_{*}}{\hbar^{2}}m\omega^{2}r^{2}. (112)

\delta\mu is then determined by the condition of a fixed number of particles \int_{r<R_{0}}d^{3}r\,\delta n(r)=0, where R_{0} is the Thomas-Fermi radius of the gas for \eta=0. In this way, taking the variation of (130) to first order in \eta around \eta=0, one can calculate the effect of a non-uniformity of a(r) on dE/d(1/a) to first order in \eta. We give here the result for the unitary gas:

 \displaystyle\frac{\delta\left[dE/d(1/a)\right]}{dE/d(1/a)|_{\eta=0}}=\eta k_{% F}R_{*}\left[\left(\frac{64}{25\pi}-\frac{63\pi}{256}\right)\xi^{-1/4}\zeta% \right.\\ \displaystyle\left.+\frac{315\pi}{2048}\,\xi^{3/4}\zeta^{-1}f^{(2)}(0)+\frac{2% }{3}\,\xi^{1/2}k_{F}a_{\rm bg}\right], (113)

where we recall that \xi=f(0), \zeta=-f^{\prime}(0) , see (24), and f^{(2)}(0) is the second order derivative of f(x) with respect to x in x=0. Note that, in this first order variation with respect to \eta, the trapping potential U(\mathbf{r}) is kept fixed so that k_{F} is also fixed.

The calculation (113) shows that, at unitarity, a measurement of the variation of the number of closed-channel molecules N_{b} as a function of U_{b}(\mathbf{r})/U(\mathbf{r}) may give access to f^{(2)}(0). It also shows that, in the unitary limit, the condition to neglect a_{\rm bg} with respect to a(\mathbf{r}) as done in the qualitative reasoning around (102), is k_{F}|a_{\rm bg}|\ll 1.

### B.4 Metastability issues for a position dependent scattering length

We now discuss whether it is energetically possible for a pair of atoms to form a dimer which could escape from the trap by tunneling. The energy variation in such a process is

 \Delta E=-\frac{\hbar^{2}}{ma^{2}(\infty)}-2\mu (114)

where \mu is the chemical potential of the trapped gas and a(\infty) is the effective scattering length at infinity, which we assume to be positive. We discuss under which conditions this process is energetically forbidden, i.e. \Delta E>0.

We assume that the gas has the usual Thomas-Fermi density profile and is weakly affected by the position dependence of B_{\rm eff}(\mathbf{r}), see subection B.3. We distinguish between two regimes for the scattering length a in the trap center: the negative-a regime where 1/a\leq 0, and the BEC regime where 0<k_{F}a<1. Under the usual assumption that the gas is confined in a small region around the bottom of the laser trap we have \mu\simeq U(\mathbf{0}) in the negative-a regime and \mu\simeq U(\mathbf{0})-\hbar^{2}/(2ma^{2}) in the BEC regime. We also make the simplifying assumption a(\infty)\gg|a_{\rm bg}| so that (95) reduces to

 \frac{1}{a(\infty)}\simeq\frac{1}{a}+2\eta\frac{m}{\hbar^{2}}R_{*}U(\mathbf{0}). (115)

In the BEC regime one easily deduces

 \frac{\Delta E}{2|U(\mathbf{0})|}\simeq 1+\mathcal{A}-\mathcal{B} (116)

where

 \displaystyle\mathcal{A} \displaystyle\equiv \displaystyle 2\eta\frac{R_{*}}{a} (117) \displaystyle\mathcal{B} \displaystyle\equiv \displaystyle 2\eta^{2}\frac{mR_{*}^{2}}{\hbar^{2}}|U(\mathbf{0})|; (118)

since we already assumed (107) we automatically have |\mathcal{A}|\ll 1. In the negative-a regime, Eq.(115) immediately gives an upper bound on 1/a(\infty), from which we get

 \frac{\Delta E}{2|U(\mathbf{0})|}\gtrsim 1-\mathcal{B}. (119)

In both regimes we conclude that under the condition

 |\mathcal{B}|\ll 1 (120)

we have \Delta E\gtrsim 2|U(\mathbf{0})| so that the process is energetically forbidden.

Finally we note that under the assumption |a|\gg|a_{\rm bg}| we have from (90,87)

 \delta\left(\frac{1}{k_{F}a}\right)\simeq\eta k_{F}R_{*}\frac{|U(\mathbf{0})|}% {k_{B}T_{F}} (121)

so that

 \mathcal{B}\simeq\left|\delta\left(\frac{1}{k_{F}a}\right)\right|^{2}\frac{k_{% B}T_{F}}{|U(\mathbf{0})|}. (122)

In a typical experiment such as the one at Rice we have k_{B}T_{F}<|U(\mathbf{0})|; thus the condition |\mathcal{B}|\ll 1 is satisfied as soon as |\delta(1/k_{F}a)|\ll 1, which is likely to hold at Rice, see (92).

## Appendix C Calculation of dE/d(1/a) in the local density approximation

Restricting to a spin balanced gas, we first consider the general case of the two-channel model (1), in which U_{b}(\mathbf{r})\neq 2U(\mathbf{r}). As explained in the appendix B, where an effective position dependent magnetic field B_{\rm eff}(\mathbf{r}) is introduced, this leads to a position dependent scattering length a(\mathbf{r}). As a consequence, in the local density approximation, here at zero temperature, the gas density n(\mathbf{r}) is given by

 \mu=U(\mathbf{r})+\mu^{\rm hom}[n(\mathbf{r}),1/a(\mathbf{r})] (123)

where \mu^{\rm hom}[n,1/a] is the chemical potential of the free space homogeneous gas of density n and scattering length a. The mean energy and particle numbers in the gas are given in the same approximation by

 \displaystyle E \displaystyle= \displaystyle\int d^{3}r\,\{U(\mathbf{r})+\epsilon^{\rm hom}[n(\mathbf{r}),1/a% (\mathbf{r})]\}\,n(\mathbf{r}) (124) \displaystyle N \displaystyle= \displaystyle\int d^{3}r\,n(\mathbf{r}) (125)

where \epsilon^{\rm hom}[n,1/a] is the mean energy per particle of the free space gas. We wish to calculate the derivative of E with respect to 1/a, where a is the scattering length in the trap center \mathbf{r}=\mathbf{0}, with the constraint that N should be fixed.

Let us perform an infinitesimal variation \delta\left(\frac{1}{a}\right) of 1/a. This leads to a variation \delta\mu of the chemical potential and \delta n(\mathbf{r}) of the density, that may be calculated from (123) but this is not required here. This also leads to a variation of 1/a(\mathbf{r}), according to (95):

 \delta\left(\frac{1}{a(\mathbf{r})}\right)=\left(\frac{1-a_{\rm bg}/a(\mathbf{% r})}{1-a_{\rm bg}/a}\right)^{2}\delta\left(\frac{1}{a}\right). (126)

The variation of the gas energy is then

 \displaystyle\delta E=\int d^{3}r\,\,\delta n(\mathbf{r})\left[U(\mathbf{r})+% \epsilon^{\rm hom}[n(\mathbf{r}),1/a(\mathbf{r})]\phantom{\frac{1}{1}}\right.% \\ \displaystyle\left.+n(\mathbf{r})\frac{\partial\epsilon^{\rm hom}}{\partial n}% [n(\mathbf{r}),1/a(\mathbf{r})]\right]\\ \displaystyle+\delta\left(\frac{1}{a(\mathbf{r})}\right)n(\mathbf{r})\frac{% \partial\epsilon^{\rm hom}}{\partial(1/a)}[n(\mathbf{r}),1/a(\mathbf{r})]. (127)

Writing \epsilon^{\rm hom}[n,1/a]=E^{\rm hom}[N,V,1/a]/N, where V is the volume and N the particle number of the homogeneous gas, we find the relation

 n\,\frac{\partial\epsilon^{\rm hom}}{\partial n}[n,1/a]=\mu^{\rm hom}[n,1/a]-% \epsilon^{\rm hom}[n,1/a]. (128)

Using this expression and (123) one finds that a simplification occurs in (127). Furthermore, the fact that the particle number is fixed imposes

 \int d^{3}r\,\delta n(\mathbf{r})=0 (129)

 \displaystyle\left(\frac{dE}{d(1/a)}\right)_{N}=\int d^{3}r\,n(\mathbf{r})% \left(\frac{1-a_{\rm bg}/a(\mathbf{r})}{1-a_{\rm bg}/a}\right)^{2}\\ \displaystyle\times\frac{\partial\epsilon^{\rm hom}}{\partial(1/a)}[n(\mathbf{% r}),1/a(\mathbf{r})]. (130)

We note that an alternative derivation of this result is possible: One can introduce the density of closed-channel molecules n_{b}^{\rm hom}[n,1/a] for a free space gas of density n and scattering length a, a quantity that can be related to the derivative of \epsilon^{\rm hom}[n,1/a] with respect to 1/a simply by applying the general formula (8) to a free space system. Then setting for the trapped system in the local density approximation N_{b}=\int d^{3}r\,n_{b}^{\rm hom}[n(\mathbf{r}),1/a(\mathbf{r})] and using (8) for the trapped system leads to (130).

In practice, in the present paper, we perform explicit calculations in the limiting cases presented in §2.2: One takes the zero-range limit, which also allows to neglect the spatial variation of a(\mathbf{r}) under conditions discussed in the appendix B. Then \epsilon^{\rm hom} is given by the relation (16) involving some universal function f(x). For the free space chemical potential we then introduce the convenient parametrization

 \mu^{\rm hom}[n,1/a]=\frac{\hbar^{2}}{2ma^{2}}u\left(\frac{1}{k_{F}^{\rm hom}a% }\right) (131)

where k_{F}^{\rm hom}=(3\pi^{2}n)^{1/3} is the ideal gas Fermi wavector, and the function u is

 u(x)=\frac{f(x)-\frac{1}{5}xf^{\prime}(x)}{x^{2}}. (132)

One also restricts to a harmonic trap for U(\mathbf{r}), see (13). The surface n(\mathbf{r})=0 is then an ellipsoid with Thomas-Fermi radii R_{\alpha} along the directions \alpha.

For a<0 we set

 \mu=\frac{\hbar^{2}q^{2}}{2m} (133)

with q>0. Then \omega_{\alpha}R_{\alpha}=\hbar q/m, and the density depends only on X\geq 0 such that X^{2}=\sum_{\alpha}r_{\alpha}^{2}/R_{\alpha}^{2}. Equation (123) then reduces to the implicit equation

 q^{2}a^{2}(1-X^{2})=u\left[\frac{1}{k_{F}^{\rm hom}(X)a}\right] (134)

that one can solve numerically, using the interpolation for f(x) given in appendix D.

For a>0, \mu^{\rm hom}[n,1/a] tends to -\hbar^{2}/(2ma^{2}) for n\to 0 so that we set

 \mu=\frac{\hbar^{2}q^{2}}{2m}-\frac{\hbar^{2}}{2ma^{2}} (135)

with q>0. Then again \omega_{\alpha}R_{\alpha}=\hbar q/m, and (123) reduces to

 q^{2}a^{2}(1-X^{2})=1+u\left[\frac{1}{k_{F}^{\rm hom}(X)a}\right] (136)

to be solved numerically in general.

One can then calculate dE/d(1/a) and k_{F} defined in (14) from the expressions valid on both sides of the resonance:

 \displaystyle\frac{dE}{d(1/a)} \displaystyle= \displaystyle\frac{3N\hbar^{2}}{10m}\frac{\int_{0}^{1}dX\,X^{2}k_{F}^{\rm hom}% (X)^{4}f^{\prime}\left[\frac{1}{k_{F}^{\rm hom}(X)a}\right]}{\int_{0}^{1}dX\,X% ^{2}k_{F}^{\rm hom}(X)^{3}} (137) \displaystyle k_{F} \displaystyle= \displaystyle\frac{\sqrt{2q}}{\pi^{1/3}}\left[4\pi\int_{0}^{1}dX\,X^{2}k_{F}^{% \rm hom}(X)^{3}\right]^{1/6} (138)

The analytical calculations performed in the weakly interacting regimes are also obtained from the above formulas.

## Appendix D Interpolation formula for the homogeneous gas energy

We explain how an interpolation formula was constructed for the function f defined in (16), which was then used to produce the solid line in Fig.1a and Fig.1b.

On the BEC side of the resonance, we applied a cubic spline interpolation to the points in Table I of Giorgini () with values of x=1/(k_{F}^{\rm hom}a) equal to 0, 1, 2, 4 and 6, using the last column of Table I (the one said to give E/N-\epsilon_{b}/2) and adding to the corresponding data the quantity -\hbar^{2}/(2ma^{2}). We added a home made point at x=8 where the value of f(x) and its first order derivative f^{\prime}(x) are obtained from the bosonic Lee-Huang-Yang expansion (2.2).

On the BCS side of the resonance we used the empirical interpolation formula of Salasnich (),

 f_{S}(x)=\alpha_{1}-\alpha_{2}\arctan[\alpha_{3}x(\beta_{1}-x)/(\beta_{2}-x)]. (139)

Constraints on \alpha_{1},\alpha_{2},\alpha_{3} and \beta_{1},\beta_{2} are obtained from the values f(0), f^{\prime}(0) and from the large |x| expansion f(x)=1+c_{1}/x+c_{2}/x^{2}+O(1/x^{3}), where the coefficients c_{1} and c_{2} are given in (22). This differs from Salasnich () where the coefficient c_{2} was not used and \beta_{2} was deduced from a fit to the Monte Carlo data of Giorgini () on the BCS side. Another difference with Salasnich () is that we took the value f^{\prime}(0)=-0.95, that we extracted from Lobo () using the relation (34). Also, we set the value of f(0) to 0.41, which is within the error bars of Giorgini (), in order to reduce the spurious discontinuity of the second order derivative of f(x) in x=0 between the interpolation formula for x<0 and the spline for x>0. All this leads to \alpha_{1}\simeq 0.41,\alpha_{2}\simeq 0.3756,\alpha_{3}\simeq 1.062,\beta_{1}% \simeq 0.9043,\beta_{2}\simeq 0.3797.

A test of this interpolation formulation is to compare its prediction for f^{\prime}(-1) with the value \simeq-0.13 that one can extract from Lobo () using (34); one gets the satisfactory result f_{S}^{\prime}(-1)\simeq-0.14.

## References

• (1) S. Jochim, M. Bartenstein, A. Altmeyer, G. Hendl, S. Riedl, C. Chin, J. Hecker Denschlag, R. Grimm, Science 302, 2101 (2003)
• (2) M. Greiner, C. Regal, D. Jin, Nature, 426, 537 (2003)
• (3) M. Zwierlein, C. Stan, C. Schunck, S. Raupach, S. Gupta, Z. Hadzibabic, W. Ketterle, Phys. Rev. Lett. 91, 250401 (2003)
• (4) T. Bourdel, L. Khaykovich, J. Cubizolles, J. Zhang, F. Chevy, M. Teichmann, L. Tarruell, S.J.J.M.F. Kokkelmans, C. Salomon, Phys. Rev. Lett. 93, 050401 (2004)
• (5) G.B. Partridge, K.E. Strecker, R.I. Kamar, M.W. Jack, R.G. Hulet, Phys. Rev. Lett. 95, 020404 (2005)
• (6) K. O’Hara, S. Hemmer, M. Gehm, S. Granade, J. Thomas, Science 298, 2179 (2002)
• (7) C. Regal, C. Ticknor, J. Bohn, D. Jin, Nature 424, 47 (2003)
• (8) T. Bourdel, J. Cubizolles, L. Khaykovich, K. Magalhaes, S. Kokkelmans, G. Shlyapnikov, C. Salomon, Phys. Rev. Lett. 91, 020402 (2003)
• (9) M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C. Chin, J. Hecker Denschlag, R. Grimm, Phys. Rev. Lett. 92, 120401 (2004)
• (10) C. Chin, M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, J. Hecker Denschlag, R. Grimm, Science 305, 1128 (2004)
• (11) J.E. Thomas, J. Kinast, A. Turlapov, Phys. Rev. Lett. 95, 120402 (2005)
• (12) A. Altmeyer, S. Riedl, C. Kohstall, M.J. Wright, R. Geursen, M. Bartenstein, C. Chin, J. Hecker Denschlag, R. Grimm, Phys. Rev. Lett. 98, 040401 (2007)
• (13) M. Zwierlein, J. Abo-Shaeer, A. Schirotzek, C. Schunck, W. Ketterle, Nature 435, 1047 (2005)
• (14) J.T. Stewart, J.P. Gaebler, C.A. Regal, D.S. Jin, Phys. Rev. Lett. 97, 220406 (2006)
• (15) L. Luo, B. Clancy, J. Joseph, J. Kinast, J.E. Thomas, Phys. Rev. Lett. 98, 080402 (2007)
• (16) G.B. Partridge, W. Li, R.I. Kamar, Y. Liao, R.G. Hulet, Science 311, 503 (2006)
• (17) Y. Shin, M. Zwierlein, C. Schunck, A. Schirotzek, W. Ketterle, Phys. Rev. Lett. 97, 030401 (2006)
• (18) E. Timmermans, P. Tommasini, M. Hussein, A. Kerman, Phys. Rep. 315, 199 (1999)
• (19) M. Holland, S.J.J.M.F. Kokkelmans, M.L. Chiofalo, R. Walser, Phys. Rev. Lett. 87, 120406 (2001)
• (20) J.N. Milstein, S.J.J.M.F. Kokkelmans, M.J. Holland, Phys. Rev. A 66, 043604 (2002)
• (21) T. Köhler, K. Góral, P.S. Julienne, Rev. Mod. Phys. 78, 1311 (2006)
• (22) D. Lee, Phys. Rev. B 73, 11511 (2006)
• (23) A. Bulgac, J.E. Drut, P. Magierski, Phys. Rev. Lett. 96, 090404 (2006)
• (24) E. Burovski, N. Prokof’ev, B. Svistunov, M. Troyer, Phys. Rev. Lett. 96, 160402 (2006)
• (25) E. Burovski, N. Prokof’ev, B. Svistunov, M. Troyer, New J. Phys. 8, 153 (2006)
• (26) O. Juillet, New J. Phys. 9, 163 (2007)
• (27) A. Bulgac, J.E. Drut, P. Magierski, Phys. Rev. A 78, 023625 (2008)
• (28) E. Burovski, E. Kozik, N. Prokof’ev, B. Svistunov, M. Troyer, Phys. Rev. Lett. 101, 090402 (2008)
• (29) Y. Castin, Lecture notes of the 2003 Les Houches Spring School Quantum Gases in Low Dimensions, M. Olshanii, H. Perrin, L. Pricoupenko eds., J. Phys. IV France 116, 89-132 (2004)
• (30) L. Pricoupenko, Y. Castin, J. Phys. A 40, 12863 (2007)
• (31) J. Carlson, S.Y. Chang, V.R. Pandharipande, K.E. Schmidt, Phys. Rev. Lett. 91, 050401 (2003)
• (32) G.E. Astrakharchik, J. Boronat, J. Casulleras, S. Giorgini, Phys. Rev. Lett. 93, 200404 (2004)
• (33) J. Carlson, S. Reddy, Phys. Rev. Lett. 95, 060401 (2005)
• (34) C. Lobo, I. Carusotto, S. Giorgini, A. Recati, S. Stringari, Phys. Rev. Lett. 97, 100405 (2006)
• (35) C. Lobo, A. Recati, S. Giorgini, S. Stringari, Phys. Rev. Lett. 97, 200403 (2006)
• (36) S. Pilati, S. Giorgini, Phys. Rev. Lett. 100, 030401 (2008)
• (37) D.S. Petrov, Phys. Rev. Lett. 93, 143201 (2004)
• (38) Y. Castin, in: Proceedings of the International School of Physics “Enrico Fermi” on Ultra-Cold Fermi gases, edited by M. Inguscio, W. Ketterle and C. Salomon (SIF, Bologna, 2007)
• (39) S. Tan, Ann. Phys. 323, 2952 (2008)
• (40) S. Tan, Ann. Phys. 323, 2971 (2008)
• (41) N.N. Bogoliubov, Lectures on quantum statistics, Vol. 1 (Gordon and Breach, New York, 1967)
• (42) X. Leyronas, R. Combescot, Phys. Rev. Lett. 99, 170402 (2007)
• (43) K. Huang, C.N. Yang, Phys. Rev. 105, 767 (1957)
• (44) T.D. Lee, C.N. Yang, Phys. Rev. 105, 1119 (1957)
• (45) D.S. Petrov, C. Salomon, G.V. Shlyapnikov, Phys. Rev. Lett. 93, 090404 (2004)
• (46) I.V. Brodsky, M.Y. Kagan, A.V. Klaptsov, R. Combescot, X. Leyronas, Phys. Rev. A 73, 032724 (2006)
• (47) A.A. Abrikosov, I.M. Khalatnikov, Sov. Phys. JETP 6, 888 (1958), [J.E.T.P. (USSR) 33, 1154 (1957)]
• (48) V.M. Galitskii, Sov. Phys. JETP 7, 104 (1958), [J.E.T.P. (USSR) 34, 151 (1958)]
• (49) L. Pricoupenko, Y. Castin, Phys. Rev. A 69, 051601 (2004)
• (50) A. Bulgac, G.F. Bertsch, Phys. Rev. Lett. 94, 070401 (2005)
• (51) R. Combescot, X. Leyronas, Phys. Rev. Lett. 93, 138901 (2004)
• (52) J. Kinast, A. Turlapov, J.E. Thomas, Phys. Rev. A 70, 051401 (2004)
• (53) M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C. Chin, J. Hecker Denschlag, R. Grimm, Phys. Rev. Lett. 92, 203201 (2004), see also: A. Altmeyer, S. Riedl, C. Kohstall, M.J. Wright, J. Hecker Denschlag, R. Grimm, arXiv:cond-mat/0611285v2
• (54) G.E. Astrakharchik, R. Combescot, X. Leyronas, S. Stringari, Phys. Rev. Lett. 95, 030404 (2005)
• (55) M. Urban, Phys. Rev. A 75, 053607 (2007)
• (56) M. Bartenstein, A. Altmeyer, S. Riedl, R. Geursen, S. Jochim, C. Chin, J. Hecker Denschlag, R. Grimm, A. Simoni, E. Tiesinga et al., Phys. Rev. Lett. 94, 103201 (2005)
• (57) G.M. Falco, H.T.C. Stoof, Phys. Rev. A 71, 063614 (2005)
• (58) R. Hulet, private communication
• (59) E. Braaten, L. Platter, Phys. Rev. Lett. 100, 205301 (2008)
• (60) D.M. Gangardt, G.V. Shlyapnikov, Phys. Rev. Lett. 90, 010401 (2003)
• (61) T. Kinoshita, T. Wenger, D.S. Weiss, Phys. Rev. Lett. 95, 190406 (2005)
• (62) M. Yasuda, F. Shimizu, Phys. Rev. Lett. 77, 3090 (1996)
• (63) M. Schellekens, R. Hoppeler, A. Perrin, J.V. Gomes, D. Boiron, A. Aspect, C.I. Westbrook, Science 310, 648 (2005)
• (64) T. Jeltes, J.M. McNamara, W. Hogervorst, W. Vassen, V. Krachmalnicoff, M. Schellekens, A. Perrin, H. Chang, D. Boiron, A. Aspect et al., Nature 445, 402 (2007)
• (65) C. Lobo, S.D. Gensemer, Phys. Rev. A 78, 023618 (2008)
• (66) E.H. Lieb, W. Liniger, Phys. Rev. 130, 1605 (1963)
• (67) J. Yvon, C. R. Acad. Sci. 246, 2858 (1958)
• (68) R. Snider, J. Chem. Phys. 32, 1051 (1960)
• (69) K. Huang, Statistical Mechanics (Wiley, New York, 1963)
• (70) M. Holzmann, Y. Castin, Eur. Phys. J. D 7, 425 (1999)
• (71) A.J. Leggett, Lecture 2 at the international quantum gases workshop (Institut Henri Poincaré, Paris, 2007), www.phys.ens.fr/~castin/files_pdf.html
• (72) C.A. Regal, M. Greiner, S. Giorgini, M. Holland, D.S. Jin, Phys. Rev. Lett. 95, 250404 (2005)
• (73) L. Tarruell, M. Teichmann, J. McKeever, T. Bourdel, J. Cubizolles, L. Khaykovich, J. Zhang, N. Navon, F. Chevy, C. Salomon, in: Proceedings of the International School of Physics “Enrico Fermi” on Ultra-Cold Fermi gases, edited by M. Inguscio, W. Ketterle and C. Salomon (SIF, Bologna, 2007)
• (74) Q. Chen, C.A. Regal, D.S. Jin, K. Levin, Phys. Rev. A 74, 011601 (2006)
• (75) M.L. Chiofalo, S. Giorgini, M. Holland, Phys. Rev. Lett. 97, 070404 (2006)
• (76) G.E. Astrakharchik, J. Boronat, J. Casulleras, S. Giorgini, Phys. Rev. Lett. 95, 230405 (2005)
• (77) M. Olshanii, V. Dunjko, Phys. Rev. Lett. 91, 090401 (2003)
• (78) Y. Castin, I. Carusotto, Optics Comm. 243, 81 (2004)
• (79) M. Punk, W. Zwerger, Phys. Rev. Lett. 99, 170404 (2007)
• (80) G. Baym, C.J. Pethick, Z. Yu, M.W. Zwierlein, Phys. Rev. Lett. 99, 190407 (2007)
• (81) F. Werner, Y. Castin, Phys. Rev. A 74, 053604 (2006)
• (82) J.E. Thomas, Phys. Rev. A 78, 013630 (2008)
• (83) S. Tan, Ann. Phys. 323, 2987 (2008)
• (84) F. Werner, Phys. Rev. A 78, 025601 (2008)
• (85) E. Braaten, D. Kang, L. Platter, arXiv:0806.2277v1
• (86) S. Zhang, A.J. Leggett, arXiv:0809.1892v1
• (87) S. Simonucci, P. Pieri, G. Strinati, Europhys. Lett. 69, 713 (2005)
• (88) I. Schmidt-Mink, W. Müller, W. Meyer, Chem. Phys. 92, 263 (1985)
• (89) R. Hulet, private communication
• (90) N. Manini, L. Salasnich, Phys. Rev. A 71, 033625 (2005)
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters