text:
# Superfluid transition temperature in a trapped gas of Fermi atoms with a Feshbach resonance

## I Introduction

The search for the superfluid phase transition in a trapped gas of Fermi atoms is one of the most exciting problems in current physics. This involves the formation of BCS Cooper-pairs made up of degenerate Fermions of two different species (“spin up” and “spin down”). Degenerate Fermi gases have been realized in 40K and 6Li. Very recently, a Feshbach resonance was observed in 40K working with the two hyperfine states $`|9/2,9/2`$ and $`|9/2,7/2`$. In the Feshbach resonance, a quasi-molecular Boson is formed, which can produce an attractive pairing interaction between two Fermi atoms of opposite spins. As a result, one expects a BCS-type superfluidity with a high phase transition temperature $`T_\mathrm{c}`$ due to this strong attractive interaction. However, such strong interaction is also known to enhance fluctuations, which can strongly suppress $`T_\mathrm{c}`$ predicted by the usual mean-field BCS theory. Thus the inclusion of these Cooper-pair channel fluctuations is crucial in considering the high-$`T_\mathrm{c}`$ superfluidity induced by an atomic Feshbach resonance.

In a previous paper, we have discussed how this fluctuation contribution affects the superfluid phase transition temperature in a uniform Fermi gas, extending the theory developed by Nozières and Schmitt-Rink: We assume that $`N_{}=N_{}`$, where $`N_\sigma `$ is the number of Fermi atoms with spin-$`\sigma `$, and that both Fermion components have the same single-particle energy. As the threshold energy of the Feshbach resonance is lowered, the character of the phase transition is found to continuously change from the BCS-type to BEC-type. In the latter regime, $`T_\mathrm{c}`$ is strongly suppressed compared with the usual mean-field BCS theory and it approaches $`T_\mathrm{c}=0.218T_\mathrm{F}`$ in the BEC limit. (Here $`T_\mathrm{F}`$ is the Fermi temperature of the non-interacting spin-up Fermi gas.) We show that the stable Bosons which appear as we enter the BEC regime are a strongly hybridized mixture of Feshbach molecules and preformed Cooper-pairs. The Feshbach component of these Bosons dominates the excitation spectrum in the BEC limit.

In this paper, we extend our previous work by including the effect of a harmonic trap potential. This extension is clearly necessary since all experiments on Fermi gases are done in some sort of trap potential. So far, the trapping potential has been examined within the mean-field BCS theory using the simple local density approximation (LDA). In this paper, we also use the LDA but go past the mean-field approximation to include the effect of particle-particle fluctuations. We show how the BCS-BEC crossover can be observed from characteristic changes in the static atomic density profile, easily measured by standard techniques.

Besides the phase transition temperature $`T_\mathrm{c}`$, it is also important to understand how strong-coupling fluctuation effects enter above $`T_\mathrm{c}`$. In particular, it is an interesting problem to clarify at what temperature the stable Feshbach molecules and preformed Cooper-pairs in the BEC regime first appear in the normal-fluid phase above $`T_\mathrm{c}`$. Our calculations show that above $`T_\mathrm{c}`$, there exists a characteristic temperature $`T^{}`$ which describes the crossover from a free Fermi gas to a molecular Bose gas.

In this paper, our treatment of the particle-particle fluctuations assumes that their coupling to the Feshbach resonance is weak and can be thus treated perturbatively. Most of our discussion is for this case, in which the BCS-BEC crossover is driven by decreasing the value of the Feshbach resonance threshold ($`2\nu `$) relative to twice the bare Fermi energy ($`2\epsilon _\mathrm{F}`$). When the Feshbach coupling strength is very large, as Kokkelmans and coworkers have found, it is not clear that our simple treatment of the particle-particle fluctuations is valid. However, in this limit, one finds that the BCS-BEC crossover occurs for values $`2\nu 2\epsilon _\mathrm{F}`$ because of the broad Feshbach resonance (see also Ref.). Effectively, one is in the strong-coupling regime discussed by Nozières and Schmitt-Rink, where the BEC phase is associated with preformed Cooper-pairs. When the threshold energy $`2\nu `$ becomes comparable to $`2\epsilon _\mathrm{F}`$, our calculations again show that the Feshbach resonance takes over as the dominant component of the BEC phase.

This paper is organized as follows: In Section II, we review the coupled Fermion-Boson model used in earlier work and discuss the LDA extension of the results we presented in Ref.. In Section III, we use this opportunity to expand on certain aspects which were not discussed in detail in Ref.  due to lack of space. The BCS-BEC crossover in a trapped system is discussed in Section IV. In Section V, we consider the character of the composite Bosons. The crossover from a free Fermi gas to a gas of long-lived molecular Bosons above $`T_\mathrm{c}`$ is discussed in Section VI. Section VII considers the case of a strong Feshbach coupling.

## II Formulation

### A Coupled Fermion-Boson gas in a trap

The gas of interacting Fermi atoms coupled to a Feshbach resonance can be described by the coupled Fermion-Boson model given by

$``$ $`=`$ $`{\displaystyle \underset{\text{p}\sigma }{}}\epsilon _\text{p}c_{\text{p}\sigma }^{}c_{\text{p}\sigma }+{\displaystyle \underset{\text{q}}{}}(E_\text{q}^0+2\nu )b_\text{q}^{}b_\text{q}`$ (1)
$``$ $`U{\displaystyle \underset{\text{p},\text{p}^{}}{}}c_\text{p}^{}c_\text{p}^{}c_\text{p}^{}c_\text{p}^{}+g_\mathrm{r}{\displaystyle \underset{\text{p},\text{q}}{}}[b_\text{q}^{}c_{\text{p}+\text{q}/2}c_{\text{p}+\text{q}/2}+\mathrm{h}.\mathrm{c}].`$ (2)
Here $`c_{\text{p}\sigma }^{}`$ is the creation operator of a Fermi atom with the kinetic energy $`\epsilon _\text{p}=p^2/2m`$. We assume that the system involves two atomic hyper-fine states which will lead to Cooper-pair Bosons and hence superfluidity. We describe these atomic states for simplicity using the pseudo-spin language $`\sigma =,`$. The quasi-molecular Boson of momentum $`𝐪`$ associated with the Feshbach resonance is described by the operator $`b_\text{q}^{}`$. $`E_\text{q}^0q^2/2M`$ is the kinetic energy of this Boson and $`2\nu `$ represents the bottom of this resonance Bose spectrum, to be referred to as the threshold energy of the Feshbach resonance. The coupling to this Feshbach resonance Boson is described by last term in eq. (2) with a strength $`g_\mathrm{r}`$, in which two Fermi atoms form one $`b`$-Boson and the $`b`$-Boson breaks into two free Fermi atoms. The Hamiltonian also includes a standard BCS pairing interaction with the coupling constant $`U<0`$, originating from non-resonant interactions between two Fermi atoms ($`\sigma =,`$). Even assuming that this direct attractive pairing interaction $`U`$ is weak, the effective pairing interaction strength can be strong due to the effect of the coupling to the Feshbach resonance. In this paper, we do not discuss the question of how good the model Hamiltonian in eq. (2) is in the context of ultracold Fermi gases, but simply refer to the literature .

Since in this model, one Feshbach-induced $`b`$-molecule consists of two Fermi atoms, we have $`M=2m`$ and also impose the conservation of the total number of particles as

$`N`$ $`=`$ $`{\displaystyle \underset{\text{p}\sigma }{}}c_{\text{p}\sigma }^{}c_{\text{p}\sigma }+2{\displaystyle \underset{\text{q}}{}}b_\text{q}^{}b_\text{q}`$ (3)
$``$ $`N_\mathrm{F}+2N_\mathrm{B}.`$ (4)
As discussed in our previous paper, this constraint is incorporated into the model Hamiltonian in eq. (2) by replacing $`\mu N`$, where $`\mu `$ is the chemical potential. The resulting Hamiltonian has the same form as eq. (2), but the kinetic energies of the Fermions and Bosons are now given by $`\epsilon _\text{p}\epsilon _\text{p}\mu `$ and $`E_\text{q}^0E_\text{q}^02\mu `$, respectively.

Now we introduce the effect of a harmonic trap. In this paper, we consider an optical trap in which the two hyperfine Fermi atom states feel the same trap potential given by

$$V_{\mathrm{trap}}(𝐫)=\frac{1}{2}\underset{j}{}m[\omega _{0x}^2x_j^2+\omega _{0y}^2y_j^2+\omega _{0z}^2z_j^2],$$
(5)
where $`𝐫_j=(x_j,y_j,z_j)`$ is the position of the $`j`$-th atom. In our model, we also assume that the Feshbach molecules also feel this potential with $`m`$ being replaced by the Boson mass $`M=2m`$. We could relax this last assumption to deal more realistically with the dipole moment of the composite molecule but we leave this extension to future studies. Treating the resulting trapped Fermion-Boson gas within the LDA amounts to replacing the chemical potential $`\mu `$ by

$$\mu (𝐫)=\mu V_{\mathrm{trap}}(𝐫).$$
(6)
The constraint on the total number of Fermi atoms is then replaced by $`N=𝑑𝐫[N_\mathrm{F}(𝐫)+2N_\mathrm{B}(𝐫)]`$, where $`N_\mathrm{F}(𝐫)`$ and $`N_\mathrm{B}(𝐫)`$ are the number density of Fermi atoms and Bose molecules, respectively.

The atomic hyperfine states may feel different trap potentials in a magnetic trap. This situation is similar to electron spins in superconducting metals under a magnetic field. In the theory of superconductivity, it is known that a magnetic field suppresses superconductivity and, under certain conditions, produces a non-uniform kind of superconducting state called the Fulde-Ferrell state. The effect of fluctuations on this kind of Zeeman splitting in atomic Fermi gases is considered in another paper.

### B Equations for $`T_\mathrm{c}`$: Thouless criterion

The superfluid phase transition temperature $`T_\mathrm{c}`$ can be conveniently determined using the Thouless criterion, in which the phase transition is characterized by a singularity in the particle-particle scattering vertex function $`\mathrm{\Gamma }`$ at $`\omega =𝐪=0`$, associated with the formation of Cooper-pairs. The four-point vertex function $`\mathrm{\Gamma }`$ within the generalized $`t`$-matrix approximation in terms of $`U`$ and $`g_\mathrm{r}`$ is given by the diagrams in Fig.1. In this figure, the first line describes ladder processes associated with the attractive interaction $`U<0`$, processes familiar in strong coupling superconductivity in electronic systems. On the other hand, the second line includes the effective Fermion-Fermion interaction mediated by the $`b`$-Boson propagator $`D_0`$ defined in eq. (11), scattering processes which only arise in atomic Fermi gases coupled to a Feshbach resonance, such as in eq. (2). The equation of the vertex function $`\mathrm{\Gamma }`$ is given by

$`\mathrm{\Gamma }`$ $`(p_+,p_{};p_+^{},p_{}^{};𝐫)=(Ug_\mathrm{r}^2D_0(q;𝐫))`$ (7)
$`+`$ $`(Ug_\mathrm{r}^2D_0(q;𝐫)){\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\text{p}^{\prime \prime },\omega _m^{\prime \prime }}{}}G_0(p_+^{\prime \prime };𝐫)G(p_{}^{\prime \prime };𝐫)\mathrm{\Gamma }(p_+^{\prime \prime },p_{}^{\prime \prime };p_+^{},p_{}^{};𝐫),`$ (8)
where $`p_+(𝐩+𝐪/2,i\omega _m+i\nu _n)`$, $`p_{}(𝐩𝐪/2,i\omega _m)`$, $`q(𝐪,i\nu _m)`$ and $`\beta 1/T`$ is the inverse of temperature; $`i\omega _m`$ and $`i\nu _n`$ are the Fermion and Boson Matsubara frequencies, respectively. $`G_0`$ and $`D_0`$ represent the bare one-particle Fermion and Boson thermal Green functions given by

$`\{\begin{array}{c}G_0(𝐩,i\omega _m;𝐫)={\displaystyle \frac{1}{i\omega _m\epsilon _\text{p}+\mu (𝐫)}},\hfill \\ D_0(𝐪,i\nu _n;𝐫)={\displaystyle \frac{1}{i\nu _nE_\text{q}^02\nu +2\mu (𝐫)}}.\hfill \end{array}`$ (11)
Since the modified chemical potential $`\mu (𝐫)`$ includes the trap potential, the LDA Green functions and vertex function $`\mathrm{\Gamma }`$ also depend on position $`𝐫`$. Multiplying both the sides of eq. (8) by $`\frac{1}{\beta }G_0(p_+;𝐫)G_0(p_{};𝐫)`$ and summing over $`\omega _m`$ and $`𝐩`$, we obtain

$$\frac{1}{\beta }\underset{\text{p}^{\prime \prime },\omega _m^{\prime \prime }}{}G_0(p_+^{\prime \prime };𝐫)G(p_{}^{\prime \prime };𝐫)\mathrm{\Gamma }(p_+^{\prime \prime },p_{}^{\prime \prime };p_+^{},p_{}^{};𝐫)=\frac{(Ug_\mathrm{r}^2D_0(q;𝐫))\mathrm{\Pi }(q;𝐫)}{1(Ug_\mathrm{r}^2D_0(q;𝐫))\mathrm{\Pi }(q;𝐫)}.$$
(12)
Here $`\mathrm{\Pi }(q;𝐫)=\mathrm{\Pi }(𝐪,i\nu _n;𝐫)`$ is the well-known particle-particle correlation function of the Cooper-pair-field operator defined by $`\widehat{\mathrm{\Delta }}(𝐪,\tau ;𝐫)_\text{p}c_{\text{p}+\text{q}/2}(\tau ;𝐫)c_{\text{p}+\text{q}/2}(\tau ;𝐫)`$ in the absence of $`U`$ and $`g_\mathrm{r}`$:

$`\mathrm{\Pi }(𝐪,i\nu _n;𝐫)`$ $`=`$ $`{\displaystyle _0^\beta }𝑑\tau e^{i\nu _n\tau }T_\tau \{\widehat{\mathrm{\Delta }}(𝐪,\tau ;𝐫)\widehat{\mathrm{\Delta }}^{}(𝐪,0;𝐫)\}`$ (13)
$`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\text{p},\omega _m}{}}G_0(𝐩+𝐪/2,i\omega _m+i\nu _n;𝐫)G_0(𝐩+𝐪/2,i\omega _m;𝐫)`$ (14)
$`=`$ $`{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{1f(\epsilon _{\text{p}+\text{q}/2}\mu (𝐫))f(\epsilon _{\text{p}\text{q}/2}\mu (𝐫))}{\epsilon _{\text{p}+\text{q}/2}+\epsilon _{\text{p}\text{q}/2}2\mu (𝐫)i\nu _n}},`$ (15)
where $`f(\epsilon )`$ is the Fermi distribution function. Substituting eq. (12) into eq. (8), we obtain

$$\mathrm{\Gamma }(𝐪,i\nu _n;𝐫)=\frac{Ug_\mathrm{r}^2D_0(𝐪,i\nu _n;𝐫)}{1[Ug_\mathrm{r}^2D_0(𝐪,i\nu _n;𝐫)]\mathrm{\Pi }(𝐪,i\nu _n;𝐫)}.$$
(16)
We shall see that the effective attractive pairing interaction is given by

$$U_{\mathrm{eff}}(𝐪,i\nu _n;𝐫)=Ug_\mathrm{r}^2D_0(𝐪,i\nu _n;𝐫).$$
(17)
According to the Thouless criterion, $`T_\mathrm{c}`$ is given as the temperature at which the vertex function $`\mathrm{\Gamma }(𝐪,i\nu _n;𝐫)`$ in eq. (16) develops a pole at $`𝐪=\nu _n=0`$. Within the LDA, since $`\mathrm{\Gamma }`$ depends on position $`𝐫`$, this condition is satisfied at different temperatures $`T_0(𝐫)`$ depending on $`𝐫`$. In this approximation, we have to regard the highest $`T_0(𝐫)`$ as the phase transition temperature of the trapped gas. Since the density of particles is maximum at the center of the trap, we have $`T_\mathrm{c}=T_0(𝐫=0)`$. Then the equation for $`T_\mathrm{c}`$ within the LDA is given by

$$1=U_{\mathrm{eff}}\mathrm{\Pi }(0,0;𝐫=0),$$
(18)
where

$$U_{\mathrm{eff}}U_{\mathrm{eff}}(𝐪=0,i\nu _n=0;𝐫=0)=U+\frac{g_\mathrm{r}^2}{2\nu 2\mu }$$
(19)
is the static part of the effective pairing interaction in the long wave length limit at $`𝐫=0`$. Equation (18) indicates that even if the direct pairing interaction $`U`$ is weak, the effective pairing interaction $`U_{\mathrm{eff}}`$ in eq. (19) can be large if $`g_\mathrm{r}`$ is large or when the chemical potential of the $`b`$-Boson ($`\mu _B2\mu `$) approaches the bottom of the the Bose band $`2\nu `$. More explicitly, we note that eq. (18) can be written in the well known BCS form

$$1=(U+g^2\frac{1}{2\nu 2\mu })\underset{\text{p}}{}\frac{\mathrm{tanh}(\epsilon _\text{p}\mu )/2T_\mathrm{c}}{2\epsilon _\text{p}2\mu }.$$
(20)
As shown in the literature, the pole which develops at $`T_\mathrm{c}`$ at $`𝐪=\nu _n=0`$ signals an instability of the normal Fermi gas. At $`T<T_\mathrm{c}`$, the amplitude of this two-particle bound state can be shown to grow exponentially with time. This instability is removed by the appearance of a new phase below $`T_\mathrm{c}`$, described by the Cooper-pair order parameter.

### C Fluctuation contributions to the chemical potential

The chemical potential $`\mu `$ in eq. (20) is determined by the equation of the total number of Fermi atoms where a $`b`$-Boson also counts as two Fermi atoms. In the weak-coupling BCS theory in the absence of the $`b`$-Boson, this equation simply gives $`\mu =\epsilon _\mathrm{F}`$ (where $`\epsilon _\mathrm{F}`$ is the Fermi energy in the absence of any $`b`$-Bosons), assuming that the temperature dependence of $`\mu `$ can be neglected. Indeed, in most metallic (weak-coupling) superconductors, the deviation of $`\mu (T_\mathrm{c})`$ from $`\epsilon _\mathrm{F}`$ is negligibly small as $`\delta \mu (T_\mathrm{c}/\epsilon _\mathrm{F})^21`$. On the other hand, when $`g_\mathrm{r}`$ is large or the Feshbach resonance strongly enhances the effective pairing interaction $`U_{\mathrm{eff}}`$, in analogy with the “strong-coupling” superconductivity discussed by Nozières and Schmitt-Rink, the chemical potential $`\mu `$ is expected to deviate remarkably from $`\epsilon _\mathrm{F}`$ due to the effect of fluctuations associated with the Cooper-pair channel. We now discuss these fluctuations.

The equation giving the number of Fermi atoms determines the chemical potential $`\mu `$. This is obtained, as usual, from first calculating the thermodynamic potential $`\mathrm{\Omega }𝑑𝐫\mathrm{\Omega }(𝐫)`$ and then using the identity $`N=\mathrm{\Omega }/\mu `$. The fluctuation contribution to the local thermodynamic potential $`\mathrm{\Omega }(𝐫)`$ is shown in terms of diagrams in Fig. 2. In Fig. 2(a), we show the usual diagrams associated with the direct pairing interaction $`U`$. Working to all orders in $`U`$, this gives

$$\delta \mathrm{\Omega }_U(𝐫)=\frac{1}{\beta }\underset{\text{q},\nu _n}{}e^{i\nu _n\delta }\mathrm{ln}[1U\mathrm{\Pi }(𝐪,i\nu _n;𝐫)].$$
(21)
This is the contribution discussed in Ref.. In Fig. 2(b), we show the new fluctuation diagrams contributing to $`\mathrm{\Omega }(𝐫)`$ associated with the Feshbach coupling interaction $`g_\mathrm{r}`$. Summing up to all orders, this gives

$$\delta \mathrm{\Omega }_{g_\mathrm{r}}(𝐫)=\frac{1}{\beta }\underset{\text{q},\nu _n}{}e^{i\nu _n\delta }\mathrm{ln}[1+g_\mathrm{r}^2D_0(𝐪,i\nu _n;𝐫)\stackrel{~}{\mathrm{\Pi }}(𝐪,i\nu _n;𝐫)],$$
(22)
where

$$\stackrel{~}{\mathrm{\Pi }}(𝐪,i\nu _n;𝐫)\frac{\mathrm{\Pi }(𝐪,i\nu _n;𝐫)}{1U\mathrm{\Pi }(𝐪,i\nu _n;𝐫)}$$
(23)
describes the particle-particle bubble including the effect of $`U`$ in the ladder approximation.

We note that the total fluctuations contribution to $`\mathrm{\Omega }(𝐫)`$ is the sum of the two contributions in eqs. (21) and (22), which can be combined to give

$`\delta \mathrm{\Omega }(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}e^{i\nu _n\delta }\mathrm{ln}\left[[1U\mathrm{\Pi }(𝐪,i\nu _n;𝐫)]\right[1+g_\mathrm{r}^2D_0(𝐪,i\nu _n;𝐫)\stackrel{~}{\mathrm{\Pi }}(𝐪,i\nu _n;𝐫)]]`$ (24)
$`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}e^{i\nu _n\delta }\mathrm{ln}[1U_{\mathrm{eff}}(𝐪,i\nu _n;𝐫)\mathrm{\Pi }(𝐪,i\nu _n;𝐫)],`$ (25)
where the effective pairing interaction $`U_{\mathrm{eff}}(𝐪,i\nu _n;𝐫)`$ is as defined in eq. (17). In other words, the only effect of the diagrams in Fig. 2(b) is to renormalize the pairing interaction as $`UU_{\mathrm{eff}}(𝐪,i\nu _n;𝐫)`$ in eq. (21). We also note that the collective fluctuations in eq. (25) (which are Bosons) are precisely those which are associated with the poles of the particle-particle vertex function $`\mathrm{\Gamma }`$ in eq. (16). The thermodynamic potential for the contribution of the free Fermions and $`b`$-Bosons is given by the sum of

$`\{\begin{array}{c}\mathrm{\Omega }_F^0(𝐫)={\displaystyle \frac{2}{\beta }}{\displaystyle \underset{\text{p},\omega _m}{}}e^{i\omega _m\delta }\mathrm{ln}G_0^1(𝐩,i\omega _m;𝐫),\hfill \\ \mathrm{\Omega }_B^0(𝐫)={\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\text{p},\nu _n}{}}e^{i\nu _n\delta }\mathrm{ln}D_0^1(𝐩,i\nu _n;𝐫).\hfill \end{array}`$ (28)
Combining all these results and using $`N=\mathrm{\Omega }/\mu `$, we find the result

$$N=𝑑𝐫\left[N_\mathrm{F}^0(𝐫)+2N_\mathrm{B}^0(𝐫)\frac{1}{\beta }\underset{\text{q},i\nu _n}{}e^{i\delta \nu _n}\frac{}{\mu }\mathrm{ln}\left[1\mathrm{U}_{\mathrm{eff}}(𝐪,\mathrm{i}\nu _\mathrm{n};𝐫)\mathrm{\Pi }(𝐪,\mathrm{i}\nu _\mathrm{n};𝐫)\right]\right],$$
(29)
where

$`\{\begin{array}{c}N_F^0={\displaystyle \frac{2}{\beta }}{\displaystyle \underset{\text{p},\omega _m}{}}e^{i\omega _m\delta }G_0(𝐪,i\omega _m;𝐫)=2{\displaystyle \underset{\text{p}}{}}f(\epsilon _\text{p}\mu (𝐫)),\hfill \\ N_B^0={\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\text{p},\nu _n}{}}e^{i\nu _n\delta }D_0(𝐪,i\nu _n;𝐫)={\displaystyle \underset{\text{p}}{}}n_B(E_\text{q}^0+2\nu 2\mu (𝐫)).\hfill \end{array}`$ (32)
Here $`n_B(E)`$ represents the Bose distribution function. Equations (18) and (29) provide us with a closed set of two coupled equations, from which we can obtain $`\mu (T,N)`$ and $`T_\mathrm{c}(N)`$. In solving the coupled equations numerically, we have to take into account carefully the convergence factor $`e^{i\nu _n\delta }`$ in the last term in eq. (29). In the appendix, we explain how to sum up the Matsubara frequencies in this term. Above $`T_\mathrm{c}`$, we need only solve eq. (29) for $`\mu (T,N)`$. Our numerical results will be discussed in Sections IV and V. In the next section, however, we first discuss the physics involved in treating the fluctuations within the approximation developed in this section. This will allow us to better understand the nature of the Bosons involved in the crossover to a BEC regime.

## III Physics of Coupled $`b`$-Bosons and Cooper Pairs

The contribution arising from the interaction in eq. (29) can be described in terms of Bose fluctuations given by the solutions of

$$1=U_{\mathrm{eff}}(𝐪,\omega +i\delta ;𝐫)\mathrm{\Pi }(𝐪,\omega +i\delta ;𝐫),$$
(33)
where we made the usual analytic continuation from the imaginary Bose Matsubara frequency to the real frequency axis as $`i\nu _n\omega +i\delta `$. While eq. (33) appears relatively simple, it describes a complex situation involving the formation of preformed Cooper-pairs and their coupling to the bare $`b`$-Bosons in eq. (2) due to the Feshbach coupling interaction $`g_\mathrm{r}`$. The preformed Cooper-pairs and the $`b`$-molecules are coupled to each other, which leads to hybridization and only one kind of composite Bose molecule. It is useful to try to disentangle these effect, although there seems no unique way of doing this in such a coupled excitation problem.

Denoting the number density of atoms given by the last term in eq. (29) as $`\delta N_{\mathrm{FL}}(𝐫)`$, we have

$$\delta N_{\mathrm{FL}}(𝐫)=\frac{1}{\beta }\underset{\text{q},\nu _n}{}e^{i\nu _n\delta }\left[2g_\mathrm{r}^2\mathrm{\Pi }D_0^2+U_{\mathrm{eff}}(q)\frac{\mathrm{\Pi }}{\mu }\right]\left(\frac{1}{1U_{\mathrm{eff}}(q)\mathrm{\Pi }}\right),$$
(34)
where for simplicity we denote $`\mathrm{\Pi }(𝐪,i\nu _n;𝐫)\mathrm{\Pi }`$, $`D_0(𝐪,i\nu _n;𝐫)D_0`$, $`U_{\mathrm{eff}}(𝐪,i\nu _n;𝐫)U_{\mathrm{eff}}(q)`$, and we have used the identity $`D_0/\mu =2D_0^2`$. The effect of the particle-particle fluctuations described by $`\mathrm{\Pi }(𝐪,i\nu _n;𝐫)`$ naturally gives rise to a renormalized $`b`$-Boson propagator $`\stackrel{~}{D}`$ with a self-energy given by

$$\mathrm{\Sigma }(𝐪,i\nu _n;𝐫)=g_\mathrm{r}^2\stackrel{~}{\mathrm{\Pi }}(𝐪,i\nu _n;𝐫),$$
(35)
where $`\stackrel{~}{\mathrm{\Pi }}`$ has been defined in eq. (23). Thus we have

$`\stackrel{~}{D}^1(𝐪,i\nu _n;𝐫)`$ $``$ $`D_0^1+g_\mathrm{r}^2\stackrel{~}{\mathrm{\Pi }}`$ (36)
$`=`$ $`i\nu _n[E_\text{q}^0+2\nu g_\mathrm{r}^2\stackrel{~}{\mathrm{\Pi }}(𝐪,i\nu _n;𝐫)]+2\mu (𝐫).`$ (37)
The self-energy $`\mathrm{\Sigma }`$ in eq. (35) describes the effect of the decay of a $`b`$-Boson into two Fermi atoms and their recombination to a $`b`$-Boson, treating the non-resonant interaction $`U`$ within the ladder approximation. This self-energy can be also obtained by summing up the diagrams shown in the third line in Fig. 1 (shaded bubble in the figure).

One may verify that the first term in the square brackets in eq. (34) renormalizes the number of $`b`$-Bosons, giving

$`2N_\mathrm{B}^0(𝐫)+{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}e^{i\nu _n\delta }{\displaystyle \frac{2g_\mathrm{r}^2\mathrm{\Pi }D_0^2}{1U_{\mathrm{eff}}(q)\mathrm{\Pi }}}`$ $`=`$ $`{\displaystyle \frac{2}{\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}e^{i\nu _n\delta }{\displaystyle \frac{(1U\mathrm{\Pi })D_0}{1U_{\mathrm{eff}}(q)\mathrm{\Pi }}}`$ (38)
$`=`$ $`{\displaystyle \frac{2}{\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}e^{i\nu _n\delta }\stackrel{~}{D}(𝐪,i\nu _n;𝐫)`$ (39)
$``$ $`2\stackrel{~}{N}_\mathrm{B}(𝐫).`$ (40)
We recall that in our interacting Boson-Fermion model Hamiltonian, we have \[see eq. (4)\]

$$N=N_\mathrm{F}+2N_\mathrm{B},$$
(41)
where, within the LDA,

$`\{\begin{array}{c}N_\mathrm{F}={\displaystyle 𝑑𝐫\underset{\text{p},\sigma }{}c_{\text{p}\sigma }^{}(𝐫)c_{\text{p}\sigma }(𝐫)},\hfill \\ N_\mathrm{B}={\displaystyle 𝑑𝐫\underset{\text{q}}{}b_\text{q}^{}(𝐫)b_\text{q}(𝐫)}\hfill \end{array}`$ (44)
give the number of Fermions and Bosons, respectively. As discussed in many body texts, $`N_\mathrm{B}`$ is defined in terms of the renormalized single particle $`b`$-Boson Green function $`\stackrel{~}{D}(𝐪,i\nu _n;𝐫)`$. In fact, one sees that $`N_\mathrm{B}=\stackrel{~}{N}_\mathrm{B}𝑑𝐫\stackrel{~}{N}_\mathrm{B}(𝐫)`$, where $`\stackrel{~}{N}_\mathrm{B}(𝐫)`$ is defined in eq. (40). This means that $`N_\mathrm{F}`$ in eq. (44) is given by

$$N_\mathrm{F}=N_\mathrm{F}^0+2N_\mathrm{C},$$
(45)
where $`2N_\mathrm{C}`$ is obtained from the second term in the square brackets in eq. (34) as

$$2N_\mathrm{C}=\frac{1}{\beta }\underset{\text{q},\nu _n}{}𝑑𝐫e^{i\nu _n\delta }\frac{\mathrm{\Pi }}{\mu }\left(\frac{U_{\mathrm{eff}}(q)}{1U_{\mathrm{eff}}(q)\mathrm{\Pi }}\right).$$
(46)
Equation (45) shows that the expression for $`2N_\mathrm{C}`$ in eq. (46) is the natural definition of the number of Cooper-pair Bosons which arise in the interacting Fermion gas with an attractive interaction given by $`U_{\mathrm{eff}}(q)`$.

The preceding analysis leads to the following expression for $`N`$ given in Ref.:

$$N=N_\mathrm{F}^0+2N_\mathrm{B}+2N_\mathrm{C},$$
(47)
where $`N_\mathrm{B}`$ is the contribution from the $`b`$-Bosons and $`N_\mathrm{C}`$ is the contribution from the Cooper-pairs, both being strongly renormalized through the effect of the Feshbach resonance coupling $`g_\mathrm{r}`$. The problem with the decomposition given in eq. (47) is that when $`g_\mathrm{r}0`$, one must remember that the $`b`$-Bosons and the preformed Cooper-pairs are strongly hybridized. As a result, they give rise to a single Boson branch described by the pole of $`\stackrel{~}{D}`$ in eq. (37), namely

$$\omega _\text{q}=E_\text{q}^0+2\nu g_\mathrm{r}^2\stackrel{~}{\mathrm{\Pi }}(𝐪,i\nu _n\omega _\text{q}+i\delta ;𝐫)2\mu (𝐫).$$
(48)
One can understand the basic physics by going back to eq. (25) and noting that

$`1U_{\mathrm{eff}}(q)\mathrm{\Pi }`$ $`=`$ $`[1U\mathrm{\Pi }][1+g_\mathrm{r}^2D_0\stackrel{~}{\mathrm{\Pi }}]`$ (49)
$`=`$ $`[1U\mathrm{\Pi }]+g_r^2\mathrm{\Pi }{\displaystyle \frac{1}{\omega (E_\text{q}^0+2\nu 2\mu (𝐫))}}.`$ (50)
This form shows clearly that the zeros of $`1U_{\mathrm{eff}}(q)\mathrm{\Pi }`$ correspond to a hybridized spectrum associated with the preformed Cooper-pair \[$`1U\mathrm{\Pi }=0`$\] and the $`b`$-Bosons \[$`\omega (E_𝐪^0+2\nu 2\mu (𝐫))=0`$\] by the coupling $`g_\mathrm{r}^2\mathrm{\Pi }`$. The zeros of eq. (50) are seen to be given by

$$[1U\mathrm{\Pi }]\left[\omega [E_𝐪^0+2\nu 2\mu (𝐫)]\right]+g_\mathrm{r}^2\mathrm{\Pi }=0.$$
(51)
One sees that the same hybridization physics is given by the poles of $`\stackrel{~}{D}`$, i.e.,

$$[\omega (E_𝐪^0+2\nu 2\mu (𝐫))]+g_\mathrm{r}^2\frac{\mathrm{\Pi }}{1U\mathrm{\Pi }}=0.$$
(52)
All these versions emphasize that the renormalized Bosons given by $`\stackrel{~}{D}^1=0`$ are the hybridized version of the bare preformed Cooper-pairs and the bare $`b`$-Bosons. One can still ask what is the relative “weighting” of these two components in the renormalized Bose spectrum given by the pole of $`\stackrel{~}{D}`$, as given by the values of $`N_\mathrm{B}`$ and $`N_\mathrm{C}`$ in eq. (47). In Ref. , we calculated both these contributions as a function of the threshold $`2\nu `$ in the case of a uniform Fermi gas.

Finally, we note that one can write

$`1U_{\mathrm{eff}}(q)\mathrm{\Pi }`$ $`=`$ $`D_0\left[D_0^1U_{\mathrm{eff}}(q)D_0^1\mathrm{\Pi }\right]`$ (53)
$`=`$ $`D_0\left[D_0^1[{\displaystyle \frac{U}{D_0}}g_\mathrm{r}^2]\mathrm{\Pi }\right]`$ (54)
$``$ $`D_0D_\mathrm{M}^1.`$ (55)
Putting this into eq. (29), one finds that the term $`2N_\mathrm{B}^0`$ is cancelled out by the $`D_0`$ factor in eq. (55) and we are left with

$$N=𝑑𝐫[N_\mathrm{F}^0(𝐫)+2N_\mathrm{M}(𝐫)].$$
(56)
Here the number density of composite molecular Bosons $`N_\mathrm{M}(𝐫)`$ is defined as

$$2N_\mathrm{M}(𝐫)\frac{1}{\beta }\underset{\text{q},\nu _n}{}e^{i\nu _n\delta }\frac{}{\mu }\mathrm{ln}D_\mathrm{M}^1,$$
(57)
where $`D_\mathrm{M}`$ is defined in eq. (55),

$$D_\mathrm{M}^1(𝐪,i\nu _n;𝐫)i\nu _n[E_\text{q}^0+2\nu +(U/D_0g_\mathrm{r}^2)\mathrm{\Pi }]+2\mu (𝐫).$$
(58)
Note that $`\mathrm{\Pi }`$ enters here, not $`\stackrel{~}{\mathrm{\Pi }}`$ as in eq. (37). This composite molecule is a true hybrid of preformed Cooper-pairs and $`b`$-Bosons and cannot be associated with either component. However, one can verify that the poles of $`D_\mathrm{M}`$ are identical to those of $`\stackrel{~}{D}`$ as defined in eq. (37), in that \[see eqs. (50) and (55)\]

$`(1U\mathrm{\Pi })\stackrel{~}{D}^1=D_0^1({\displaystyle \frac{U}{D_0}}g_\mathrm{r}^2)\mathrm{\Pi }=D_\mathrm{M}^1.`$ (59)
## IV Transition Temperature in the BCS-BEC crossover with a Feshbach resonance

### A Renormalization for $`U`$ and $`g_\mathrm{r}`$

In this section, we consider the case of an isotopic harmonic trap ($`\omega _{0x}=\omega _{0y}=\omega _{0z}\omega _0`$) for simplicity. It is convenient to take $`\epsilon _\mathrm{F}`$ as the unit of energy (for one Fermi species). As for the unit of length, we take the radius $`R_\mathrm{F}`$ of a non-interacting Fermi gas at $`T=0`$ defined by $`\epsilon _\mathrm{F}(R_\mathrm{F})\epsilon _\mathrm{F}m\omega _0^2R_\mathrm{F}^2/2=0`$, or $`R_\mathrm{F}^2=2\epsilon _\mathrm{F}/m\omega _0^2`$. We also take, as the unit of the number of particles, the total number of Fermi atoms given by

$$N=\frac{4}{\pi }_0^{R_\mathrm{F}}r^2𝑑r_0^{p_\mathrm{F}(r)}p^2𝑑p=\frac{(m\omega _0)^3R_\mathrm{F}^6}{24},$$
(60)
where $`p_\mathrm{F}(r)=\sqrt{2m\epsilon _\mathrm{F}(r)}`$ is the Fermi momentum at $`r`$ within the LDA. With this normalization, since $`N`$ and $`R_\mathrm{F}`$ always appear with the coupling constants $`U`$ and $`g_\mathrm{r}`$ as $`\stackrel{~}{U}UN/R_\mathrm{F}^3`$ and $`\stackrel{~}{g}_\mathrm{r}^2g_\mathrm{r}^2N/R_\mathrm{F}^3`$, they can be absorbed into these coupling constants. (In the following, we simply write $`\stackrel{~}{U}`$ and $`\stackrel{~}{g}_\mathrm{r}`$ as $`U`$ and $`g_\mathrm{r}`$.) In our numerical calculations, we simply introduce a Gaussian cutoff $`e^{\epsilon _\text{p}^2/\omega _\mathrm{c}^2}`$, with $`\omega _\mathrm{c}=2\epsilon _\mathrm{F}`$, in the momentum summation involved in the particle-particle response function $`\mathrm{\Pi }(𝐪,i\nu _n;𝐫)`$ defined in eq. (15) and thus also in eq. (20).

The interactions $`U`$ and $`g_\mathrm{r}`$ used in this paper should be regarded as renormalized quantities effectively involving the effects of screening by (pseudo-spin and density) fluctuations and renormalization in terms of high energy processes: Gor’kov et al. showed within the mean-field BCS theory that in the dilute limit, $`T_\mathrm{c}`$ is suppressed by an effective interaction mediated by fluctuations described by particle-hole (p-h) bubble diagrams. This effect was also studied recently by Combescot in connection with the possibility of high $`T_\mathrm{c}`$ superfluidity in 6Li. This reduction of $`T_\mathrm{c}`$ can be physically understood as a screening effect by (pseudo-) spin and density fluctuations, which weakens the pairing interaction between Fermi atoms. One should understand that $`U`$ and $`g_\mathrm{r}`$ in our model in eq. (2) include this p-h screening.

As for the renormalization of $`U`$ and $`g_\mathrm{r}`$ by high energy processes, we note that our Fermi atomic gas has no physical cutoff-energy, in contrast to superconductivity in metals, where the lattice phonon Debye frequency acts the effective cutoff in the BCS gap equation. However, a prescription about the renormalization of the high energy interactions is implicit when we introduce a cutoff $`\omega _c`$ in the momentum summations. To see this, it is convenient to introduce an energy cutoff $`E_c`$. This cutoff energy is assumed to be much higher than $`\omega _c`$ and $`T`$, and it may be arbitrarily large. We can then divide the summation in the equation for $`T_\mathrm{c}`$ in eq. (20) into a low energy part $`\epsilon _\text{p}=[0,\omega _c]`$ and high energy part $`\epsilon _\text{p}=[\omega _c,E_c]`$. In this decomposition, when $`E_c\omega _c|\mu |>T_\mathrm{c}`$, we can take $`\mathrm{tanh}(\epsilon _\text{p}\mu )/2T_\mathrm{c}1`$ in the high energy part. Since this replacement corresponds to ignoring the influence of the Fermi distribution \[recall that in eq. (20), $`\mathrm{tanh}x/2T=12f(x)`$\], the high energy region can be described as two atoms interacting in a vacuum. The effect of the surrounding gas and Fermi statistics is absent. Then the decomposition of eq. (20) is reduced to

$$1=U_{\mathrm{eff}}\underset{[0,\omega _c]}{}\frac{\mathrm{tanh}(\epsilon _\text{p}\mu )/2T_\mathrm{c}}{2\epsilon _\text{p}2\mu }+U_{\mathrm{eff}}\underset{[\omega _c,E_c]}{}\frac{1}{2\epsilon _\text{p}2\mu },$$
(61)
where $`U_{\mathrm{eff}}=U+g_\mathrm{r}^2/(2\nu 2\mu )`$. This equation can be rewritten as the form

$`1`$ $`=`$ $`{\displaystyle \frac{U_{\mathrm{eff}}}{1U_{\mathrm{eff}}_{[\omega _c,E_c]}\frac{1}{2\epsilon _\text{p}2\mu }}}{\displaystyle \underset{[0,\omega _c]}{}}{\displaystyle \frac{\mathrm{tanh}(\epsilon _\text{p}\mu )/2T_\mathrm{c}}{2\epsilon _\text{p}2\mu }}`$ (62)
$``$ $`U_{\mathrm{eff}}^R(\omega _c){\displaystyle \underset{[0,\omega _c]}{}}{\displaystyle \frac{\mathrm{tanh}(\epsilon _\text{p}\mu )/2T_\mathrm{c}}{2\epsilon _\text{p}2\mu }}.`$ (63)
In this equation for $`T_\mathrm{c}`$, high energy two-particle scattering processes in the region $`[\omega _c,E_c]`$ has been absorbed into the renormalized coupling $`U_{\mathrm{eff}}^R(\omega _c)=U^R(\omega _c)+g_{\mathrm{r}}^{R}{}_{}{}^{2}(\omega _c)/(2\nu ^R(\omega _c)2\mu )`$ describing the low energy physics in the region $`[0,\omega _c]`$. The renormalized interaction $`U^R(\omega _c)`$, Feshbach coupling $`g_\mathrm{r}^R(\omega _c)`$ and threshold energy $`2\nu ^R(\omega _c)`$ are thus given by

$`\{\begin{array}{c}U^R(\omega _c)={\displaystyle \frac{U}{1U_{[\omega _c,E_c]}\frac{1}{2\epsilon _\text{p}2\mu }}},\hfill \\ g_\mathrm{r}^R(\omega _c)={\displaystyle \frac{g_\mathrm{r}}{1U_{[\omega _c,E_c]}\frac{1}{2\epsilon _\text{p}2\mu }}},\hfill \\ 2\nu ^R(\omega _c)=2\nu g_r^2{\displaystyle \frac{\underset{[\omega _c,E_c]}{}\frac{1}{2\epsilon _\text{p}2\mu }}{1U_{[\omega _c,E_c]}\frac{1}{2\epsilon _\text{p}2\mu }}}.\hfill \end{array}`$ (67)
Physically, $`U^R(\omega _c)`$ is the non-resonant Fermion-Fermion interaction enhanced by high energy scattering processes and is diagrammatically described by Fig. 3 (a). The renormalized Feshbach coupling $`g_\mathrm{r}^R(\omega _c)`$ involves the three point vertex correction coming from the high-energy processes in Fig. 3(b). Finally, the renormalized threshold energy $`2\nu ^R(\omega _c)`$ originates from the self-energy correction caused by the break of a Bose molecule into two Fermi atoms with high momenta, which can be described by the diagrams in Fig. 3(c).

We should determine the renormalized quantities experimentally. In this regard, the renormalized quantities for the low energy limit ($`\omega _c=0`$) can be introduced. One can write the renormalized parameters in eq. (67) in terms of their low energy values ($`\omega _c=0`$) as follows:

$`\{\begin{array}{c}U^R(\omega _c)={\displaystyle \frac{U^R(0)}{1+U^R(0)_{[0,\omega _c]}\frac{1}{2\epsilon _\text{p}2\mu }}},\hfill \\ g_\mathrm{r}^R(\omega _c)={\displaystyle \frac{g_\mathrm{r}^R(0)}{1+U^R(0)_{[0,\omega _c]}\frac{1}{2\epsilon _\text{p}2\mu }}},\hfill \\ 2\nu ^R(\omega _c)=2\nu ^R(0)+g_r^R(0)^2{\displaystyle \frac{\underset{[0,\omega _c]}{}\frac{1}{2\epsilon _\text{p}2\mu }}{1+U^R(0)_{[0,\omega _c]}\frac{1}{2\epsilon _\text{p}2\mu }}}.\hfill \end{array}`$ (71)
Written in terms of the renormalized variables given in eq. (71), the BCS gap equation in eq. (63) no longer involves $`E_c`$ explicitly. Although we simply write the couplings as $`U`$, $`g_\mathrm{r}`$ and the threshold energy as $`2\nu `$ in our model, we can regard them as renormalized quantities incorporating the high energy two-particle scattering processes, as discussed above.

The present way of renormalization can reproduce the cutoff-free theories in Refs.. When we write eq. (63) using the renormalized variables in terms of the low energy value $`U_{\mathrm{eff}}^R(\omega _\mathrm{c}=0)`$ as given by eq. (71), we obtain

$$1=\frac{U_{\mathrm{eff}}^R(0)}{1+U_{\mathrm{eff}}^R(0)_{[0,\omega _c]}\frac{1}{2\epsilon _\text{p}2\mu }}\underset{[0,\omega _c]}{}\frac{\mathrm{tanh}(\epsilon _\text{p}\mu )/2T_\mathrm{c}}{2\epsilon _\text{p}2\mu }.$$
(72)
This equation can be in turn rewritten in the form

$`1`$ $`=`$ $`U_{\mathrm{eff}}^R(0){\displaystyle \underset{[0,\omega _c]}{}}\left[{\displaystyle \frac{\mathrm{tanh}(\epsilon _\text{p}\mu )/2T_\mathrm{c}}{2\epsilon _\text{p}2\mu }}{\displaystyle \frac{1}{2\epsilon _\text{p}2\mu }}\right]`$ (73)
$``$ $`U_{\mathrm{eff}}^R(0){\displaystyle \underset{[0,\mathrm{}]}{}}\left[{\displaystyle \frac{\mathrm{tanh}(\epsilon _\text{p}\mu )/2T_\mathrm{c}}{2\epsilon _\text{p}2\mu }}{\displaystyle \frac{1}{2\epsilon _\text{p}2\mu }}\right].`$ (74)
Since the two terms in the square brackets almost cancel out each other in the high energy region $`[\omega _c,\mathrm{}]`$ (we assume that $`\omega _cT`$), we can eliminate the finite upper limit $`\omega _c`$ in the summation. This cutoff-free expression is similar to the one used in Ref.. The cutoff-free BCS equation for $`T_\mathrm{c}`$ used by Randeria et al. can be also obtained from eq. (20). Writing it as

$$1=U_{\mathrm{eff}}\underset{[0,E_c]}{}\left[\frac{\mathrm{tanh}(\epsilon _\text{p}\mu )/T_\mathrm{c}}{2\epsilon _\text{p}\mu }\frac{1}{2\epsilon _\text{p}}+\frac{1}{2\epsilon _\text{p}}\right],$$
(75)
we can turn this into

$$1=\overline{U}_{\mathrm{eff}}\underset{[0,E_c]}{}\left[\frac{\mathrm{tanh}(\epsilon _\text{p}\mu )/T_\mathrm{c}}{2\epsilon _\text{p}\mu }\frac{1}{2\epsilon _\text{p}}\right],$$
(76)
where the renormalized interaction $`\overline{U}_{\mathrm{eff}}=\overline{U}+\overline{g}_\mathrm{r}^2/(2\overline{\nu }2\mu )`$ is defined by

$`\{\begin{array}{c}\overline{U}={\displaystyle \frac{U}{1U_{[0,E_c]}\frac{1}{2\epsilon _\text{p}}}},\hfill \\ \overline{g}_\mathrm{r}=g_\mathrm{r}{\displaystyle \frac{1}{1U_{[0,E_c]}\frac{1}{2\epsilon _\text{p}}}},\hfill \\ 2\overline{\nu }=2\nu g_\mathrm{r}^2{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{\underset{[0,E_c]}{}\frac{1}{2\epsilon _\text{p}}}{1U_{[0,E_c]}\frac{1}{2\epsilon _\text{p}}}}.\hfill \end{array}`$ (80)
All these expressions for the equation for $`T_\mathrm{c}`$ give physically the same result, assuming that one takes the appropriate values for renormalized parameters.

With regard to the renormalization of the coupling constants, recent work by Kokkelmans et al. has discussed the coupled Fermion-Boson model given in eq. (2) and incorporated the full two-body scattering theory into an improved fashion. Milstein et al. have more recently derived identical equations to eqs. (4) and (6) in Ref.  (these correspond to eqs. (20) and (29) for a trapped Fermi gas), which include the effect of particle-particle fluctuations and the appearance of Cooper-pairs. They have extended our analysis in Ref.  by renormalizing these equations in a way which treats the two-body scattering and bound states correctly. Their work thus leads to explicit numerical values for the various parameters in the renormalized model corresponding to eq. (2).

In the rest of this Section and Sections V and VI, we discuss our results for $`T_\mathrm{c}`$ and the region above, using $`U=0.3\epsilon _\mathrm{F}`$ and $`g_\mathrm{r}=0.6\epsilon _\mathrm{F}`$. Our diagrammatic approximation for the fluctuations, as summarized in Figs. 1 and 2, assumes that the interactions are weak. In Section VII, we discuss the analogous results predicted by our equations (20) and (29) using a value $`g_\mathrm{r}\epsilon _\mathrm{F}`$, such as found in Refs. . In this case, the Feshbach resonance is very broad and thus can have a huge effect on $`T_\mathrm{c}`$ even if the threshold $`2\nu 2\epsilon _\mathrm{F}`$.

### B Crossover behavior at $`T_\mathrm{c}`$ within LDA

Our numerical solution for $`T_\mathrm{c}`$ and $`\mu (T_\mathrm{c},N)`$ from eqs. (20) and (29) for a trapped Fermi gas gives results which are very similar to those for a uniform Fermi gas obtained in Ref. . Figure 4 shows the BCS-BEC crossover in a trapped Fermi gas. As the threshold energy $`2\nu `$ is lowered, Fig. 4(a) shows that $`T_\mathrm{c}`$ deviates from the mean-field BCS theory (‘BCS’ in the figure) and approaches the BEC phase transition temperature of a gas of $`N/2`$ free Bosons of mass $`M=2m`$ (‘BEC’ in the figure). The chemical potential $`\mu `$ also changes remarkably in this crossover as shown in Fig. 4(b). As expected from the mean-field BCS theory, we obtain $`\mu \epsilon _\mathrm{F}`$ for $`2\nu >2\epsilon _\mathrm{F}`$ (BCS regime) while $`\mu `$ deviates from $`\epsilon _\mathrm{F}`$ and approaches $`\nu `$ for $`\nu <0`$ (BEC regime). The latter result reflects that the $`b`$-Bosons become the dominant contribution in the BEC regime. As a result, the phase transition occurs when the chemical potential of the $`b`$-Boson $`\mu _\mathrm{B}`$ ($`=2\mu `$) reaches the bottom of the Bose band ($`2\nu `$) at the center of the trap (within our LDA). We note that we can rewrite the equation for $`T_\mathrm{c}`$ in eq. (18) in the form

$$2\nu =2\mu +g_\mathrm{r}^2\stackrel{~}{\mathrm{\Pi }}(0,0,𝐫=0).$$
(81)
This equation reduces to $`2\nu =2\mu `$ $`(=\mu _\mathrm{B})`$ in the BEC limit, where the Fermi atoms and hence the particle-particle fluctuations described by $`\mathrm{\Pi }`$ are absent. Since, in this BEC regime, the dominant excitations are $`N/2`$ free Bosons, we find the usual result for a trapped free Bose gas $`T_\mathrm{c}=\omega _0[N/2\zeta (3)]^{1/3}`$, where $`\omega _0`$ is the trap frequency and $`\zeta (x)`$ is the usual zeta function (we set $`\mathrm{}=1`$). Using the fact that $`T_\mathrm{F}=(6N_{})^{1/3}\omega _0`$ is the Fermi temperature for a trapped free Fermi gas (with $`N_{}=N/2`$), we can write this expression for $`T_\mathrm{c}`$ in terms of $`T_\mathrm{F}`$ to give

$$T_\mathrm{c}=\left(\frac{1}{6\zeta (3)}\right)^{1/3}\epsilon _\mathrm{F}=0.518T_\mathrm{F}.$$
(82)
Figure 4(a) shows that this $`T_\mathrm{c}`$ is obtained in the BEC regime as the limiting maximum value in the trapped Fermi gas. We recall that Fig. 4 is for $`U=0.3\epsilon _\mathrm{F}`$ and $`g_\mathrm{r}=0.6\epsilon _\mathrm{F}`$.

We note that this upper limit for $`T_\mathrm{c}`$ is much higher than that of a uniform Fermi gas in the absence of a trap. In the uniform system, $`T_\mathrm{c}`$ in the BEC regime is given by

$$T_\mathrm{c}=2\left(\frac{1}{6\sqrt{\pi }\zeta (3/2)}\right)^{2/3}\epsilon _\mathrm{F}=0.218T_\mathrm{F},$$
(83)
where the Fermi temperature in the uniform system is given by $`T_\mathrm{F}=(6\pi ^2N_{})^{2/3}/2m`$. One reason for this different maximum $`T_\mathrm{c}`$’s between the two cases comes from the different definition of $`T_\mathrm{F}`$. In addition, one should note that the harmonic trap potential enhances $`T_\mathrm{c}`$ through the effective density of states (DOS) of Bose molecules. When we carry out the spatial integration in the equation for the total number of trapped particles in the BEC regime \[$`N=2𝑑𝐫N_\mathrm{B}(𝐫)`$\], this equation is reduced to

$$1=3_0^{\mathrm{}}𝑑EE^2\frac{1}{e^{\overline{\beta }E}1}(T=T_\mathrm{c}),$$
(84)
where $`\overline{\beta }T_\mathrm{F}/T`$ and the energy in the integral is normalized in units of $`\epsilon _\mathrm{F}`$. The corresponding equation in the uniform system is given by

$$1=3\sqrt{2}_0^{\mathrm{}}𝑑E\sqrt{E}\frac{1}{e^{\overline{\beta }E}1}.$$
(85)
This shows that the low energy region of the DOS has larger weight in the uniform gas (DOS $`\sqrt{E}`$) than in the trapped one (DOS $`E^2`$). As a result, the transition temperature (relative to $`T_\mathrm{F}`$ for each case) must be higher in the trapped gas than in the uniform gas.

## V Long-lived Composite Molecular Bosons

### A Damping of $`b`$-Boson and condition for stable molecules at $`T_\mathrm{c}`$

The decay of $`b`$-Boson into two Fermi atoms associated with the Feshbach resonance is described by the imaginary part of the (analytic-continued) self-energy of the renormalized Bose Green function in eq. (37)

$$\gamma Im\mathrm{\Sigma }(𝐪,i\nu _n\omega +i\delta ;𝐫)=g_\mathrm{r}^2Im\stackrel{~}{\mathrm{\Pi }}(𝐪,\omega +i\delta ;𝐫).$$
(86)
Since $`Im\mathrm{\Pi }(𝐪,\omega +i\delta ;𝐫)`$ is given by

$$Im\mathrm{\Pi }(𝐪,\omega +i\delta ;𝐫)=\frac{m\sqrt{m}}{4\pi \beta \sqrt{E_\text{q}^0}}\mathrm{\Theta }(\omega +2\mu (𝐫)E_\text{q}^0)\mathrm{ln}\frac{\mathrm{cosh}\frac{\beta }{2}(\frac{\omega }{2}+\sqrt{(\omega +2\mu (𝐫)E_\text{q}^0)E_\text{q}^0})}{\mathrm{cosh}\frac{\beta }{2}(\frac{\omega }{2}\sqrt{(\omega +2\mu (𝐫)E_\text{q}^0)E_\text{q}^0})},$$
(87)
where $`\mathrm{\Theta }(x)`$ is the step function, the damping rate $`\gamma `$ is finite only in the region $`\omega \omega _{\mathrm{th}}E_\text{q}^02\mu (𝐫)`$. This step function restriction is clear when we note that the denominator of eq. (15) can be written as $`[i\nu _n+2\mu (𝐫)E_\text{q}^0]2\epsilon _\text{p}`$. At $`𝐪=0`$, eq. (87) can be simplified to

$$Im\mathrm{\Pi }(𝐪=0,\omega +i\delta ;𝐫)=\frac{m\sqrt{m}}{4\pi }\mathrm{\Theta }(\omega +2\mu (𝐫))\sqrt{\omega +2\mu (𝐫)}\mathrm{tanh}\frac{\beta \omega }{4}.$$
(88)
Since the energy to dissociate one $`b`$-Boson (with center of mass momentum $`𝐪`$) into two Fermi atoms is given by $`[\epsilon _{\text{p}+\text{q}/2}\mu (𝐫)]+[\epsilon _{\text{p}+\text{q}/2}\mu (𝐫)]`$, the threshold energy $`\omega _{\mathrm{th}}`$ is found to correspond to the minimum energy to destroy one $`b`$-Boson with $`𝐪=0`$.

The damping rate $`\gamma `$ leads to finite lifetime of $`b`$-Boson, or equivalently gives a finite width to the peak in the spectral density of the (renormalized) Bose excitation spectrum given by

$$\rho _\mathrm{B}(𝐪,\omega ;𝐫)\frac{1}{\pi }Im\stackrel{~}{D}(𝐪,i\nu _n\omega +i\delta ;𝐫).$$
(89)
In the BCS regime, where $`\gamma `$ is always finite for $`\omega >0`$ (see the inset in Fig. 5(a)), the quasi-particle $`b`$-Bose spectral density shows a broad peak, as shown in Fig. 5(a). This $`b`$-Boson thus has a finite lifetime, which is given by the inverse of the peak width, due to decay into two Fermi atoms by the coupling to the Feshbach resonance. As one approaches the BEC regime, $`\omega _{\mathrm{th}}>0`$ is realized (as shown in the insets in Figs. 5(b) and 5(c)). In this case, when the $`b`$-Boson pole given by eq. (48),

$$\omega _\text{q}=\left[E_\text{q}^02\mu (𝐫)\right]+\left[2\nu g_\mathrm{r}^2\stackrel{~}{\mathrm{\Pi }}(𝐪,\omega _\text{q}+i\delta ;𝐫)\right],$$
(90)
appears below $`\omega _{\mathrm{th}}`$, its lifetime is infinite. It thus appears as a $`\delta `$-function peak in the $`b`$-Boson excitation spectral density $`\rho _\mathrm{B}`$ given by eq. (89). We show this case in Figs. 5(b) and 5(c). In Fig. 5(b), the excitation spectrum is still accompanied by a incoherent part in high energy region ($`\omega >\omega _{\mathrm{th}}`$), in addition to the undamped bound state pole below $`\omega _{\mathrm{th}}`$. On the other hand, the incoherent part is seen to be almost absent in the BEC regime in Fig. 5(c), evidence that a stable $`b`$-molecule dominates the two-particle excitation spectrum.

The condition for such a stable $`b`$-Boson is simply that

$$2\stackrel{~}{\nu }2\nu g_\mathrm{r}^2\stackrel{~}{\mathrm{\Pi }}(𝐪,\omega _\text{q}+i\delta ;𝐫)<0[E_\text{q}^02\mu (𝐫)>0].$$
(91)
This means that the renormalized threshold energy defined by $`2\stackrel{~}{\nu }`$ in eq. (91) must be negative for the stable undamped $`b`$-Bose molecules to appear. In the region $`\omega =𝐪=0`$, this condition can be approximately written as $`2\stackrel{~}{\nu }2\nu g_\mathrm{r}^2\stackrel{~}{\mathrm{\Pi }}(0,0;𝐫)<0`$. Since $`\mathrm{\Pi }(0,0,𝐫)`$ can be shown to be positive and also decreases as $`𝐫`$ increases, the condition $`2\stackrel{~}{\nu }<0`$ is first satisfied at the center of the trap.

We can divide the number $`N_\mathrm{B}`$ of renormalized $`b`$-Bosons in eq. (47) into the contribution from stable long-lived part ($`N_\mathrm{B}^{\gamma =0}`$) and a Feshbach resonating contribution with a finite lifetime (denoted by $`N_\mathrm{B}^{\gamma >0}`$). In order to evaluate the contribution from the stable part $`N_\mathrm{B}^{\gamma =0}`$, it is convenient to introduce the spectral representation of the renormalized Bose Green function as

$$\stackrel{~}{D}(𝐪,i\nu _n;𝐫)=\frac{1}{\pi }_{\mathrm{}}^{\mathrm{}}dz\frac{1}{i\nu _nz}Im\stackrel{~}{D}(𝐪,i\nu _nz+i\delta ;𝐫).$$
(92)
Using this expression, the summation in terms of the Boson Matsubara frequency in $`N_\mathrm{B}`$ in eq. (47) can be carried out to give

$$N_\mathrm{B}=\frac{1}{\pi }\underset{\text{q}}{}𝑑𝐫_{\mathrm{}}^{\mathrm{}}𝑑zn_B(z)Im\stackrel{~}{D}(𝐪,z+i\delta ;𝐫).$$
(93)
Using this, the contribution from the stable poles (with no imaginary part) is given by

$`N_\mathrm{B}^{\gamma =0}`$ $`=`$ $`{\displaystyle 𝑑𝐫\underset{\text{q}}{\overset{\mathrm{poles}}{}}_{\mathrm{}}^{\mathrm{}}𝑑z\delta (z[E_\text{q}2\mu (𝐫)][2\mu g_\mathrm{r}^2\stackrel{~}{\mathrm{\Pi }}(𝐪,z+i\delta ;𝐫)])}`$ (94)
$`=`$ $`{\displaystyle 𝑑𝐫\underset{\text{q}}{\overset{\mathrm{poles}}{}}Z(𝐪;𝐫)n_\mathrm{B}(\omega _𝐪)}.`$ (95)
Here, $`Z(𝐪;𝐫)`$ arises from what is called mass renormalization and is given by

$$Z(𝐪;𝐫)^11+g_\mathrm{r}^2\frac{\stackrel{~}{\mathrm{\Pi }}(𝐪,\omega _𝐪,𝐫)}{\omega _\text{q}}.$$
(96)
The sum in eq. (95) is over the solutions $`\omega _\text{q}`$ satisfying eq. (90).

Besides this Feshbach component of the two-particle bound state given by eq. (95), there is a component from Cooper-pair associated with $`N_\mathrm{C}`$ as defined in eq. (46). One can rewrite eq. (46) in the form

$$2N_\mathrm{C}=\frac{1}{\beta }\underset{\text{q},\nu _n}{}𝑑𝐫\frac{\stackrel{~}{D}}{D_0}\frac{U_{\mathrm{eff}}(q)}{1U\mathrm{\Pi }}\frac{\mathrm{\Pi }}{\mu }.$$
(97)
We next use the fact that

$$\frac{\stackrel{~}{\mathrm{\Pi }}}{\mu }=\left(\frac{1}{1U\mathrm{\Pi }}\right)^2\frac{\mathrm{\Pi }}{\mu },$$
(98)
and that for $`\omega =\omega _\text{q}`$ \[given by eq. (90)\], one has the identities $`D_0^1=g_\mathrm{r}^2\stackrel{~}{\mathrm{\Pi }}`$ and $`U_{\mathrm{eff}}(q)=\mathrm{\Pi }^1`$. With these relations, the contribution $`N_\mathrm{C}^{\gamma =0}`$ of undamped Bosons associated with Cooper pairing can be reduced to the expression

$$N_\mathrm{C}^{\gamma =0}=\frac{g_\mathrm{r}^2}{2}𝑑𝐫\underset{\text{q}}{\overset{\mathrm{poles}}{}}\frac{\stackrel{~}{\mathrm{\Pi }}(𝐪,\omega _\text{q};𝐫)}{\mu }Z(𝐪;𝐫)n_\mathrm{B}(\omega _\text{q}),$$
(99)
where $`Z(𝐪;𝐫)`$ has been defined in eq. (96). There is also a contribution from scattering states and hence we have $`N_\mathrm{C}=N_\mathrm{C}^{\gamma =0}+N_\mathrm{C}^{\mathrm{sc}}`$, where $`N_\mathrm{C}^{\mathrm{sc}}`$ represents the scattering contribution associated with damped particle-particle (or Cooper channel) fluctuations.

Figure 4(c) shows the change of the character of the excitations in the BCS-BEC crossover in the trap. As $`\nu `$ is lowered and approaches $`\epsilon _\mathrm{F}`$, the Feshbach resonance molecular state starts to play a role and the Feshbach damped molecules ($`N_\mathrm{B}^{\gamma >0}`$) appear. These Bosons are replaced by the stable molecules ($`N_\mathrm{B}^{\gamma =0}`$) around $`\nu =0.18\epsilon _\mathrm{F}`$ for the parameters used in Fig. 4(c). At the same time, a contribution from preformed Cooper-pairs ($`N_\mathrm{C}^{\gamma =0}`$) appear. In the BEC regime, stable Feshbach molecules become dominant and $`N_\mathrm{B}^{\gamma =0}`$ approaches $`0.5`$ ($`=N/2`$). The contribution from preformed Cooper-pairs are seen to again disappear in the BEC limit $`\nu <0`$. The increase in the total number of molecules is, however, slower than that found in the uniform system. This is simply because the stable molecules initially only begin to appear at the center of the trap, rather than throughout the system as in the uniform gas case..

Figure 6 shows the number of composite Bosons $`N_\mathrm{M}𝑑𝐫N_\mathrm{M}(𝐫)`$ as a function of $`\nu `$ \[where $`N_\mathrm{M}(𝐫)`$ is defined in eq. (57)\]. Although each component of Boson ($`N_\mathrm{B}^{\gamma =0}`$, $`N_\mathrm{B}^{\gamma >0}`$, $`N_\mathrm{C}^{\gamma =0}`$ and $`N_\mathrm{C}^{\mathrm{sc}}`$) shows a singular behavior around $`\nu =0.18\epsilon _\mathrm{F}`$ in Fig. 4(c), $`N_\mathrm{M}`$ ($`N_\mathrm{B}^{\gamma =0}+N_\mathrm{B}^{\gamma >0}+N_\mathrm{C}^{\gamma =0}+N_\mathrm{C}^{\mathrm{sc}}`$) is found to increase smoothly as the threshold energy $`2\nu `$ is lowered, as expected.

As we have discussed in detail in Section III, the strong hybridization induced by the Feshbach coupling interaction $`g_\mathrm{r}`$ means that there is only a single two-particle Bose branch, whose energy is given by the solution of eq. (90). Thus $`N_\mathrm{B}^{\gamma =0}`$ and $`N_\mathrm{C}^{\gamma =0}`$ only represent the relative spectral weights of the two components (corresponding to long-lived $`b`$-Bosons and stable preformed Cooper-pairs) of this composite Boson. This decomposition is only meant as a way of understanding the physics of the BCS-BEC crossover which arises in the model given by eq. (2). This interpretation is, however, very useful to understand the difference between the present BCS-BEC crossover and that in strong-coupling superconductivity with no Feshbach resonance: In the latter case, the preformed Cooper-pairs are the dominant excitations in the strong coupling BEC regime. In the present case, for the values of the parameters $`U=0.3\epsilon _\mathrm{F}`$ and $`g_\mathrm{r}=0.6\epsilon _\mathrm{F}`$ used in Fig. 4, although preformed Cooper-pairs make a contribution, the stable Feshbach molecules dominate the excitation spectrum in the BEC limit ($`\nu <0`$).

### B The binding energy of composite Bosons

The binding energy of the long-lived composite Boson is the energy to break this molecule into two Fermi atoms. This dissociation energy can be evaluated from the energy difference between the bound state energy $`\omega _𝐪`$ in eq. (90) and the excitation spectrum of the free Fermi gas measured from the chemical potential $`\mu `$. In particular, the binding energy of the zero-energy ($`\omega _{\text{q}=0}=0`$) two-particle bound state with $`𝐪=0`$ which appears at $`T_\mathrm{c}`$ is given by

$`E_{\mathrm{binding}}=2[\epsilon _{\text{q}=0}\mu (𝐫=0)]\omega _{\text{q}=0}`$ $`=`$ $`2|\mu |`$ (100)
$`=`$ $`2|\stackrel{~}{\nu }(𝐪=0,\omega _\text{q}=0;𝐫=0)|,`$ (101)
where we have used $`\mu (𝐫=0)<0`$ and eq. (91) to obtain expressions on the right hand side. This result means that $`E_{\mathrm{binding}}`$ is equal to the absolute value of the energy of a renormalized $`b`$-molecule at the bottom of the renormalized Bose spectrum.

The threshold energy $`\omega _{\mathrm{th}}(𝐪=0)`$ of the two-particle spectrum is defined as where $`Im\mathrm{\Pi }(𝐪=0,\omega ;𝐫=0)`$ given in eq. (87) becomes finite (see Fig. 5). We recall that $`Im\mathrm{\Pi }`$ describes the excitation spectrum of the Cooper-pair channel fluctuations. This frequency $`\omega _{\mathrm{th}}(𝐪=0)`$ is the same as the threshold energy of the continuum spectrum of the renormalized $`b`$-Boson self-energy in eq. (35), since $`\gamma =Im\mathrm{\Sigma }=g_\mathrm{r}^2Im\stackrel{~}{\mathrm{\Pi }}(i\nu _n\omega +i\delta )Im\mathrm{\Pi }(\omega +i\delta )`$. When a stable $`b`$-Boson at $`\omega _{\text{q}=0}=0`$ is excited to a state with energy $`E\omega _{\mathrm{th}}`$, this excited $`b`$-Boson will decay into two Fermi atoms lading to a finite lifetime for this excited state. We conclude that $`\omega _{\mathrm{th}}(𝐪=0)=E_{\mathrm{binding}}`$, i.e., the particle-particle threshold energy $`\omega _{\mathrm{th}}(𝐪=0)`$ is equal to the binding energy of a stable two-particle bound state.

Since the chemical potential $`\mu `$ approaches $`\nu `$ in the BEC limit as shown in Fig. 4(b), we find from eq. (101) that $`E_{\mathrm{binding}}2|\nu |`$ in this regime. This energy $`2|\nu |`$ is just equal to the energy needed to transfer a free $`b`$-Boson at the bottom of the Bose energy band ($`E_{\text{q}=0}^0+2\nu =2|\nu |`$) to two Fermi atoms with zero energy ($`2\epsilon _{\text{q}=0}=0`$). This result is consistent with the fact that the preformed Cooper-pair component in the bound state is absent in the BEC limit, leaving the $`b`$-molecule component.

The binding energy of the preformed Cooper-pair in the strong-coupling BEC regime has been also discussed in superconductivity, and it is useful to understand how the present results are related to this previous work. To see this, let us recall the cutoff-free equation for $`T_\mathrm{c}`$ given by eq. (76). Since $`\mu 0`$ in the BEC limit, one can approximately take $`\mathrm{tanh}(\epsilon _\text{p}\mu )/2T=1`$. Summing over $`𝐩`$ in eq. (76), we obtain

$$1=\left(\overline{U}+\overline{g}_\mathrm{r}^2\frac{1}{2\overline{\nu }2\mu }\right)\frac{m\sqrt{2m|\mu |}}{4\pi }.$$
(102)
When we write the renormalized attractive interaction as $`\overline{U}_{\mathrm{eff}}\overline{U}+\overline{g}_\mathrm{r}^2/(2\overline{\nu }2\mu )4\pi a_s/m`$, where $`a_s`$ is an $`s`$-wave scattering length, we find $`|\mu |=1/2ma_s^2`$, which is formally the same expression as that obtained in strong-coupling superconductivity discussions. In the present model, where the pairing interaction mediated by the Feshbach resonance is dominant in the BEC regime, we can neglect the non-resonant part $`\overline{U}`$ in $`\overline{U}_{\mathrm{eff}}`$. Then in the case of $`\overline{\nu }<0`$, the solution of eq. (102) is

$$\frac{\mu }{\overline{\nu }}=1+\frac{m^3\overline{g}_\mathrm{r}^4}{64\pi ^2\overline{\nu }}\left[\sqrt{1+\frac{128\pi ^2|\overline{\nu }|}{m^3\overline{g}_\mathrm{r}^4}}1\right].$$
(103)
The right hand side of this equation is reduced to unity in the BEC limit, $`\overline{\nu }\mathrm{}`$, namely, the chemical potential for the $`b`$-Boson ($`2\mu `$) approaches the (renormalized) threshold energy $`2\overline{\nu }`$, as one expects.

The relation between $`\mu `$ and the binding energy of a bound state in a two-particle system can be also obtained by extending the previous discussions given in the superconductivity case. When we consider two Fermi atoms in a vacuum and employ relative coordinates, the scattering $`t`$-matrix for one Fermi atom with reduced mass $`m/2`$ and energy $`\omega `$ is given by

$`T(\omega )`$ $`=`$ $`\left(U+g_\mathrm{r}^2{\displaystyle \frac{1}{\omega 2\nu }}\right)+\left(U+g_\mathrm{r}^2{\displaystyle \frac{1}{\omega 2\nu }}\right){\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{1}{\omega 2\epsilon _\text{p}}}\left(U+g_\mathrm{r}^2{\displaystyle \frac{1}{\omega 2\nu }}\right)+\mathrm{}`$ (104)
$`=`$ $`{\displaystyle \frac{U+g_\mathrm{r}^2\frac{1}{2\nu \omega }}{1\left(U+g_\mathrm{r}^2\frac{1}{2\nu \omega }\right)_\text{p}\frac{1}{2\epsilon _\text{p}\omega }}}.`$ (105)
The $`𝐪=0`$ two-particle bound state is given by the pole of $`T(\omega )`$, namely

$$1=U_{\mathrm{eff}}^{}\underset{\text{p}}{}\frac{1}{2\epsilon _\text{p}\omega },$$
(106)
where $`U_{\mathrm{eff}}^{}U+g_\mathrm{r}^2/(2\nu \omega )`$. When we renormalize the interaction in eq. (106) as we have done to obtain eq. (76), eq. (106) is reduced to

$$1=\overline{U}_{\mathrm{eff}}^{}\underset{\text{p}}{}\left[\frac{1}{2\epsilon _\text{p}\omega }\frac{1}{2\epsilon _\text{p}}\right],$$
(107)
where $`\overline{U}_{\mathrm{eff}}^{}\overline{U}+\overline{g}_\mathrm{r}^2/(2\overline{\nu }\omega )`$, in which the renormalized variables are given in eq. (80). This equation is the same as eq. (105) in the BEC limit ($`\mu 0`$, where $`\mathrm{tanh}(\epsilon _\text{p}\mu )/2T1`$) when we replace $`\omega 2\mu `$. Thus the energy of the bound state is related to the chemical potential by the relation $`\omega =2|\mu |`$, which is consistent with our earlier result in eq. (101).

## VI Strong Coupling Effect above $`T_\mathrm{c}`$

### A BCS-BEC crossover effect on the atomic density profile

As expected, the change of the character of the particles from the Fermion to Boson in the BCS-BEC crossover shown in Fig. 4(c) strongly affects the density profile of atoms in the trap. This is shown in Fig. 7(a). In the BCS regime, the density profile is spread out, as shown by result for $`\nu =\epsilon _\mathrm{F}`$. In this regime, the density profile is mainly composed of Fermi atoms, and the Pauli exclusion principle between Fermi atoms effectively acts as a repulsive interaction. This pushes the atoms away from the center of the trap. As the system approaches the BEC regime ($`\nu <0`$), the number of stable Bosons gradually increases and thus the repulsive effect of the Pauli exclusion principle becomes less important. As expected for Bosons, the particles are now seen to cluster near the center of the trap as $`\nu `$ is lowered, as shown by Fig. 7(a) \[see also the inset in Fig. 7(a) for $`\nu =0`$\]. These results indicate that the easily measured atomic density profile may be a very useful signature in looking for the BCS-BEC crossover phenomenon experimentally, especially as one gets close to $`T_\mathrm{c}`$.

We note that the atomic density profile shows a steep decrease near the trap center at $`T_\mathrm{c}`$ (Fig. 7(a)), but this is absent at $`T=1.5T_\mathrm{c}`$ (Fig. 7(b)). The cusp-like behavior at $`T_\mathrm{c}`$ originates from the Boson component ($`N_\mathrm{B}`$ and $`N_\mathrm{C}`$) in the density, as shown in the inset in Fig. 7(a). In the extreme BEC limit, since the gas is described by $`N/2`$ free Bose gas in a trap, the density of atoms is proportional to $`g_{3/2}(z(𝐫))`$, where $`g_{\frac{3}{2}}(z)_{n=1}^{\mathrm{}}z^n/n^{3/2}`$ and the fugacity is $`z(𝐫)e^{(2\mu 2\nu 2V_{\mathrm{trap}}(𝐫))/T}`$. At $`T_\mathrm{c}`$, since the BEC is realized when the chemical potential of Bosons ($`2\mu `$) reaches the bottom of the Boson excitation spectrum ($`2\nu `$), we find $`z(𝐫)1`$ at the trap center. From the well known behavior of $`g_{3/2}(z)`$ at $`z=1`$, the atomic density profile shows a finite slope at $`𝐫=0`$ at $`T_\mathrm{c}`$. In contrast, this sharp cusp at the trap center is absent for temperatures above $`T_\mathrm{c}`$, where $`2\mu <2\nu `$ and hence the fugacity $`z(𝐫)`$ is no longer close to unity. We also show in the inset in Fig. 7(b) the density profile above $`T_\mathrm{c}`$ when the temperature is fixed as $`T=0.75T_\mathrm{F}`$. In this case, the density profile becomes more spatially diffuse as $`\nu `$ is increased. This is because of the increasing dominance of the Fermions and the resulting enhancement of the effect of the Pauli exclusion principle.

### B Crossover from a Fermi gas to a Bose gas above $`T_\mathrm{c}`$

It is useful to clarify when and how the stable Bosons ($`N_\mathrm{B}^{\gamma =0}`$ and $`N_\mathrm{C}^{\gamma =0}`$) at $`T_\mathrm{c}`$ shown in Fig. 4(c) disappear as the temperature is increased above $`T_\mathrm{c}`$. Figure 8(a) shows the temperature dependence of the number of particles above $`T_\mathrm{c}`$ (for $`\nu =0`$). The number of stable Bosons ($`N_\mathrm{B}^{\gamma =0}`$ and $`N_\mathrm{C}^{\gamma =0}`$), as well as the finite lifetime Bosons ($`N_\mathrm{B}^{\gamma >0}`$) and the scattering contribution ($`N_\mathrm{C}^{\mathrm{sc}}`$), are found to all decrease at higher temperatures. This result is also obtained when $`\nu <0`$. Thus only the free Fermi atoms contribute at higher temperatures well above $`T_\mathrm{c}`$.

The reason that the Fermi atom contribution dominates at high temperatures has a very simple explanation in the difference of the chemical potentials between a Fermi atom ($`\mu _\mathrm{F}=\mu `$) and a $`b`$-Boson ($`\mu _\mathrm{B}=2\mu `$). As the temperature is increased, the chemical potential $`\mu `$ must decrease so as to conserve the total number of atoms \[$`N`$=(Fermi atoms)+2$`\times `$(Bose molecules)\], as shown in Fig. 8(b). Thus for $`TT_\mathrm{c}`$, the Boson excitation energy $`E_\text{q}^0+2\nu 2\mu `$ always increases faster than the excitation energy of the Fermions $`\epsilon _\text{p}\mu `$, simply due to the factor two in front of $`\mu `$ in the Bose excitation energy. As a result, the occupation of the Bose band always approaches zero for $`TT_\mathrm{c}`$, even though the Bosons can be the dominant excitation at $`T_\mathrm{c}`$ in the BEC regime.

However, Fig. 8 still shows that there is a significant region above $`T_\mathrm{c}`$ where the number of long-lived Bosons is substantial. It is thus useful to define a “crossover temperature” $`T^{}`$ which separates the Fermi gas and the Bose gas regimes, even though the increase in the number of stable Bosons occurs continuously. We define $`T^{}`$ as the temperature at which $`N_\mathrm{B}^{\gamma =0}=0.05`$, or when 10$`\%`$ of the Fermi atoms combine to form stable $`b`$-molecules (Fig. 9). Figure 9 shows that $`T^{}`$ rapidly increases when $`\nu `$ is small or negative, for which case stable Bosons appear at $`T_\mathrm{c}`$ ($`\nu /\epsilon _\mathrm{F}<0.18`$, for the parameters used). Even for $`\nu 0`$, $`T^{}\epsilon _\mathrm{F}`$ is considerably larger than $`T_\mathrm{c}`$. Since this “Bose gas regime” lying between $`T_\mathrm{c}`$ and $`T^{}`$ is fairly wide, one might expect to observe the formation of the stable molecules somewhat above $`T_\mathrm{c}`$ if we can decrease the threshold energy $`2\nu `$ enough, namely to values where the superfluid transition at $`T_\mathrm{c}`$ would be of the BEC type (see Fig. 4).

Figure 10 compares the temperature dependence of the number of particles and the chemical potential (lines in the figure) for $`\nu =\epsilon _\mathrm{F}`$ with the case of a free ($`U=g_\mathrm{r}=0`$) Fermion-Boson mixture (filled circles and squares). In the non-interacting Fermion-Boson mixture, we only impose the constraint that the total number of atoms is conserved \[$`N`$=(Fermi atoms)+2$`\times `$(Bose molecules)\]. Figure 10 shows that deep in the BEC regime ($`\nu =\epsilon _\mathrm{F}`$), this constraint alone determines the temperature dependence of the number of particles and the chemical potential. As far as these quantities are concerned, the effect of including the interaction $`g_\mathrm{r}`$ and $`U`$ is clearly negligible.

## VII The case of strong Feshbach coupling

In all of the numerical results presented so far, we have considered the case of a weak Feshbach coupling described by $`g_\mathrm{r}=0.6\epsilon _\mathrm{F}`$. In this regard, recent calculations have argued that strong coupling ($`g_\mathrm{r}\epsilon _\mathrm{F}`$) was a more realistic description of Feshbach resonances observed in ultracold gases. In this section, we investigate what arises in this strong Feshbach coupling regime, although the present treatment of fluctuations implicitly assumes a weak Feshbach coupling. We present some numerical results for $`T_\mathrm{c}`$ in the limit of a large Feshbach coupling (broad resonance) for a uniform gas. In this regard, we note that our analysis in the present paper and in Ref. has shown that the BCS-BEC crossover in a trap is very similar to that in a uniform gas, at least for the case of weak Feshbach coupling.

Figure 11 show the BCS-BEC crossover in the case of $`g_\mathrm{r}=20\epsilon _\mathrm{F}\epsilon _\mathrm{F}`$. Figures 11(a) and 11(b) show that the crossover behavior already occurs around $`\nu =150\epsilon _\mathrm{F}`$, which is much higher than $`\epsilon _\mathrm{F}`$. Results similar to this are presented in Ref.. Since the bottom of the Bose excitation spectrum $`2\nu `$ is then still much higher than the Fermion band, the Feshbach related $`b`$-Bosons described by $`N_\mathrm{B}=N_\mathrm{B}^{\gamma =0}+N_\mathrm{B}^{\gamma >0}`$ are almost absent in this crossover regime, as shown in Fig. 11(c). The high-energy $`b`$-Bosons only contribute to virtual scattering processes involving Fermions, which mediate the pairing interaction between Fermi atoms. Thus, when the Feshbach coupling $`g_\mathrm{r}`$ is very strong, the character of the BCS-BEC crossover becomes identical to that in the case of strong-coupling superconductivity as discussed by Nozières and Schmitt-Rink. The fact that the crossover to the BEC region occurs at such large values of $`2\nu `$ is easy to understand. Since the chemical potential $`\mu `$ is always of order of $`\epsilon _\mathrm{F}`$ or less, at large values of $`2\nu `$ and $`g_\mathrm{r}`$, we have $`U_{\mathrm{eff}}=U+g_\mathrm{r}^2/(2\nu 2\mu )g_\mathrm{r}^2/2\nu `$. Since the strong-coupling effect should be most pronounced when $`U_{\mathrm{eff}}`$ becomes comparable to $`\epsilon _\mathrm{F}`$, the BCS-BEC crossover will occur at $`U_{\mathrm{eff}}\epsilon _\mathrm{F}`$, which gives $`2\nu g_\mathrm{r}^2/\epsilon _\mathrm{F}`$. For $`g_\mathrm{r}=20\epsilon _\mathrm{F}`$, this predicts that the crossover will occur at $`\nu 200\epsilon _\mathrm{F}`$, consistent with the numerical results in Fig. 11(a).

That the case of a very broad Feshbach resonance ($`g_\mathrm{r}\epsilon _\mathrm{F}`$) is equivalent to the case studied in Ref. is also shown clearly by the fact that stable preformed Cooper-pairs are the dominant excitation \[see Fig. 11(c)\] in the crossover region at $`\nu 150\epsilon _\mathrm{F}`$. The Feshbach $`b`$-molecule component only becomes dominant when $`\nu <0`$. This result is in contrast to our results for a weak Feshbach coupling $`g_\mathrm{r}=0.6\epsilon _\mathrm{F}`$ as shown in Fig. 4. In that case, the $`b`$-molecule becomes the dominant component at the BCS-BEC crossover.

Figure 12 shows the total number of composite Bosons $`N_\mathrm{M}N_\mathrm{B}^{\gamma =0}+N_\mathrm{B}^{\gamma >0}+N_\mathrm{C}^{\gamma =0}+N_\mathrm{C}^{\mathrm{sc}}`$ as a function of $`\nu `$, for $`g_\mathrm{r}=20\epsilon _\mathrm{F}`$. Although the individual Cooper-pair components $`N_\mathrm{C}^{\gamma =0}`$ and $`N_\mathrm{C}^{\mathrm{sc}}`$ show \[see Fig. 11(c)\] singular behavior in the BCS-BEC crossover regime (at $`\nu 150\epsilon _\mathrm{F}`$), $`N_\mathrm{M}`$ itself increases smoothly as the threshold energy is decreased.

Since the BCS-BEC crossover around $`\nu =150\epsilon _\mathrm{F}`$ shown in Fig. 11 is intrinsically the same phenomenon as discussed by Nozières and Schmitt-Rink, the origin of the small peak in $`T_\mathrm{c}`$ shown in Fig. 11(a) is the analogue of the one formed also in the case of strong-coupling superconductivity. However, Haussmann showed that the slight peak in $`T_\mathrm{c}`$ obtained at the BCS-BEC crossover in strong-coupling superconductivity disappears when one includes higher order scattering processes beyond the simple $`t`$-matrix approximation, such as used in the present paper. Thus, we regard the small peak $`T_\mathrm{c}`$ in Fig. 11(a) as probably an artifact of our simple $`t`$-matrix approximation. As in the weak Feshbach coupling case shown in Fig. 4, we expect that the maximum transition temperature $`T_\mathrm{c}`$ will again be given by the BEC expressions in eqs. (82) and (83).

## VIII Summary

To summarize, we have extended our recent work and investigated the effect of particle-particle fluctuations in conjunction with an atomic Feshbach resonance on the superfluid phase transition in a gas of Fermi atoms in a harmonic trap, using a simple LDA approach. The BCS-BEC crossover predicts a maximum transition temperature of $`T_\mathrm{c}=0.518T_\mathrm{F}`$, somewhat larger than $`T_\mathrm{c}=0.218T_\mathrm{F}`$ in a uniform gas. In the BEC regime, the excitation spectrum is dominated by the stable long-lived Feshbach-related molecules in the case of a small Feshbach coupling parameter $`g_\mathrm{r}`$. This is in contrast with results for the BCS-BEC crossover with no Feshbach resonance, where the preformed Cooper-pairs are the dominant bound states. For a very broad Feshbach resonance ($`g_\mathrm{r}\epsilon _\mathrm{F}`$) as studied in Refs., we find that the crossover region is also dominated by preformed Cooper-pairs.

The BCS-BEC crossover is shown to strongly affect the atomic density profile in the trap, which can be measured. In the BCS regime, where Fermi atoms are dominant excitations, the spatial distribution of the atoms is spread out from the trap center due to the Pauli exclusion principle. The atoms increasingly cluster at the center of the trap as the threshold energy of $`b`$-Bose excitation spectrum is lowered and the the number of stable composite Bosons increases.

In the BEC regime, we show that stable molecules can exist above $`T_\mathrm{c}`$. As the temperature increases, however, the number of these Bosons decreases. We introduce a crossover temperature $`T^{}`$ defined as the temperature above $`T_\mathrm{c}`$ at which 10$`\%`$ of the Fermi atoms have formed long-lived $`b`$-molecules. Plotting $`T^{}`$ as a function of the threshold $`2\nu `$ gives a sort of “phase diagram” in the normal phase which may be useful in experimental searches for the BCS-BEC crossover.

The experimental observation of Cooper-pairs in trapped Fermi gases will be a very exciting milestone in physics. What makes this topic even more interesting is the possibility that we can observe the BCS-BEC crossover in such systems. As we have shown in this paper, this crossover involves very interesting physics and, moreover, it leads to characteristic changes in the properties of the trapped gas. The changes in the atomic density profile may provide a crucial experimental signature for the existence of Cooper pairing and superfluidity.

We will discuss the BCS-BEC crossover in the superfluid region below $`T_\mathrm{c}`$ in a future paper. In particular, we will discuss how the collective modes vary as we pass through the crossover region.

###### Acknowledgements.

The authors would like to acknowledge useful discussion with M. Holland and J. Milstein concerning Refs.. Y. O would like to thank Professor S. Takada for valuable discussions, as well as financial support from the Japanese Government while a visiting professor at the University of Toronto. A.G. acknowledges support from NSERC of Canada.

## A Calculation of fluctuation contribution to the number of Fermi atoms

The density of atoms associated with fluctuations given by $`\delta N_{\mathrm{FL}}(𝐫)`$ in eq. (29) can be split into two contributions \[using eq. (25)\],

$$\delta N_{\mathrm{FL}}(𝐫)=\frac{1}{\beta }\underset{\text{q},\nu _n}{}e^{i\nu _n\delta }\frac{U}{1U\mathrm{\Pi }}\frac{\mathrm{\Pi }}{\mu }\frac{1}{\beta }\underset{\text{q},\nu _n}{}e^{i\nu _n\delta }\frac{}{\mu }\mathrm{ln}\left[1+g_\mathrm{r}^2D_0\stackrel{~}{\mathrm{\Pi }}\right].$$
(A1)
Since the second term involving $`\frac{}{\mu }\mathrm{ln}[1+g_\mathrm{r}^2D_0\stackrel{~}{\mathrm{\Pi }}]`$ goes as $`\nu _n^2`$ in the large $`\nu _n`$ limit, the convergence factor $`e^{i\nu _n\delta }`$ has no effect. However, this convergence factor cannot be dropped in the first term in eq. (A1), because the quantity $`\frac{U}{1U\mathrm{\Pi }}\frac{\mathrm{\Pi }}{\mu }`$ only goes as $`\nu _n^1`$ for large $`\nu _n`$. In numerical calculations, it is difficult to sum up the Matsubara frequencies directly taking into account the convergence factor involving a infinitesimally small number $`\delta `$. One method to avoid this difficulty is to employ a spectral representation as carried out by Nozières and Schmitt-Rink. However, this method requires high numerical accuracy in executing the resulting energy integration. We used a different approach.

To calculate the first term in eq. (A1) \[$`\delta N_{\mathrm{FL}}^{(1)}(𝐫)`$\] without employing a spectral representation, we compare $`\delta N_{\mathrm{FL}}^{(1)}(𝐫)`$ to the expression without the convergence factor,

$$\delta N_{\mathrm{FL}}^{(2)}(𝐫)\frac{1}{\beta }\underset{\text{q},\nu _n}{}\frac{U}{1U\mathrm{\Pi }}\frac{\mathrm{\Pi }}{\mu }.$$
(A2)
This quantity is easy to calculate numerically. When we transform the summation over the discrete Matsubara frequencies into an integral in the complex plane in the usual way, we obtain

$`\{\begin{array}{c}\delta N_{\mathrm{FL}}^{(1)}(𝐫)={\displaystyle \frac{1}{2\pi i}}{\displaystyle \underset{\text{q}}{}}{\displaystyle _{C_0}}𝑑ze^{\delta z}n_\mathrm{B}(z){\displaystyle \frac{U}{1U\mathrm{\Pi }(i\nu _nz)}}{\displaystyle \frac{\mathrm{\Pi }(i\nu _nz)}{\mu }},\hfill \\ \delta N_{\mathrm{FL}}^{(2)}(𝐫)={\displaystyle \frac{1}{2\pi i}}{\displaystyle \underset{\text{q}}{}}{\displaystyle _{C_0}}𝑑zn_\mathrm{B}(z){\displaystyle \frac{U}{1U\mathrm{\Pi }(i\nu _nz)}}{\displaystyle \frac{\mathrm{\Pi }(i\nu _nz)}{\mu }},\hfill \end{array}`$ (A5)
where the path $`C_0`$ is shown in Fig. 13. we are allowed to add the circular contours $`C_1`$ and $`C_2`$ for the integration in $`\delta N_{\mathrm{FL}}^{(1)}(𝐫)`$ without changing the result, because the contributions from $`C_1`$ and $`C_2`$ vanish in the large radius limit ($`R\mathrm{}`$) due to, respectively, the factors $`n_\mathrm{B}(z)`$ and $`e^{\delta z}`$. The integration in $`\delta N_{\mathrm{FL}}^{(1)}(𝐫)`$ is thus reduced to summing up the poles inside the closed contours $`C_0+C_1`$ and $`C_0+C_2`$ shown in Fig. 13.

We next show how $`\delta N_{\mathrm{FL}}^{(1)}(𝐫)`$ is related to $`\delta N_{\mathrm{FL}}^{(2)}(𝐫)`$ defined in eq. (A2). First, we note that when dealing with $`\delta N_{\mathrm{FL}}^{(2)}(𝐫)`$, the integration along the path $`C_2`$ gives a finite contribution when we add the paths $`C_1`$ and $`C_2`$ due to the absence of the convergence factor. However, if we write the integration in $`\delta N_{\mathrm{FL}}^{(2)}(𝐫)`$ as

$$_{C_0}𝑑z=_{C_0+C_1+C_2}𝑑z_{C2}𝑑z,$$
(A6)
we find that the first term can be calculated by summing up the contribution from the poles inside the closed paths $`C_0+C_1`$ and $`C_0+C_2`$, and the result is identical to $`\delta N_{\mathrm{FL}}^{(1)}(𝐫)`$. The last term in eq. (A6) \[$`\delta N_{\mathrm{FL}}^{\mathrm{corr}}(𝐫)`$\] can be evaluated to give, in the limit $`R\mathrm{}`$,

$`\delta N_{\mathrm{FL}}^{\mathrm{corr}}(𝐫)`$ $`=`$ $`{\displaystyle \frac{U}{2\pi }}{\displaystyle \underset{\text{p},\text{q}}{}}e^{\epsilon _\text{p}^2/\omega _c^2}{\displaystyle _{C_2}}𝑑\theta \left[{\displaystyle \frac{f(\epsilon _{\text{p}+\text{q}/2}\mu (𝐫))}{\mu }}+{\displaystyle \frac{f(\epsilon _{\text{p}\text{q}/2}\mu (𝐫))}{\mu }}\right]`$ (A7)
$`=`$ $`U{\displaystyle \underset{\text{p},\text{q}}{}}e^{\epsilon _\text{p}^2/\omega _c^2}{\displaystyle \frac{f(\epsilon _{\text{p}+\text{q}/2}\mu (𝐫))}{\mu }}.`$ (A8)
Here, $`e^{\epsilon _\text{p}^2/\omega _c^2}`$ is the Gaussian cutoff introduced in Section IV, and we have used eq. (15) and changed the integration variable using $`z=Re^{i\theta }`$. Thus we have shown that $`\delta N_{\mathrm{FL}}^{(1)}(𝐫)`$ involving the convergence factor is related to $`N_{\mathrm{FL}}^{(2)}(𝐫)`$ without the convergence factor by

$$\delta N_{\mathrm{FL}}^{(1)}(𝐫)=\delta N_{\mathrm{FL}}^{(2)}(𝐫)\delta N_{\mathrm{FL}}^{\mathrm{corr}}(𝐫).$$
(A9)
Using this expression, we can calculate $`\delta N_{\mathrm{FL}}^{(1)}(𝐫)`$ numerically by evaluating the simpler expression $`\delta N_{\mathrm{FL}}^{(2)}(𝐫)`$.

The essence of this prescription may be easily understood from the following example: Consider the two quantities:

$`\{\begin{array}{c}N_\mathrm{B}^{(1)}(𝐪)={\displaystyle \underset{\nu _n}{}}e^{i\nu _n\delta }{\displaystyle \frac{1}{i\nu _nE_\text{q}}},\hfill \\ N_\mathrm{B}^{(2)}(𝐪)={\displaystyle \underset{\nu _n}{}}{\displaystyle \frac{1}{i\nu _nE_\text{q}}}.\hfill \end{array}`$ (A12)
In this case, the frequency sum in $`N_\mathrm{B}^{(1)}(𝐪)`$ can be easily done analytically and reduced to the Bose distribution function $`n_\mathrm{B}(E_\text{q})`$. The sum in $`N_\mathrm{B}^{(2)}`$ gives

$$N_\mathrm{B}^{(2)}(𝐪)=\underset{\nu _n}{}\frac{E_\text{q}}{\nu _n^2+E_\text{q}^2}=\frac{1}{2}\mathrm{coth}\frac{\beta E_\text{q}}{2}.$$
(A13)
The contribution from $`C_2`$ in the last term in eq. (A6) is given by

$$\delta N_\mathrm{B}^{\mathrm{corr}}\frac{1}{2\pi i}_{C_2:R\mathrm{}}𝑑z\frac{1}{zE_\text{q}}=\frac{1}{2}.$$
(A14)
Thus we can easily verify that $`N_\mathrm{B}^{(1)}=N_\mathrm{B}^{(2)}\delta N_\mathrm{B}^{\mathrm{corr}}`$, as in eq. (A9).