format: text

text:
# Superfluidity and collective modes in a uniform gas of Fermi atoms with a Feshbach resonance

## I Introduction

One of the most challenging topics in current physics is the realization of superfluidity in a trapped atomic gas composed of two Fermi hyperfine states. A considerable experimental effort has already been made to cool down Fermi atom gases, such as 6Li and 40K. The temperature can now be lowered to $`T0.2T_\mathrm{F}`$, where the Fermi gas should be highly degenerate and the observation of superfluid behavior seems imminent.

As a promising mechanism of superfluidity with a high phase transition temperature $`T_\mathrm{c}`$, making use of an atomic Feshbach resonance has attracted much attention. The Feshbach resonance describes quasi-molecular bosons, which can mediate a pairing interaction between Fermi atoms. This pairing interaction is tunable by the threshold energy $`2\nu `$ of the Feshbach resonance, and can become strong as $`2\nu `$ is decreased relative to twice the Fermi energy of the atoms. Using this strong paring interaction, one can hope to achieve a high value of $`T_\mathrm{c}`$. Experimentally, the threshold energy $`2\nu `$ can be controlled by a weak applied magnetic field. Very recently, this tunable interaction was observed in a Fermi gas of 40K.

In our recent work, we pointed out the importance of fluctuations in the Cooper-channel in considering a high-$`T_\mathrm{c}`$ superfluidity originating from the strong pairing interaction associated with the Feshbach resonance. We extended the strong-coupling theory developed by Nozières and Schmitt-Rink, to include the effects of a Feshbach resonance and the associated quasi-molecular bosons. We showed that these particle-particle fluctuations strongly suppress $`T_\mathrm{c}`$ from the value expected within the simple mean-field BCS theory. In addition, the character of the phase transition was shown to continuously change from the BCS-type to a BEC of composite bosons (consisting of preformed Cooper-pairs and long-lived Feshbach molecules) as the threshold energy $`2\nu `$ is lowered. Our strong-coupling theory thus gave an upper limit of $`T_\mathrm{c}=0.518T_\mathrm{F}`$ for a Fermi gas in an isotropic harmonic trap potential and $`T_\mathrm{c}=0.218T_\mathrm{F}`$ for a uniform gas. These values are simply the BEC transition temperatures, expressed in terms of the Fermi temperature $`T_\mathrm{F}`$ of one of the Fermi components.

In this paper, we investigate the BCS-BEC crossover in the superfluid state, extending our previous work to the superfluid region below $`T_\mathrm{c}`$. Going past the previous BCS mean-field approximation, we include strong fluctuations around the BCS mean-field solution. We clarify how the order parameter described by the Cooper-pair amplitude $`\mathrm{\Delta }=U_\text{p}c_\text{p}c_\text{p}`$ in the weak-coupling BCS theory changes to the BEC order parameter related to the number of condensed bosons $`\varphi _\mathrm{m}=b_{\text{q}=0}`$ in the BCS-BEC crossover. Here $`c_{\text{p}\sigma }`$ is the annihilation operator of a Fermi atom in one of two hyperfine states ($`\sigma =,`$) and $`b_\text{q}`$ is the annihilation operator of the boson molecule associated with the Feshbach resonance.

In the field of trapped ultra-cold Fermi gases, a crucial issue is to determine a clear unambiguous signature for superfluidity. Another important problem is how to experimentally track the system in the BCS-BEC crossover region. In this regard, the Goldstone collective mode is a very useful quantity because it is deeply related to the spontaneous breakdown of the gauge symmetry associated with the superfluid phase transition. The Goldstone mode is known as the Anderson-Bogoliubov mode in the BCS state, while it is the Bogoliubov phonon in the BEC phase. In this paper, we discuss how these collective modes change from one to the other as we go through the BCS-BEC crossover. We show that the velocity of the Goldstone phonon $`v_\varphi `$ strongly depends on $`2\nu `$, and thus it offers a way of observing the BCS-BEC crossover phenomenon by tuning the threshold energy $`2\nu `$ in a cigar shaped trap (where the gas is fairly uniform in the axial direction). We also show that the Goldstone mode appears as a resonance in the spectrum of the density-density correlation function.

The present paper only considers the superfluid phase in a uniform two-component Fermi gas with an attractive interaction. In Ref. , we discussed the same model at and above the superfluid transition temperature for a trapped gas, using the local density approximation (LDA). The extension of the present work to an inhomogeneous trapped gas will be considered in the future. As well known, the low energy collective modes in a trapped atomic gas (Fermi or Bose) are strongly altered by the trap potential. However, the equivalent of the Goldstone phonon modes we discuss in this paper also arise in trapped two component Fermi gases. These low energy collective modes have been extensively discussed in the recent literature for the weak-coupling BCS limit. In this paper, we discuss the physics of the Goldstone phonon mode as a function of the threshold energy $`2\nu `$. A similar analysis remains to be done for the collective modes of a trapped Fermi gas in the BCS-BEC crossover region.

This paper is organized as follows. In Section II, we present our coupled fermion-boson model. We explain how to include the strong-coupling (fluctuation) effect originating from the Feshbach resonance in Section III. In Section IV, we consider the strong-coupling theory presented in Section III in terms of Gaussian fluctuations. In Section V, we consider the Goldstone mode. We first derive correlation functions describing Cooper-pair fluctuations, as well as a renormalized boson Green’s function for quasi-molecules associated with the Feshbach resonance, under the Hartree-Fock random phase approximation (HF-RPA). The Goldstone mode is then obtained from their poles. In Section VI, we discuss the BCS-BEC crossover below $`T_\mathrm{c}`$ based on our numerical results. In this section, we also discuss the BCS-BEC crossover behavior of the order parameter and the Goldstone phonon mode, as a function of $`2\nu `$. In Section VII, we discuss the coupling of the Goldstone mode with the density fluctuations in a gas of Fermi atoms.

## II Coupled fermion-boson model

We consider a gas of Fermi atoms composed of two atomic hyperfine states, coupled to molecular two particle state. We describe the two hyperfine states using a pseudo-spin variable $`\sigma `$ ($`=,`$). The coupled fermion-boson model Hamiltonian is given by

$``$ $`=`$ $`{\displaystyle \underset{\text{p}\sigma }{}}\epsilon _\text{p}c_{\text{p}\sigma }^{}c_{\text{p}\sigma }U{\displaystyle \underset{\text{p},\text{p}^{},\text{q}}{}}c_{\text{p}+\frac{\text{q}}{2}}^{}c_{\text{p}+\frac{\text{q}}{2}}^{}c_{\text{p}^{}+\frac{\text{q}}{2}}c_{\text{p}^{}+\frac{\text{q}}{2}}`$ (1)
$`+`$ $`{\displaystyle \underset{\text{q}}{}}[\epsilon _{B\text{q}}+2\nu ]b_\text{q}^{}b_\text{q}+g_\mathrm{r}{\displaystyle \underset{\text{p},\text{q}}{}}[b_\text{q}^{}c_{\text{p}+\frac{\text{q}}{2}}c_{\text{p}+\frac{\text{q}}{2}}+\mathrm{h}.\mathrm{c}.].`$ (2)
Here a Fermi atom and a quasi-molecular boson associated with the Feshbach resonance are, respectively, described by the destruction operators $`c_{\text{p}\sigma }`$ and $`b_\text{q}`$. The kinetic energy of a Fermi atom is $`\epsilon _\text{p}p^2/2m`$, and $`\epsilon _{B\text{q}}+2\nu q^2/2M+2\nu `$ is the excitation spectrum of the $`b`$-molecular bosons. Here $`2\nu `$ represents the lowest excitation energy of $`b`$-bosons, also referred to as the threshold energy of the Feshbach resonance. The last term in Eq. (2) describes the Feshbach resonance with a coupling constant $`g_\mathrm{r}`$, which describes how a $`b`$-molecule can dissociate into two Fermi atoms and how two Fermi atoms can form one $`b`$-boson. The Hamiltonian in Eq. (2) also includes an attractive fermion-fermion interaction $`U`$ ($`<0`$) arising from non-resonant processes.

Since one $`b`$-Bose molecule consists of two Fermi atoms, the boson mass $`M=2m`$ and the conservation of the total number of particles $`N`$ imposes the relation

$$N=\underset{\text{p}\sigma }{}c_{\text{p}\sigma }^{}c_{\text{p}\sigma }+2\underset{\text{q}}{}b_\text{q}^{}b_\text{q}.$$
(3)
We incorporate this latter constraint into the model Hamiltonian in Eq. (2) using a chemical potential, $`H\mu N`$. The resulting grand-canonical Hamiltonian $``$ has the same form as Eq. (2), except that the kinetic energies of Fermi atoms and $`b`$-bosons are replaced as $`\epsilon _\text{p}\xi _\text{p}\epsilon _\text{p}\mu `$ and $`\epsilon _{B\text{q}}+2\nu \xi _{B\text{q}}\epsilon _{B\text{q}}+2\nu 2\mu `$, respectively. In the latter replacement, the factor $`2`$ in $`2\mu `$ reflects that one $`b`$-boson consists of two Fermi atoms.

In this paper, we investigate strong-coupling effects in the superfluid phase, as well as the Goldstone mode associated with superfluidity in the BCS-BEC crossover region. As a start, we consider a simple uniform Fermi gas and leave the effect of a trapping potential to future work. In this regard, we have shown in Refs.. that while a trap potential enhances the transition temperature $`T_\mathrm{c}`$ in the BEC regime, the qualitative behavior of $`T_\mathrm{c}`$ in the BCS-BEC crossover is not very different from a uniform Fermi gas. Within weak-coupling BCS theory, several papers have discussed collective excitations in a trapped Fermi gas with attractive interactions. (See, for example, Refs. .)

When the Feshbach coupling term is absent in Eq. (2), the fermions and $`b`$-bosons are decoupled from each other. In this limit, a BCS superfluid phase transition of Fermi atoms and BEC transition of $`b`$-bosons can occur, at different temperatures. These two superfluid phases are, respectively, described by independent order parameters

$`\{\begin{array}{c}\mathrm{\Delta }U_\text{p}c_\text{p}c_\text{p},\hfill \\ \varphi _\mathrm{m}b_{\text{q}=0}.\hfill \end{array}`$ (6)
On the other hand, when the Feshbach resonance term is present ($`g_\mathrm{r}0`$), we find the following identity in the equilibrium state:

$`0=i{\displaystyle \frac{d\varphi _\mathrm{m}}{dt}}=i{\displaystyle \frac{db_0}{dt}}=[b_0,]=(2\nu 2\mu )\varphi _\mathrm{m}+{\displaystyle \frac{g_\mathrm{r}}{U}}\mathrm{\Delta }.`$ (7)
Eq. (7) gives

$$\varphi _\mathrm{m}=\frac{g_\mathrm{r}}{2\nu 2\mu }\frac{\mathrm{\Delta }}{U}.$$
(8)
This last result shows that the BEC order parameter $`\varphi _\mathrm{m}`$ and the Cooper-pair order parameter $`\mathrm{\Delta }`$ are no longer independent, due to the hybridization induced by the Feshbach resonance $`g_\mathrm{r}`$. Both $`\mathrm{\Delta }`$ and $`\varphi _\mathrm{m}`$ are finite in the superfluid phase, and there is a unique superfluid phase transition in this coupled fermion-boson model.

For later convenience, we define the following composite order parameter

$$\stackrel{~}{\mathrm{\Delta }}\mathrm{\Delta }g_\mathrm{r}\varphi _\mathrm{m}.$$
(9)
We will find that $`\stackrel{~}{\mathrm{\Delta }}`$ determines the excitation energy gap in the spectrum of fermion quasi-particles below $`T_\mathrm{c}`$ in the BCS-BEC crossover regime.

## III Strong-coupling effect on superfluidity

#### 1 Review on the strong-coupling theory for $`T_\mathrm{c}`$

In this section, we review the strong-coupling theory for $`T_\mathrm{c}`$ discussed in our previous papers. This formulation is extended to the region below $`T_\mathrm{c}`$ in the next subsection.

In previous work, we extended the strong-coupling theory developed by Nozières and Schmitt-Rink to the coupled fermion-boson model in Eq. (2). The equation for $`T_\mathrm{c}`$ was obtained by using the Thouless criterion, which states that the superfluid phase transition occurs when the particle-particle vertex function $`\mathrm{\Gamma }(𝐪,\omega )`$ describing the Cooper-channel develops a pole at $`𝐪=\omega =0`$. Within the $`t`$-matrix approximation in terms of $`U`$ and $`g_\mathrm{r}`$ described diagrammatically in Fig. 1(a), this equation for $`T_\mathrm{c}`$ is given by

$$1=U_{\mathrm{eff}}\underset{\text{p}}{}\frac{\mathrm{tanh}2\xi _\text{p}/2T_\mathrm{c}}{2\xi _\text{p}},$$
(10)
where

$$U_{\mathrm{eff}}U+g_\mathrm{r}^2\frac{1}{2\nu 2\mu }$$
(11)
is an effective pairing interaction. The last term in Eq. (11) describes the interaction mediated by $`b`$-bosons associated with the Feshbach resonance. Eq. (10) is formally identical to the equation for $`T_\mathrm{c}`$ in ordinary weak-coupling BCS theory. However, the chemical potential $`\mu `$ in the kinetic energy $`\xi _\text{p}=\epsilon _\text{p}\mu `$ of the Fermi atoms can deviate strongly from $`\epsilon _\mathrm{F}`$ as one approaches the BEC regime (where $`\epsilon _\mathrm{F}`$ is the bare Fermi energy of one spin component). This contrasts with simple BCS theory, where one finds that $`\mu \epsilon _\mathrm{F}`$.

The chemical potential $`\mu `$ is determined by the equation for the total number of Fermi atoms, using the identity $`N=\mathrm{\Omega }/\mu `$. We include the effect of fluctuations in the Cooper-channel \[the first line in Fig. 1(b)\] as well as the Feshbach resonance \[the second line in Fig. 1(b)\] in the thermodynamic potential $`\mathrm{\Omega }`$. The resulting equation relating $`\mu `$ and $`N`$ is

$$N=N_\mathrm{F}^0+2N_\mathrm{B}^0\frac{1}{\beta }\underset{\text{q},\nu _n}{}e^{i\delta \nu _n}\frac{}{\mu }\mathrm{ln}\left[1[Ug_\mathrm{r}^2D_0(𝐪,i\nu _n)]\mathrm{\Pi }(𝐪,i\nu _n)\right].$$
(12)
Here $`\beta =1/T`$ is the inverse of temperature, $`N_\mathrm{F}^02_\text{p}f(\xi _\text{p})`$ and $`N_\mathrm{B}^0_\text{q}n_\mathrm{B}(\xi _{B\text{q}})`$, where $`f(\epsilon )`$ and $`n_\mathrm{B}(\epsilon )`$ represent the Fermi and Bose distribution functions, respectively. The last term in Eq. (12) describes the fluctuation contribution to $`N`$. The $`b`$-boson Green’s function is given by

$$D_0(𝐫,i\nu _n)\frac{1}{i\nu _n\xi _{B\text{q}}},$$
(13)
where the Bose Matsubara frequency is $`\nu _n=2n\pi T(n=0,\pm 1,\pm 2,\mathrm{})`$. $`\mathrm{\Pi }(𝐪,i\nu _n)`$ is the correlation function of the Cooper-pair field operator $`\widehat{F}_\text{q}_\text{p}c_{\text{p}+\frac{\text{q}}{2}}c_{\text{p}+\frac{\text{q}}{2}}`$, given by

$`\mathrm{\Pi }(𝐪,i\nu _n)`$ $``$ $`{\displaystyle _0^\beta }𝑑\tau e^{i\nu _n\tau }T_\tau \{\widehat{F}_\text{q}(\tau )\widehat{F}_\text{q}^{}(0)\}`$ (14)
$`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\text{p},\omega _m}{}}G_0(𝐩+{\displaystyle \frac{𝐪}{2}},i\omega _m+i\nu _n)G_0(𝐩+{\displaystyle \frac{𝐪}{2}},i\omega _m)`$ (15)
$`=`$ $`{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{1f(\xi _{\text{p}+\frac{\text{q}}{2}})f(\xi _{\text{p}\frac{\text{q}}{2}})}{\xi _{\text{p}+\frac{\text{q}}{2}}+\xi _{\text{p}\frac{\text{q}}{2}}i\nu _n}}.`$ (16)
Physically, $`\mathrm{\Pi }(𝐪,i\nu _n)`$ describes fluctuations of Cooper-pairs in the normal phase. $`G_0(𝐪,i\omega _m)`$ is a fermion thermal Green’s function defined by

$$G_0(𝐪,i\omega _m)=\frac{1}{i\omega _m\xi _\text{p}},$$
(17)
where the Fermi Matsubara frequency is $`\omega _m=(2m+1)\pi T(m=0,\pm 1,\pm 2,\mathrm{})`$.

The coupled equations (10) and (12) determine $`T_\mathrm{c}`$ and $`\mu `$ self-consistently. In calculating these equations, a cutoff is necessary to make the momentum summation converge. In Refs., we simply introduced a Gaussian cutoff $`e^{\epsilon _\text{p}^2/\omega _c^2}`$, which is also used in this paper.

The self-consistent solution ($`T_\mathrm{c}`$ and $`\mu `$) of these coupled equations is summarized in Fig.2. When $`\nu \epsilon _\mathrm{F}`$, since the chemical potential is at most $`\mu <\epsilon _\mathrm{F}`$, the Feshbach-induced contribution to the pairing interaction $`g_\mathrm{r}^2/(2\nu 2\mu )`$ in Eq. (11) is small. In this regime, Fig. 2(a) shows that the superfluid phase transition is well described by the weak-coupling BCS theory with a (weak) pairing interaction $`U`$. In addition, we see that $`\mu \epsilon _\mathrm{F}`$, as shown in Fig. 2(b). However, the chemical potential $`\mu `$ gradually deviates from $`\epsilon _\mathrm{F}`$ as the threshold energy $`2\nu `$ is lowered towards $`2\epsilon _\mathrm{F}`$ and below. In particular, one finds $`\mu `$ approaches $`\nu `$ when $`\nu <0`$. In this regime, the Feshbach-induced pairing interaction $`g_\mathrm{r}^2/(2\nu 2\mu )`$ in Eq. (11) is large, and $`T_\mathrm{c}`$ deviates significantly from the prediction of weak-coupling BCS theory, as shown in Fig. 2(a). Since $`2\nu `$ is the lowest excitation energy of $`b`$-bosons and their chemical potential is $`2\mu `$, the situation $`2\nu =2\mu `$, realized in the limit of large negative values of $`\nu /\epsilon _\mathrm{F}`$, is equivalent to the condition of BEC in a non-interacting Bose gas. Indeed, Fig. 2(a) shows that $`T_\mathrm{c}`$ corresponds precisely to the transition temperature of a free Bose gas of $`N/2`$ atoms when $`\nu /\epsilon _\mathrm{F}<1`$. These results show that the BCS-BEC crossover occurs in the region around zero threshold ($`2\nu =0`$), at least for small values of $`g_\mathrm{r}`$. This crossover phenomenon can be simply controlled by the threshold energy $`2\nu `$ of the Feshbach resonance.

#### 2 Strong-coupling theory below $`T_\mathrm{c}`$

In order to formulate the analogous strong-coupling theory below $`T_\mathrm{c}`$, we separate out the fluctuations of Cooper-pairs and condensed $`b`$-bosons around their mean-field values denoted by $`\mathrm{\Delta }`$ and $`\varphi _\mathrm{m}`$, respectively. For this purpose, we write $`\widehat{F}_\text{q}=\widehat{F}_{\text{q}=0}\delta _{\text{q},0}+\delta \widehat{F}_\text{q}`$ and $`\widehat{b}_\text{q}=\varphi _m\delta _{\text{q},0}+\delta b_\text{q}`$, where $`\widehat{F}_\text{q}`$ is defined before Eq. (16). Separating out the fluctuation contribution to the Hamiltonian in Eq. (2), we obtain

$``$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}^2}{U_{\mathrm{eff}}}}+{\displaystyle \underset{\text{p}}{}}\xi _\text{p}+{\displaystyle \underset{\text{p}}{}}\widehat{\mathrm{\Psi }}_\text{p}^{}[\xi _\text{p}\tau _3\stackrel{~}{\mathrm{\Delta }}\tau _1]\widehat{\mathrm{\Psi }}_\text{p}+{\displaystyle \underset{\text{q}}{}}\xi _{B\text{q}}b_\text{q}^{}b_\text{q}`$ (18)
$`+`$ $`{\displaystyle \frac{g_\mathrm{r}}{2}}{\displaystyle \underset{\text{q}}{}}[b_\text{q}^{}\rho _\text{q}^{}+b_\text{q}\rho _\text{q}^+]{\displaystyle \frac{U}{4}}{\displaystyle \underset{\text{q}}{}}[\rho _{1\text{q}}\rho _{1\text{q}}+\rho _{2\text{q}}\rho _{2\text{q}}].`$ (19)
Here we have introduced the Nambu field operator for Fermi atoms as $`\widehat{\mathrm{\Psi }}_\text{p}^{}(c_\text{p}^{},c_\text{p})`$ and the corresponding $`2\times 2`$-Pauli matrices $`\tau _i`$ ($`i=1,2,3`$) acting on the particle-hole space. The order parameter $`\stackrel{~}{\mathrm{\Delta }}`$ is defined by Eqs. (6) and (9), which we can take to be real and proportional to $`\tau _1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)`$ without loss of generality. In Eq. (19) and the subsequent discussion, we use the generalized density operators $`\rho _\text{q}^\pm \rho _{1\pm \text{q}}\pm i\rho _{2\pm \text{q}}`$, where

$$\rho _{j\text{q}}\underset{\text{p}}{}\mathrm{\Psi }_{\text{p}+\frac{\text{q}}{2}}^{}\tau _j\mathrm{\Psi }_{\text{p}\frac{\text{q}}{2}}.$$
(20)
We note that $`\rho _{3\text{q}}=_{\text{p},\sigma }c_{\text{p}+\frac{\text{q}}{2},\sigma }^{}c_{\text{p}\frac{\text{q}}{2},\sigma }`$ is the ordinary density fluctuation operator. Similarly, one has

$`\{\begin{array}{c}\rho _{1\text{q}}=_\text{p}[c_{\text{p}\frac{\text{q}}{2}}c_{\text{p}\frac{\text{q}}{2}}+c_{\text{p}+\frac{\text{q}}{2}}^{}c_{\text{p}+\frac{\text{q}}{2}}^{}],\hfill \\ \rho _{2\text{q}}=i_\text{p}[c_{\text{p}\frac{\text{q}}{2}}c_{\text{p}\frac{\text{q}}{2}}c_{\text{p}+\frac{\text{q}}{2}}^{}c_{\text{p}+\frac{\text{q}}{2}}^{}],\hfill \end{array}`$ (23)
and hence

$`\{\begin{array}{c}\rho _\text{q}^+=2_\text{p}c_{\text{p}+\frac{\text{q}}{2}}^{}c_{\text{p}+\frac{\text{q}}{2}}^{},\hfill \\ \rho _\text{q}^{}=2_\text{p}c_{\text{p}\frac{\text{q}}{2}}c_{\text{p}\frac{\text{q}}{2}}.\hfill \end{array}`$ (26)
The operators $`\rho _{1\text{q}}`$ and $`\rho _{2\text{q}}`$ describe, respectively, the amplitude and phase fluctuations of Cooper-pair field fluctuation operator $`\widehat{F}_\text{q}`$. In Eq. (19), the fermion-fermion interaction is seen to be neatly expressed as the sum of interactions between the amplitude fluctuations ($`\frac{U}{4}_\text{q}\rho _{1\text{q}}\rho _{1\text{q}}`$) and phase fluctuations ($`\frac{U}{4}_\text{q}\rho _{2\text{q}}\rho _{2\text{q}}`$). The Feshbach resonance is also expressed as an interaction between the $`b`$-bosons and the fluctuations described by $`\rho _\text{q}^\pm `$. In Eq. (19), we have simply written $`\rho _{1\text{q}=0}\rho _{1\text{q}=0}\rho _{1\text{q}=0}`$ and $`b_{\text{q}=0}\varphi _\mathrm{m}b_{\text{q}=0}`$.

Within the mean-field approximation described by the third term in Eq. (19), the fermion thermal Green’s function is conveniently discussed in terms of a $`2\times 2`$ matrix Green’s function

$$\widehat{G}(𝐩,i\omega _m)=\frac{1}{i\omega _m\xi _\text{p}\tau _3+\stackrel{~}{\mathrm{\Delta }}\tau _1}=\frac{i\omega _m+\xi _\text{p}\tau _3\stackrel{~}{\mathrm{\Delta }}\tau _1}{\omega _m^2+E_\text{p}^2}.$$
(27)
Here $`E_\text{p}\sqrt{\xi _\text{p}^2+\stackrel{~}{\mathrm{\Delta }}^2}`$ is the energy spectrum of fermion quasi-particles below $`T_\mathrm{c}`$, which we shall call the BCS-Bogoliubov quasi-particle spectrum. Eq. (27) reduces to the mean-field BCS matrix Green’s function when the composite order parameter $`\stackrel{~}{\mathrm{\Delta }}`$ is replaced by $`\mathrm{\Delta }`$. The off-diagonal (static) self-energy $`\widehat{\mathrm{\Sigma }}_\mathrm{F}\stackrel{~}{\mathrm{\Delta }}\tau _1`$ comes from the mean-field term $`\stackrel{~}{\mathrm{\Delta }}\tau _1`$ appearing in the Hamiltonian in Eq. (19). This self-energy corresponds to the mean-field diagrams shown in Fig.3,

$`\{\begin{array}{c}\widehat{\mathrm{\Sigma }}_{(\mathrm{a})}=\tau _1{\displaystyle \frac{U}{2\beta }}{\displaystyle \underset{\text{p},\omega _m}{}}\mathrm{Tr}[\tau _1\widehat{G}(𝐩,i\omega _m)]=\tau _1U{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{2E_\text{p}}}\mathrm{tanh}{\displaystyle \frac{\beta }{2}}E_\text{p},\hfill \\ \widehat{\mathrm{\Sigma }}_{(\mathrm{b})}=\tau _+{\displaystyle \frac{g_\mathrm{r}^2}{4\beta }}D_0(0,0){\displaystyle \underset{\text{p},\omega _m}{}}\mathrm{Tr}[\tau _{}\widehat{G}(𝐩,i\omega _m)]=\tau _+{\displaystyle \frac{g_\mathrm{r}^2}{2\nu 2\mu }}{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{4E_\text{p}}}\mathrm{tanh}{\displaystyle \frac{\beta }{2}}E_\text{p},\hfill \\ \widehat{\mathrm{\Sigma }}_{(\mathrm{c})}=\tau _{}{\displaystyle \frac{g_\mathrm{r}^2}{4\beta }}D_0(0,0){\displaystyle \underset{\text{p},\omega _m}{}}\mathrm{Tr}[\tau _+\widehat{G}(𝐩,i\omega _m)]=\tau _{}{\displaystyle \frac{g_\mathrm{r}^2}{2\nu 2\mu }}{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{4E_\text{p}}}\mathrm{tanh}{\displaystyle \frac{\beta }{2}}E_\text{p},\hfill \end{array}`$ (31)
where $`\tau _\pm \tau _1\pm i\tau _2`$. In Eq. (31), $`\widehat{\mathrm{\Sigma }}_{(\mathrm{b})}`$ and $`\widehat{\mathrm{\Sigma }}_{(\mathrm{c})}`$ include the pairing interaction mediated by a Feshbach $`b`$-molecule described by the propagator $`D_0`$. $`\widehat{\mathrm{\Sigma }}_{(\mathrm{a})}`$ comes from the (weak) non-resonant interaction $`U`$. The matrix self-energy $`\widehat{\mathrm{\Sigma }}_\mathrm{F}=\widehat{\mathrm{\Sigma }}_{(\mathrm{a})}+\widehat{\mathrm{\Sigma }}_{(\mathrm{b})}+\widehat{\mathrm{\Sigma }}_{(\mathrm{c})}`$ sums up to give

$$\widehat{\mathrm{\Sigma }}_\mathrm{F}=\tau _1U_{\mathrm{eff}}\underset{\text{p}}{}\frac{\stackrel{~}{\mathrm{\Delta }}}{2E_\text{p}}\mathrm{tanh}\frac{\beta }{2}E_\text{p}=\tau _1\stackrel{~}{\mathrm{\Delta }}.$$
(32)
In the last expression, we have used the gap equation in Eq. (37).

Using eq. (27), we see that

$`\mathrm{\Delta }U{\displaystyle \underset{\text{p}}{}}c_\text{p}c_\text{p}={\displaystyle \frac{U}{2\beta }}{\displaystyle \underset{\text{p},\omega _m}{}}\mathrm{Tr}[\tau _1\widehat{G}(𝐩,i\omega _m)]=U{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{2E_\text{p}}}\mathrm{tanh}{\displaystyle \frac{\beta }{2}}E_\text{p}.`$ (33)
The order parameters $`\varphi _\mathrm{m}`$ and $`\mathrm{\Delta }`$ can be obtained in terms of $`\stackrel{~}{\mathrm{\Delta }}`$ by using Eqs. (8) and (9), to give

$`\{\begin{array}{c}\mathrm{\Delta }={\displaystyle \frac{U}{U_{\mathrm{eff}}}}\stackrel{~}{\mathrm{\Delta }},\hfill \\ \varphi _\mathrm{m}={\displaystyle \frac{g_\mathrm{r}}{2\nu 2\mu }}{\displaystyle \frac{1}{U_{\mathrm{eff}}}}\stackrel{~}{\mathrm{\Delta }}.\hfill \end{array}`$ (36)
The equation for the composite order parameter $`\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }g_\mathrm{r}\varphi _\mathrm{m}`$ can then be rewritten in the form

$$\stackrel{~}{\mathrm{\Delta }}=U_{\mathrm{eff}}\underset{\text{p}}{}\frac{\stackrel{~}{\mathrm{\Delta }}}{2E_\text{p}}\mathrm{tanh}\frac{\beta }{2}E_\text{p},$$
(37)
where $`U_{\mathrm{eff}}`$ is defined in Eq. (11). This self-consistent equation has the same form as the BCS gap equation if we replace $`\stackrel{~}{\mathrm{\Delta }}\mathrm{\Delta }`$ and $`\mu \epsilon _\mathrm{F}`$. We also note that Eq. (37) reduces to the $`T_\mathrm{c}`$ equation in Eq. (10) when $`\stackrel{~}{\mathrm{\Delta }}0`$.

It is important to emphasize that the $`2\times 2`$-matrix single-particle Green’s function in Eq. (27) only includes self-energy effects arising from the (off-diagonal) static mean-fields produced by the Cooper-pairs and the Bose-condensed $`b`$-molecules. In the approximation we use in this paper, the frequency-dependent fermion self-energies associated with the order parameter collective modes are not included in Eq. (27). However, Eq. (27) does implicitly involve the self-consistent renormalized values of $`\stackrel{~}{\mathrm{\Delta }}`$ and $`\mu `$, as determined by the order parameter fluctuations (for further discussion of this kind of approximation, see Ref.). An improved theory of the BCS-BEC crossover would be based on including the fermion self-energies arising from coupling to collective modes.

The chemical potential $`\mu `$ is determined from the equation for the total number of particles $`N`$. As in our discussion of the strong-coupling theory for $`T_\mathrm{c}`$, we work with the thermodynamic potential $`\mathrm{\Omega }`$ consisting of a static mean-field part ($`\mathrm{\Omega }_{\mathrm{MF}}`$) and a fluctuation part ($`\delta \mathrm{\Omega }`$) originating from the particle-particle Cooper-channel, as modified by the Feshbach resonance. The self-consistent equation for $`N`$ is given using the identity $`N=\mathrm{\Omega }/\mu `$. The mean-field part is easily obtained from the first four terms on the right hand side in Eq. (19),

$$\mathrm{\Omega }_{\mathrm{MF}}=\frac{\stackrel{~}{\mathrm{\Delta }}^2}{U_{\mathrm{eff}}}+\underset{\text{p}}{}(\xi _\text{p}E_\text{p})2T\underset{\text{p}}{}\mathrm{ln}\left[1+e^{\beta E_\text{p}}\right]+T\underset{\text{q}}{}\mathrm{ln}\left[1e^{\beta \xi _{B\text{q}}}\right].$$
(38)
Our approximation for the fluctuation contribution $`\delta \mathrm{\Omega }`$ below $`T_\mathrm{c}`$, corresponding to the contribution at $`T_\mathrm{c}`$ shown in Fig. 1(b), is given diagrammatically in Fig. 4(a). Among the interaction terms in Eq. (19), we first carry out perturbative expansion in terms of $`\frac{U}{4}_\text{q}\rho _{2\text{q}}\rho _{2\text{q}}`$, which describes the interaction between the phase fluctuations of $`\mathrm{\Delta }`$. Summing up the loop-type diagrams in Fig. 4(a), we obtain the phase fluctuation contribution to $`\delta \mathrm{\Omega }`$

$$\delta \mathrm{\Omega }_2=\frac{1}{2\beta }\underset{\text{q},\nu _n}{}\mathrm{ln}\left[1+\frac{U}{2}\mathrm{\Pi }_{22}^0(𝐪,i\nu _n)\right].$$
(39)
In this $`\nu _n`$-summation and in the following equations, we omit the important convergence factor $`e^{i\nu _n\delta }`$, for simplicity of notation. The generalized density correlation function $`\mathrm{\Pi }_{22}^0(𝐪,i\nu _n)`$ is defined by

$`\mathrm{\Pi }_{ij}(𝐪,i\nu _n)`$ $`=`$ $`{\displaystyle _0^\beta }d\tau e^{i\nu _n\tau }T_\tau \{\rho _{i\text{q}}(\tau )\rho _{j\text{q}}(0)\}(i,j=1,2,3)`$ (40)
$`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\text{p},\omega _m}{}}\mathrm{Tr}\left[\tau _i\widehat{G}(𝐩+{\displaystyle \frac{𝐪}{2}},i\omega _m+i\nu _n)\tau _j\widehat{G}(𝐩{\displaystyle \frac{𝐪}{2}},i\omega _m)\right],`$ (41)
where the second line ($`\mathrm{\Pi }_{ij}^0`$) is the approximation neglecting the effect of the interactions $`U`$ and $`g_\mathrm{r}`$.

Eq. (41) also defines other correlation functions, which will be important, such as $`\mathrm{\Pi }_{11}`$ and $`\mathrm{\Pi }_{12}`$. Physically, $`\mathrm{\Pi }_{11}`$ and $`\mathrm{\Pi }_{22}`$ describe, respectively, the amplitude and phase fluctuations of Cooper-pairs. $`\mathrm{\Pi }_{33}`$ describes density fluctuations in the gas of Fermi atoms. $`\mathrm{\Pi }_{ij}`$ with $`ij`$ describes a coupling between fluctuations, e.g., $`\mathrm{\Pi }_{12}`$ is a coupling of amplitude fluctuations with the phase fluctuations (amplitude-phase coupling).

Summing up the Matsubara frequencies in Eq. (41), we obtain

$`\mathrm{\Pi }_{11}^0`$ $`=`$ $`{\displaystyle \underset{\text{p}}{}}\left(1{\displaystyle \frac{\xi _{\text{p}+\frac{\text{q}}{2}}\xi _{\text{p}\frac{\text{q}}{2}}\stackrel{~}{\mathrm{\Delta }}^2}{E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}}}}\right){\displaystyle \frac{E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}}}{(E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}})^2+\nu _n^2}}\left[f(E_{\text{p}+\frac{\text{q}}{2}})f(E_{\text{p}\frac{\text{q}}{2}})\right]`$ (42)
$``$ $`{\displaystyle \underset{\text{p}}{}}\left(1+{\displaystyle \frac{\xi _{\text{p}+\frac{\text{q}}{2}}\xi _{\text{p}\frac{\text{q}}{2}}\stackrel{~}{\mathrm{\Delta }}^2}{E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}}}}\right){\displaystyle \frac{E_{\text{p}+\frac{\text{q}}{2}}+E_{\text{p}\frac{\text{q}}{2}}}{(E_{\text{p}+\frac{\text{q}}{2}}+E_{\text{p}\frac{\text{q}}{2}})^2+\nu _n^2}}\left[1f(E_{\text{p}+\frac{\text{q}}{2}})f(E_{\text{p}\frac{\text{q}}{2}})\right],`$ (43)
$`\mathrm{\Pi }_{22}^0`$ $`=`$ $`{\displaystyle \underset{\text{p}}{}}\left(1{\displaystyle \frac{\xi _{\text{p}+\frac{\text{q}}{2}}\xi _{\text{p}\frac{\text{q}}{2}}+\stackrel{~}{\mathrm{\Delta }}^2}{E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}}}}\right){\displaystyle \frac{E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}}}{(E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}})^2+\nu _n^2}}\left[f(E_{\text{p}+\frac{\text{q}}{2}})f(E_{\text{p}\frac{\text{q}}{2}})\right]`$ (44)
$``$ $`{\displaystyle \underset{\text{p}}{}}\left(1+{\displaystyle \frac{\xi _{\text{p}+\frac{\text{q}}{2}}\xi _{\text{p}\frac{\text{q}}{2}}+\stackrel{~}{\mathrm{\Delta }}^2}{E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}}}}\right){\displaystyle \frac{E_{\text{p}+\frac{\text{q}}{2}}+E_{\text{p}\frac{\text{q}}{2}}}{(E_{\text{p}+\frac{\text{q}}{2}}+E_{\text{p}\frac{\text{q}}{2}})^2+\nu _n^2}}\left[1f(E_{\text{p}+\frac{\text{q}}{2}})f(E_{\text{p}\frac{\text{q}}{2}})\right],`$ (45)
$`\mathrm{\Pi }_{33}^0`$ $`=`$ $`{\displaystyle \underset{\text{p}}{}}\left(1+{\displaystyle \frac{\xi _{\text{p}+\frac{\text{q}}{2}}\xi _{\text{p}\frac{\text{q}}{2}}\stackrel{~}{\mathrm{\Delta }}^2}{E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}}}}\right){\displaystyle \frac{E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}}}{(E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}})^2+\nu _n^2}}\left[f(E_{\text{p}+\frac{\text{q}}{2}})f(E_{\text{p}\frac{\text{q}}{2}})\right]`$ (46)
$``$ $`{\displaystyle \underset{\text{p}}{}}\left(1{\displaystyle \frac{\xi _{\text{p}+\frac{\text{q}}{2}}\xi _{\text{p}\frac{\text{q}}{2}}\stackrel{~}{\mathrm{\Delta }}^2}{E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}}}}\right){\displaystyle \frac{E_{\text{p}+\frac{\text{q}}{2}}+E_{\text{p}\frac{\text{q}}{2}}}{(E_{\text{p}+\frac{\text{q}}{2}}+E_{\text{p}\frac{\text{q}}{2}})^2+\nu _n^2}}\left[1f(E_{\text{p}+\frac{\text{q}}{2}})f(E_{\text{p}\frac{\text{q}}{2}})\right].`$ (47)
The correlation functions $`\mathrm{\Pi }_{ij}^0`$ ($`ij`$) describing the coupling of different operators are given by

$`\mathrm{\Pi }_{12}^0`$ $`=`$ $`{\displaystyle \underset{\text{p}}{}}\left({\displaystyle \frac{\xi _{\text{p}+\frac{\text{q}}{2}}}{E_{\text{p}+\frac{\text{q}}{2}}}}{\displaystyle \frac{\xi _{\text{p}\frac{\text{q}}{2}}}{E_{\text{p}\frac{\text{q}}{2}}}}\right){\displaystyle \frac{\nu _n}{(E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}})^2+\nu _n^2}}\left[f(E_{\text{p}+\frac{\text{q}}{2}})f(E_{\text{p}\frac{\text{q}}{2}})\right]`$ (48)
$``$ $`{\displaystyle \underset{\text{p}}{}}\left({\displaystyle \frac{\xi _{\text{p}+\frac{\text{q}}{2}}}{E_{\text{p}+\frac{\text{q}}{2}}}}+{\displaystyle \frac{\xi _{\text{p}\frac{\text{q}}{2}}}{E_{\text{p}\frac{\text{q}}{2}}}}\right){\displaystyle \frac{\nu _n}{(E_{\text{p}+\frac{\text{q}}{2}}+E_{\text{p}\frac{\text{q}}{2}})^2+\nu _n^2}}\left[1f(E_{\text{p}+\frac{\text{q}}{2}})f(E_{\text{p}\frac{\text{q}}{2}})\right],`$ (49)
$`\mathrm{\Pi }_{23}^0`$ $`=`$ $`\stackrel{~}{\mathrm{\Delta }}\nu _n{\displaystyle \underset{\text{p}}{}}\left({\displaystyle \frac{1}{E_{\text{p}+\text{q}/2}}}{\displaystyle \frac{1}{E_{\text{p}\text{q}/2}}}\right){\displaystyle \frac{1}{(E_{\text{p}+\text{q}/2}E_{\text{p}\text{q}/2})^2+\nu _n^2}}\left[f(E_{\text{p}+\text{q}/2})f(E_{\text{p}\text{q}/2})\right]`$ (50)
$`+`$ $`\stackrel{~}{\mathrm{\Delta }}\nu _n{\displaystyle \underset{\text{p}}{}}\left({\displaystyle \frac{1}{E_{\text{p}+\text{q}/2}}}+{\displaystyle \frac{1}{E_{\text{p}\text{q}/2}}}\right){\displaystyle \frac{1}{(E_{\text{p}+\text{q}/2}+E_{\text{p}\text{q}/2})^2+\nu _n^2}}\left[1f(E_{\text{p}+\text{q}/2})f(E_{\text{p}\text{q}/2})\right],`$ (51)
$`\mathrm{\Pi }_{13}^0`$ $`=`$ $`\stackrel{~}{\mathrm{\Delta }}{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{\xi _{\text{p}+\text{q}/2}+\xi _{\text{p}\text{q}/2}}{E_{\text{p}+\text{q}/2}E_{\text{p}\text{q}/2}}}{\displaystyle \frac{E_{\text{p}+\text{q}/2}E_{\text{p}\text{q}/2}}{(E_{\text{p}+\text{q}/2}E_{\text{p}\text{q}/2})^2+\nu _n^2}}\left[f(E_{\text{p}+\text{q}/2})f(E_{\text{p}\text{q}/2})\right]`$ (52)
$``$ $`\stackrel{~}{\mathrm{\Delta }}{\displaystyle \underset{\text{p}}{}}{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{\xi _{\text{p}+\text{q}/2}+\xi _{\text{p}\text{q}/2}}{E_{\text{p}+\text{q}/2}E_{\text{p}\text{q}/2}}}{\displaystyle \frac{E_{\text{p}+\text{q}/2}+E_{\text{p}\text{q}/2}}{(E_{\text{p}+\text{q}/2}+E_{\text{p}\text{q}/2})^2+\nu _n^2}}\left[1f(E_{\text{p}+\text{q}/2})f(E_{\text{p}\text{q}/2})\right],`$ (53)
with $`\mathrm{\Pi }_{21}^0=\mathrm{\Pi }_{12}^0`$, $`\mathrm{\Pi }_{32}^0=\mathrm{\Pi }_{23}^0`$ and $`\mathrm{\Pi }_{31}^0=\mathrm{\Pi }_{13}^0`$. The density-density correlation function $`\mathrm{\Pi }_{33}^0`$ and the related coupling correlation functions $`\mathrm{\Pi }_{13}^0`$ and $`\mathrm{\Pi }_{23}^0`$ will be used in Section VII, where we show that the Goldstone mode describing the collective phase oscillation of Cooper-pairs has spectral weight in the density correlation function. This weight comes from the coupling to the amplitude-density ($`\mathrm{\Pi }_{13}^0`$) and phase-density ($`\mathrm{\Pi }_{23}^0`$) correlation.

The factor $`E_{\text{p}+\frac{\text{q}}{2}}E_{\text{p}\frac{\text{q}}{2}}`$ in the denominator of the first line in Eqs. (43)-(53) describes scattering between excitations with momenta $`𝐩\pm 𝐪/2`$ in the same quasi-particle band $`E_\text{p}`$. For this reason, the first line in Eqs. (43)-(53) is referred to as the intraband term. Since the thermal excitations of fermion quasi-particles are absent at $`T=0`$, the intraband term vanishes at $`T=0`$. On the other hand, the second line in Eqs. (43)-(53) is finite even at $`T=0`$. The factor $`E_{\text{p}+\frac{\text{q}}{2}}+E_{\text{p}\frac{\text{q}}{2}}`$ in the denominator describes interband scattering between $`E_{\text{p}+\frac{\text{q}}{2}}`$ and $`E_{\text{p}\frac{\text{q}}{2}}`$, and thus the second line in Eqs. (43)-(53) is called the interband term. The intraband term is known to give rise to Landau damping of collective modes below the excitation gap $`2\stackrel{~}{\mathrm{\Delta }}`$, while the damping due to the interband term only exists above $`2\stackrel{~}{\mathrm{\Delta }}`$. We note that the gap equation in Eq. (37) can be neatly expressed in terms of $`\mathrm{\Pi }_{22}^0`$ in Eq. (45), namely

$$1+\frac{U_{\mathrm{eff}}}{2}\mathrm{\Pi }_{22}^0(0,0)=0.$$
(54)
The fluctuation contribution ($`\delta \mathrm{\Omega }_1`$) involving the amplitude fluctuations of Cooper-pairs is similarly obtained by summing up the loop-type diagrams in Fig. 4(a) associated with the amplitude-amplitude interaction $`\frac{U}{4}_\text{q}\rho _{1\text{q}}\rho _{1\text{q}}`$. This gives

$$\delta \mathrm{\Omega }_1=\frac{1}{2\beta }\underset{\text{q},\nu _n}{}\mathrm{ln}\left[1+\frac{U}{2}\mathrm{\Pi }_{11}^0(𝐪,i\nu _n)\right].$$
(55)
The phase and amplitude fluctuations are coupled with each other through the amplitude-phase coupling $`\mathrm{\Pi }_{12}^0`$ in Eq. (49). This additional coupling effect can be formally incorporated into $`\delta \mathrm{\Omega }_1`$ by replacing $`\mathrm{\Pi }_{11}^0(𝐪,i\nu _n)`$ in Eq. (55) by $`\overline{\mathrm{\Pi }}_{11}(𝐪,i\nu _n)`$, where

$$\overline{\mathrm{\Pi }}_{11}(𝐪,i\nu _n)\mathrm{\Pi }_{11}^0(𝐪,i\nu _n)+\mathrm{\Pi }_{12}^0(𝐪,i\nu _n)\frac{U/2}{1+(U/2)\mathrm{\Pi }_{22}^0(𝐪,i\nu _n)}\mathrm{\Pi }_{21}^0(𝐪,i\nu _n).$$
(56)
The second term describes the amplitude-phase coupling effect through the coupling correlation functions $`\mathrm{\Pi }_{12}^0`$ and $`\mathrm{\Pi }_{21}^0`$. Eq. (56) is obtained by summing up the fluctuation diagrams shown in Fig. 4(b).

In summary, the fluctuation contribution $`\delta \mathrm{\Omega }_U`$ to the thermodynamic potential involving only the non-resonant interaction $`U`$ is the sum of Eqs. (39) and (56) with $`\mathrm{\Pi }_{11}^0\overline{\mathrm{\Pi }}_{11}`$. The sum is given by

$`\delta \mathrm{\Omega }_U`$ $`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}\mathrm{ln}\left[\left[1+{\displaystyle \frac{U}{2}}\overline{\mathrm{\Pi }}_{11}(𝐪,i\nu _n)\right]\left[1+{\displaystyle \frac{U}{2}}\mathrm{\Pi }_{22}^0(𝐪,i\nu _n)\right]\right]`$ (57)
$`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}\mathrm{ln}\left[\left[1+{\displaystyle \frac{U}{2}}\mathrm{\Pi }_{11}^0(𝐪,i\nu _n)\right]\left[1+{\displaystyle \frac{U}{2}}\mathrm{\Pi }_{22}^0(𝐪,i\nu _n)\right]\left({\displaystyle \frac{U}{2}}\right)^2\mathrm{\Pi }_{12}^0(𝐪,i\nu _n)\mathrm{\Pi }_{21}^0(𝐪,i\nu _n)\right]`$ (58)
$`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}\mathrm{ln}\mathrm{det}\left[1+{\displaystyle \frac{U}{2}}\widehat{\mathrm{\Pi }}^0(𝐪,i\nu _n)\right]`$ (59)
$`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}\mathrm{Tr}\mathrm{ln}\left[1+{\displaystyle \frac{U}{2}}\widehat{\mathrm{\Pi }}^0(𝐪,i\nu _n)\right],`$ (60)
where we have used the well known identity $`det\widehat{A}=e^{\mathrm{Tr}[\mathrm{ln}\widehat{A}]}`$ in the last expression. $`\widehat{\mathrm{\Pi }}^0(𝐪,i\nu _n)`$ is a $`2\times 2`$-matrix density correlation function, defined by

$`\widehat{\mathrm{\Pi }}^0(𝐪,i\nu _n)=\left(\begin{array}{cc}\mathrm{\Pi }_{11}^0& \mathrm{\Pi }_{12}^0\\ \mathrm{\Pi }_{21}^0& \mathrm{\Pi }_{22}^0\end{array}\right).`$ (63)
Finally we consider the fluctuation contributions $`\delta \mathrm{\Omega }_{\mathrm{FR}}`$ from the Feshbach resonance term $`\frac{1}{2}g_\mathrm{r}_\text{q}[b_\text{q}^{}\rho _\text{q}^{}+b_\text{q}\rho _\text{q}^+]`$ in Eq. (19). First we sum up the diagrams described in Fig. 4(a), where the dashed line now represents the $`b`$-boson Green’s function $`D_0(𝐪,i\nu _n)`$ in Eq. (13), to give

$$\delta \mathrm{\Omega }_{\mathrm{FR}}=\frac{1}{2\beta }\underset{\text{q},\nu _n}{}\mathrm{Tr}\mathrm{ln}\left[1\frac{g_\mathrm{r}^2}{4}\widehat{D}^0(𝐪,i\nu _n)\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)\right].$$
(64)
Here $`\widehat{D}^0(𝐪,i\nu _n)`$ is a $`2\times 2`$-matrix $`b`$-boson Green’s function defined by

$`\widehat{D}^0(𝐪,\nu _n)={\displaystyle \frac{1}{i\nu _n\tau _3\xi _{B\text{q}}}}={\displaystyle \frac{i\nu _n\tau _3+\xi _{B\text{q}}}{\nu _n^2+\xi _{B\text{q}}^2}},`$ (65)
and $`\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)`$ is a $`2\times 2`$-matrix correlation function (neglecting the effect of $`U`$ and $`g_\mathrm{r}`$) defined by

$`\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)={\displaystyle _0^\beta }𝑑\tau e^{i\nu _n\tau }T_\tau \left[\left(\begin{array}{cc}\rho _\text{q}^{}(\tau )\rho _\text{q}^+(0)& \rho _\text{q}^{}(\tau )\rho _\text{q}^{}(0)\\ \rho _\text{q}^+(\tau )\rho _\text{q}^+(0)& \rho _\text{q}^+(\tau )\rho _\text{q}^{}(0)\end{array}\right)\right].`$ (68)
Using the definition $`\rho _\text{q}^\pm =\rho _{1\pm \text{q}}\pm i\rho _{2\pm \text{q}}`$, we find that $`\widehat{\mathrm{\Xi }}^0`$ can be expressed in terms of the matrix elements of $`\widehat{\mathrm{\Pi }}^0`$ as follows:

$`\widehat{\mathrm{\Xi }}^0`$ $`=`$ $`2\widehat{W}^1\widehat{\mathrm{\Pi }}^0\widehat{W}`$ (69)
$`=`$ $`\left(\begin{array}{cc}\mathrm{\Pi }_{11}^0+\mathrm{\Pi }_{22}^0+i(\mathrm{\Pi }_{12}^0\mathrm{\Pi }_{21}^0)& \mathrm{\Pi }_{11}^0\mathrm{\Pi }_{22}^0\\ \mathrm{\Pi }_{11}^0\mathrm{\Pi }_{22}^0& \mathrm{\Pi }_{11}^0+\mathrm{\Pi }_{22}^0i(\mathrm{\Pi }_{12}^0\mathrm{\Pi }_{21}^0)\end{array}\right),`$ (72)
where $`\widehat{W}`$ is the unitary matrix

$`\widehat{W}={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{cc}1& 1\\ i& i\end{array}\right).`$ (75)
Next we renormalize the fluctuations in Eq. (64) by including the effects of the non-resonant interaction $`U`$ on the correlation function $`\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)`$, working within the HF-RPA. This results in $`\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)`$ in Eq. (64) being replaced by $`\widehat{\mathrm{\Xi }}_U`$, where

$$\widehat{\mathrm{\Xi }}_U(𝐪,i\nu _n)=\left[1+\frac{U}{4}\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)\right]^1\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n).$$
(76)
The resulting expression for $`\delta \mathrm{\Omega }_{\mathrm{FR}}`$ involves the fluctuation effects related to both the non-resonant interaction $`U`$ and the Feshbach resonance coupling parameter $`g_\mathrm{r}`$, namely

$`\delta \mathrm{\Omega }_{\mathrm{FR}}`$ $`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}\mathrm{Tr}\mathrm{ln}\left[1{\displaystyle \frac{g_\mathrm{r}^2}{4}}\widehat{D}^0(𝐪,i\nu _n)\widehat{\mathrm{\Xi }}_U(𝐪,i\nu _n)\right]`$ (77)
$`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}\mathrm{Tr}\mathrm{ln}\left\{1{\displaystyle \frac{g_\mathrm{r}^2}{4}}\widehat{D}^0(𝐪,i\nu _n)\left[1+{\displaystyle \frac{U}{4}}\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)\right]^1\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)\right\}.`$ (78)
The total fluctuation contribution $`\delta \mathrm{\Omega }`$ is given by the sum of Eqs. (60) and (78). Recalling the definition in Eq.(72) and the relation $`\mathrm{Tr}\mathrm{ln}[1+\frac{U}{2}\widehat{\mathrm{\Pi }}^0]=\mathrm{Tr}\mathrm{ln}[\widehat{W}\{1+\frac{U}{4}\widehat{\mathrm{\Xi }}^0\}\widehat{W}^1]=\mathrm{Tr}\mathrm{ln}[1+\frac{U}{4}\widehat{\mathrm{\Xi }}^0]`$, this sum can be written as

$`\delta \mathrm{\Omega }`$ $``$ $`\delta \mathrm{\Omega }_U+\delta \mathrm{\Omega }_{\mathrm{FR}}`$ (79)
$`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}\mathrm{Tr}\mathrm{ln}\left\{1+{\displaystyle \frac{1}{4}}\left[Ug_\mathrm{r}^2\widehat{D}^0(𝐪,i\nu _n)\right]\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)\right\}.`$ (80)
Putting everything together, the total thermodynamic potential $`\mathrm{\Omega }=\mathrm{\Omega }_{\mathrm{MF}}+\delta \mathrm{\Omega }`$ is

$`\mathrm{\Omega }`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}^2}{U_{\mathrm{eff}}}}+{\displaystyle \underset{\text{p}}{}}(\xi _\text{p}E_\text{p})2T{\displaystyle \underset{\text{p}}{}}\mathrm{ln}\left[1+e^{\beta E_\text{p}}\right]+T{\displaystyle \underset{\text{q}}{}}\mathrm{ln}\left[1e^{\beta \xi _\text{q}^B}\right]`$ (81)
$`+`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}\mathrm{Tr}\mathrm{ln}\left\{1+{\displaystyle \frac{1}{4}}\left[Ug_\mathrm{r}^2\widehat{D}^0(𝐪,i\nu _n)\right]\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)\right\}.`$ (82)
The equation for the total number of particles is then obtained from $`N=\mathrm{\Omega }/\mu `$.

In taking the derivative with respect to $`\mu `$, we note that the order parameter $`\stackrel{~}{\mathrm{\Delta }}`$ also depends on the chemical potential $`\mu `$, in addition to $`\xi _\text{p}\epsilon _\text{p}\mu `$ and $`\xi _{B\text{q}}\epsilon _{B\text{q}}+2\nu 2\mu `$. Thus one needs to calculate

$`{\displaystyle \frac{\mathrm{\Omega }}{\stackrel{~}{\mathrm{\Delta }}}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{\mu }}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_{\mathrm{MF}}}{\stackrel{~}{\mathrm{\Delta }}}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{\mu }}+{\displaystyle \frac{\delta \mathrm{\Omega }}{\stackrel{~}{\mathrm{\Delta }}}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{\mu }}`$ (83)
$`=`$ $`{\displaystyle \frac{\delta \mathrm{\Omega }}{\stackrel{~}{\mathrm{\Delta }}}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{\mu }}.`$ (84)
In the last line, we have used the fact $`\mathrm{\Omega }_{\mathrm{MF}}/\stackrel{~}{\mathrm{\Delta }}=0`$, which holds when $`\stackrel{~}{\mathrm{\Delta }}`$ satisfies the gap equation (37). However, as will be shown in Section IV, since $`\delta \mathrm{\Omega }/\stackrel{~}{\mathrm{\Delta }}`$ in Eq. (84) is a higher order correction within the perturbative approximation we are using, we can neglect the contribution from Eq. (84). Thus the dependence of $`\stackrel{~}{\mathrm{\Delta }}`$ on $`\mu `$ only leads to higher order corrections. The resulting equation for $`N`$ is

$`N`$ $`=`$ $`2\varphi _\mathrm{m}^2+{\displaystyle \underset{\text{p}}{}}\left[1{\displaystyle \frac{\xi _\text{p}}{E_\text{p}}}\mathrm{tanh}{\displaystyle \frac{\beta }{2}}E_\text{p}\right]+2{\displaystyle \underset{\text{q}}{}}n_\mathrm{B}(\xi _\text{q}^B)`$ (85)
$``$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}{\displaystyle \frac{}{\mu }}\mathrm{Tr}\mathrm{ln}\left\{1+{\displaystyle \frac{1}{4}}\left[Ug_\mathrm{r}^2\widehat{D}^0(𝐪,i\nu _n)\right]\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)\right\}.`$ (86)
Here it is understood that the $`\mu `$-derivative in the last term only acts on the chemical potential involved in $`\xi _\text{p}`$ and $`\xi _{B\text{q}}`$. Eqs. (37) and (86) give us the required self-consistent coupled equations for $`\stackrel{~}{\mathrm{\Delta }}`$ and $`\mu `$ in the superfluid phase below $`T_\mathrm{c}`$. We will discuss our numerical self-consistent solutions of the coupled equations (37) and (86) in Section VI.

When one neglects the fluctuation contribution given in the last term in Eq. (86), we obtain the mean-field expression obtained in Ref. in the context of high-$`T_\mathrm{c}`$ superconductivity. We also mention that Eq. (86) reproduces Eq. (12) for the normal phase $`TT_\mathrm{c}`$, where $`\mathrm{\Delta }=\varphi _\mathrm{m}=0`$. To see this, we note that the second term on the right hand side in Eq. (86) reduces to $`2N_\mathrm{F}^0`$ in Eq. (12). In addition, since the phase and amplitude fluctuations are indistinguishable when $`\mathrm{\Delta }=0`$, we find $`\mathrm{\Pi }_{11}^0=\mathrm{\Pi }_{22}^0`$, and thus $`\widehat{\mathrm{\Xi }}_0`$ in Eq. (72) becomes diagonal at $`T_\mathrm{c}`$. Noting that $`\mathrm{\Xi }_{11}^0(𝐪,i\nu _n,T=T_\mathrm{c})=4\mathrm{\Pi }(𝐪,i\nu _n)`$ and $`\mathrm{\Xi }_{22}^0(𝐪,i\nu _n,T=T_\mathrm{c})=4\mathrm{\Pi }(𝐪,i\nu _n)`$, where $`\mathrm{\Pi }(𝐪,i\nu _n)`$ is given in Eq. (16), we find that the last term in Eq. (86) ($`\delta N`$) reproduces the last term in Eq. (12). More explicitly, we have

$`\delta N(T_\mathrm{c})`$ $`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}{\displaystyle \frac{}{\mu }}\mathrm{ln}[(1[Ug_\mathrm{r}^2D_{11}^0(𝐪,i\nu _n)]\mathrm{\Pi }(𝐪,i\nu _n))`$ (88)
$`\times (1[Ug_\mathrm{r}^2D_{22}^0(𝐪,i\nu _n)]\mathrm{\Pi }(𝐪,i\nu _n))]`$
$`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}{\displaystyle \frac{}{\mu }}\mathrm{ln}\left\{1[Ug_\mathrm{r}^2D_0(𝐪,i\nu _n)]\mathrm{\Pi }(𝐪,i\nu _n)\right\},`$ (89)
where we have used $`D_{11}^0(𝐪,i\nu _n)=D_0(𝐪,i\nu _n)`$ and $`D_{22}^0(𝐪,i\nu _n)=D_0(𝐪,i\nu _n)`$ in the last line. Thus Eq. (86) is equivalent to Eq. (12) at $`T_\mathrm{c}`$. Our present strong-coupling theory giving the coupled equations (37) and (86) for the superfluid phase is seen to smoothly go over to our previous discussion at $`T_\mathrm{c}`$ and above.

Each term in Eq. (86) has a simple physical meaning, which is useful to discuss. The first term

$$2N_\mathrm{B}^\mathrm{c}2\varphi _\mathrm{m}^2=2b_{\text{q}=0}^2$$
(90)
gives twice the number of Bose-condensed $`b`$-bosons. The second term

$$N_\mathrm{F}\underset{\text{p}}{}[1\frac{\xi _\text{p}}{E_\text{p}}\mathrm{tanh}\frac{\beta }{2}E_\text{p}]$$
(91)
describes the number of Fermi quasi-particles. This expression can be directly obtained from $`N_\mathrm{F}=_{\text{p},\sigma }c_{\text{p}\sigma }^{}c_{\text{p}\sigma }`$ in Eq. (3). To understand the physical meanings of the last two terms in Eq. (86), it is convenient to divide the $`\mu `$-derivative in the last term into the derivative acting on $`\xi _{B\text{q}}=q^2/2M+2\nu 2\mu `$ in the $`b`$-boson Green’s function $`\widehat{D}^0(𝐪,i\nu _n)`$ and the derivative acting on $`\xi _\text{p}=\epsilon _\text{p}\mu `$ in $`\mathrm{\Xi }^0(𝐪,i\nu _n)`$ ($`/\mu _\mathrm{F}`$). Using the identity

$$2\underset{\text{q}}{}n_\mathrm{B}(\xi _{B\text{q}})=\underset{\text{q}}{}\left[1+\frac{1}{\beta }\underset{\nu _n}{}\mathrm{Tr}[\widehat{D}^0(𝐪,i\nu _n)]\right],$$
(92)
we can write Eq. (86) as

$$N=2N_\mathrm{B}^\mathrm{c}+N_\mathrm{F}+2N_\mathrm{B}^\mathrm{n}+2N_\mathrm{C},$$
(93)
where $`N_\mathrm{B}^\mathrm{n}`$ and $`N_\mathrm{C}`$ are defined by

$`\{\begin{array}{c}2N_\mathrm{B}^\mathrm{n}{\displaystyle \underset{\text{q}}{}}\left[1+{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\nu _n}{}}\mathrm{Tr}[\widehat{D}(𝐪,i\nu _n)]\right],\hfill \\ 2N_\mathrm{C}{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}{\displaystyle \frac{}{\mu _\mathrm{F}}}\mathrm{Tr}\mathrm{ln}\left\{1+{\displaystyle \frac{1}{4}}\left[Ug_\mathrm{r}^2\widehat{D}^0(𝐪,i\nu _n)\right]\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)\right\}.\hfill \end{array}`$ (96)
Here $`\widehat{D}(𝐪,i\nu _n)`$ is a renormalized $`2\times 2`$-matrix thermal $`b`$-boson Green’s function defined by

$$\widehat{D}(𝐪,i\nu _n)=\frac{1}{i\nu _n\tau _3\xi _{B\text{q}}\widehat{\mathrm{\Sigma }}(𝐪,i\nu _n)}.$$
(97)
The $`b`$-boson matrix self-energy

$$\widehat{\mathrm{\Sigma }}(𝐪,i\nu _n)\frac{g_\mathrm{r}^2}{4}\widehat{\mathrm{\Xi }}_U(𝐪,i\nu _n)=\frac{g_\mathrm{r}^2}{4}\left[1+\frac{U}{4}\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)\right]^1\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n),$$
(98)
describes the Feshbach coupling with Fermi atoms.

The renormalized $`b`$-boson Green’s function in Eq. (97) has the same form as that obtained by Kostyrko and Ranninger calculated within the HF-RPA. Comparing $`N_\mathrm{B}^\mathrm{n}`$ in Eq. (96) with Eq. (92), we may interpret $`N_\mathrm{B}^\mathrm{n}`$ as the number of non-condensed $`b`$-bosons, as renormalized by the Feshbach resonance. In analogy to our previous discussion of the strong-coupling theory for $`T_\mathrm{c}`$, $`N_\mathrm{C}`$ in Eq. (96) may be understood as the fluctuation contribution to $`N`$ from Cooper-pairs associated with the (dynamical) pairing interaction given by

$$\widehat{U}_{\mathrm{eff}}(𝐪,i\nu _n)Ug_\mathrm{r}^2\widehat{D}^0(𝐪,i\nu _n).$$
(99)
We note that $`\widehat{U}_{\mathrm{eff}}(𝐪=0,i\nu _n=0)=U_{\mathrm{eff}}\widehat{\mathrm{𝟏}}`$, where $`U_{\mathrm{eff}}=U+g_\mathrm{r}^2/(2\nu 2\mu )`$ is the effective pairing interaction appearing in the gap equation in Eq. (37).

The renormalized $`b`$-boson Green’s function in Eq. (97) can be also written in the form

$`\widehat{D}(𝐪,i\nu _n)`$ $`=`$ $`{\displaystyle \frac{1}{i\stackrel{~}{\nu }_n\tau _3\stackrel{~}{\xi }_{B\text{q}}\kappa \tau _1}}={\displaystyle \frac{i\stackrel{~}{\nu }_n\tau _3+\stackrel{~}{\xi }_{B\text{q}}\kappa \tau _1}{\stackrel{~}{\nu }_n^2+\stackrel{~}{\xi }_{B\text{q}}^2\kappa ^2}},`$ (100)
where the renormalized parameters are given by

$`\{\begin{array}{c}i\stackrel{~}{\nu }_ni\nu _ni{\displaystyle \frac{g_\mathrm{r}^2}{4}}[\mathrm{\Pi }_{U12}\mathrm{\Pi }_{U21}],\hfill \\ \stackrel{~}{\xi }_{B\text{q}}\xi _{B\text{q}}+{\displaystyle \frac{g_\mathrm{r}^2}{4}}[\mathrm{\Pi }_{U11}+\mathrm{\Pi }_{U22}],\hfill \\ \kappa {\displaystyle \frac{g_\mathrm{r}^2}{4}}[\mathrm{\Pi }_{U11}\mathrm{\Pi }_{U22}].\hfill \end{array}`$ (104)
The correlation functions $`\mathrm{\Pi }_{Uij}`$ ($`i,j=1,2`$) are the matrix elements of $`\widehat{\mathrm{\Pi }}_U`$ defined by

$$\widehat{\mathrm{\Pi }}_U(𝐪,i\nu _n)=\left[1+\frac{U}{2}\widehat{\mathrm{\Pi }}^0(𝐪,i\nu _n)\right]^1\widehat{\mathrm{\Pi }}^0(𝐪,i\nu _n).$$
(105)
More explicitly, we have

$`\{\begin{array}{c}\mathrm{\Pi }_{U11}={\displaystyle \frac{\overline{\mathrm{\Pi }}_{11}}{1+\frac{U}{2}\overline{\mathrm{\Pi }}_{11}}},\hfill \\ \mathrm{\Pi }_{U22}={\displaystyle \frac{\overline{\mathrm{\Pi }}_{22}}{1+\frac{U}{2}\overline{\mathrm{\Pi }}_{22}}},\hfill \\ \mathrm{\Pi }_{U12}={\displaystyle \frac{\mathrm{\Pi }_{12}^0}{[1+\frac{U}{2}\mathrm{\Pi }_{11}^0][1+\frac{U}{2}\mathrm{\Pi }_{22}^0]\frac{U^2}{4}\mathrm{\Pi }_{12}^0\mathrm{\Pi }_{21}^0}},\hfill \\ \mathrm{\Pi }_{U21}(𝐪,i\nu _n)=\mathrm{\Pi }_{U12}(𝐪,i\nu _n).\hfill \end{array}`$ (110)
Here $`\overline{\mathrm{\Pi }}_{11}`$ is defined in Eq. (56), and $`\overline{\mathrm{\Pi }}_{22}`$ is similarly obtained from Eq. (56) by interchanging $`12`$.

When we take $`𝐪=i\nu _n=0`$, the denominator of Eq. (100) reduces to

$`\stackrel{~}{\xi }_{B\text{q}=0}^2\kappa (0,0)^2=[1+{\displaystyle \frac{U_{\mathrm{eff}}}{2}}\mathrm{\Pi }_{22}^0(0,0)][1+{\displaystyle \frac{U_{\mathrm{eff}}}{2}}\mathrm{\Pi }_{11}^0(0,0)]{\displaystyle \frac{2\nu 2\mu }{1+\frac{U}{2}\mathrm{\Pi }_{22}^0(0,0)}}{\displaystyle \frac{2\nu 2\mu }{1+\frac{U}{2}\mathrm{\Pi }_{11}^0(0,0)}}.`$ (111)
The expression in Eq. (111) clearly vanishes when the gap equation in Eq. (54) is satisfied. This means that the excitation spectrum of the renormalized $`b`$-bosons described by $`\widehat{D}`$ is always gapless at $`𝐪=0`$ for $`T<T_\mathrm{c}`$. This is a desired result, because the Bose-condensation of $`b`$-bosons characterized by $`\varphi _\mathrm{m}`$ should be accompanied by the appearance of a Bogoliubov phonon (Goldstone) mode having a (uniform system) gapless dispersion. Thus the strong-coupling theory presented in this section correctly includes a gapless spectrum for the symmetry breaking Goldstone mode. This result is not obtained if we use a static mean-field theory, neglecting the fluctuation contribution given by the last term of Eq. (86). Within such a mean-field theory, the $`b`$-boson excitation spectrum is given by $`\epsilon _{B\text{q}}+2\nu 2\mu `$, which always has a finite excitation gap at $`𝐪=0`$. This gap is given by $`\delta E2\nu 2\mu =\frac{g_\mathrm{r}^2}{2}\mathrm{\Pi }_{U22}(0,0)>0`$, assuming that the gap equation in Eq. (54) is satisfied. We note that $`\mathrm{\Pi }_{U22}(0,0)`$ is defined in Eq. (110), with $`\overline{\mathrm{\Pi }}_{22}(0,0)=\mathrm{\Pi }_{22}^0(0,0)`$ since $`\mathrm{\Pi }_{12}^0(0,0)=0`$.

The gapless behavior of the renormalized $`b`$-boson can be also verified formally by noting that the Hugenholtz-Pines theorem is verified, namely

$$\mathrm{\Sigma }_{11}(0,0)\mathrm{\Sigma }_{12}(0,0)=\frac{g_\mathrm{r}^2}{2}\frac{\mathrm{\Pi }_{22}^0(0,0)}{1+\frac{U}{2}\mathrm{\Pi }_{22}^0(0,0)}=2\mu 2\nu .$$
(112)
Here we have used the gap equation given in Eq. (54).

In order to discuss BCS-BEC crossover in the superfluid phase, one needs to solve the coupled equations (37) and (86) numerically. Results will be discussed in Section VI. Here we briefly discuss limiting cases which can be treated analytically.

When the threshold energy $`2\nu `$ of the $`b`$-boson excitation spectrum is very large (BCS limit), $`2\nu 2\mu `$ can be realized because the chemical potential $`\mu `$ is at most the order of $`\epsilon _\mathrm{F}`$. In this case, the Fermi atoms are dominant particles while the effects of $`b`$-molecules described by $`\varphi _\mathrm{m}`$, $`n_\mathrm{B}(\xi _{B\text{q}})`$ and $`\widehat{D}^0`$ in Eq. (86) can be neglected. In addition, since the last term in Eq. (86) is small when one is dealing with a weak non-resonant interaction $`U`$, one can also drop this term. Thus we see that Eq. (86) reduces to $`N=N_\mathrm{F}`$ (where $`N_\mathrm{F}`$ is defined in Eq. (91)), which simply gives $`\mu \epsilon _\mathrm{F}`$ as far as $`\epsilon _\mathrm{F}T`$. When this result is substituted into the gap equation in Eq. (37) with $`U_{\mathrm{eff}}U`$ and $`\stackrel{~}{\mathrm{\Delta }}\mathrm{\Delta }`$ (note that $`2\nu 2\mu `$.), one obtains the usual BCS gap equation for the Cooper-pair order parameter $`\mathrm{\Delta }`$.

In the opposite limit $`2\nu <2\epsilon _\mathrm{F}`$ (BEC limit), since the $`b`$-boson branch has an energy lower than the two fermion band energy, most Fermi atoms will combine to form $`b`$-molecules and hence the fermion correlation functions $`\mathrm{\Pi }_{ij}`$ becomes less important. Then the gap equation in Eq. (54) can be rewritten as \[see also Eq. (112)\]

$$2\mu =2\nu +\frac{g_\mathrm{r}^2}{2}\frac{\mathrm{\Pi }_{22}^0(0,0)}{1+\frac{U}{2}\mathrm{\Pi }_{22}^0(0,0)}2\nu .$$
(113)
This says that the chemical potential has the energy of the bottom of the $`b`$-boson excitation spectrum. Substituting $`2\mu =2\nu `$ into Eq. (86) with $`\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)=0`$, we obtain

$$\frac{N}{2}=\varphi _\mathrm{m}+\underset{\text{q}}{}n_\mathrm{B}(\frac{q^2}{2M}).$$
(114)
This is just the equation for BEC in an non-interacting uniform Bose gas with $`N/2`$ bosons of mass $`M`$. Thus the coupled Eqs. (37) and (86) reproduce both the BCS-phase and BEC-phase for two limiting values of the Feshbach molecular resonance threshold $`2\nu `$.

We again remind the reader about the many body approximation our whole discussion is based on. As we noted after Eq. (37), we have included the fluctuations of the order parameters $`\mathrm{\Delta }`$ and $`\varphi _\mathrm{m}`$ in determining self-consistently the values of $`N`$ and $`\mu `$. However, we have not included the explicit self-energy of the fermions due to these dynamic fluctuations. For further discussion, see Refs..

## IV Relation to a Gaussian fluctuation theory

It is interesting see what type of fluctuations are taken into account in the strong-coupling theory presented in Section III. In this Section, we show that the coupled Eqs. (37) and (86) describe Gaussian fluctuations around the mean-field order parameters $`\mathrm{\Delta }`$ and $`\varphi _\mathrm{m}`$. To see this, we employ the functional integral formalism because it treats fluctuations in a very clear fashion. Indeed, using this formalism, Randeria and co-workers showed that the strong-coupling superconductivity theory developed by Nozières and Schmitt-Rink can be understood as a Gaussian approximation for the fluctuations in the Cooper-channel. Very recently, Milstein et al. also employed this approach to the coupled fermion-boson model in Eq. (2) and reproduced the equations for $`T_\mathrm{c}`$ given in Eqs. (10) and (12).

The partition function for the coupled fermion-boson model in Eq. (2) in the functional integral formalism is given by

$$Z=𝒟\mathrm{\Psi }_\sigma ^{}𝒟\mathrm{\Psi }_\sigma 𝒟\mathrm{\Phi }^{}𝒟\mathrm{\Phi }e^S.$$
(115)
The action $`S`$ has the form

$`S`$ $`=`$ $`{\displaystyle _0^\beta }d\tau {\displaystyle }d𝐫[{\displaystyle \underset{\sigma }{}}\mathrm{\Psi }_\sigma ^{}(x)[_\tau +{\displaystyle \frac{\widehat{𝐩}^2}{2m}}\mu ]\mathrm{\Psi }_\sigma (x)U\mathrm{\Psi }_{}^{}(x)\mathrm{\Psi }_{}^{}(x)\mathrm{\Psi }_{}(x)\mathrm{\Psi }_{}(x)`$ (116)
$`+`$ $`\mathrm{\Phi }^{}(x)[_\tau +{\displaystyle \frac{\widehat{𝐩}^2}{2M}}+2\nu 2\mu ]\mathrm{\Phi }(x)+g_\mathrm{r}[\mathrm{\Phi }^{}(x)\mathrm{\Psi }_{}(x)\mathrm{\Psi }_{}(x)+\mathrm{\Phi }(x)\mathrm{\Psi }_{}^{}(x)\mathrm{\Psi }_{}^{}(x)]],`$ (117)
where $`\widehat{𝐩}/i`$ and $`x(𝐫,\tau )`$. $`\mathrm{\Psi }_\sigma (x)`$ is a Grassmann variable for Fermi atoms and $`\mathrm{\Psi }_\sigma ^{}(x)`$ is its conjugate. $`\mathrm{\Phi }(x)`$ is a complex scalar Bose field describing $`b`$-molecules. When we introduce the Cooper-pair field $`\mathrm{\Delta }`$ using the Stratonovich-Hubbard transformation, the partition function in Eq. (115) is transformed to

$`Z={\displaystyle 𝒟\mathrm{\Psi }_\sigma ^{}𝒟\mathrm{\Psi }_\sigma 𝒟\mathrm{\Phi }^{}𝒟\mathrm{\Phi }𝒟\mathrm{\Delta }^{}𝒟\mathrm{\Delta }e^{\overline{S}}},`$ (118)
where the action $`\overline{S}`$ is now given by

$`\overline{S}`$ $`=`$ $`{\displaystyle _0^\beta }𝑑\tau {\displaystyle 𝑑𝐫\left[\frac{|\mathrm{\Delta }(x)|^2}{U}\widehat{\mathrm{\Psi }}^{}(x)\widehat{G}^1(x)\widehat{\mathrm{\Psi }}(x)\frac{1}{2}\widehat{\mathrm{\Phi }}^{}(x)\widehat{D}^0(x)^1\widehat{\mathrm{\Phi }}(x)\right]}`$ (119)
$`+`$ $`\beta {\displaystyle \underset{\text{p}}{}}\xi _\text{p}{\displaystyle \frac{\beta }{2}}{\displaystyle \underset{\text{q}}{}}\xi _{B\text{q}}.`$ (120)
Here we have introduced the Nambu representation for fermions as $`\widehat{\mathrm{\Psi }}^{}(\mathrm{\Psi }_{}^{},\mathrm{\Psi }_{})`$. The corresponding $`2\times 2`$-matrix thermal Green’s function $`\widehat{G}(x)`$ in the $`(𝐫,\tau )`$-representation is given by

$`\widehat{G}^1(x)=_\tau \left({\displaystyle \frac{\widehat{𝐩}^2}{2m}}\mu \right)\tau _3+\widehat{\stackrel{~}{\mathrm{\Delta }}}(x),`$ (121)
where

$`\widehat{\stackrel{~}{\mathrm{\Delta }}}(x)\left(\begin{array}{cc}0& \stackrel{~}{\mathrm{\Delta }}(x)\\ \stackrel{~}{\mathrm{\Delta }}(x)^{}& 0\end{array}\right)=\left(\begin{array}{cc}0& \mathrm{\Delta }(x)g_\mathrm{r}\mathrm{\Phi }(x)\\ \mathrm{\Delta }(x)^{}g_\mathrm{r}\mathrm{\Phi }^{}(x)& 0\end{array}\right).`$ (126)
We have also introduced a matrix notation for the $`b`$-bosons as $`\widehat{\mathrm{\Phi }}^{}(x)(\mathrm{\Phi }^{}(x),\mathrm{\Phi }(x))`$, as well as the corresponding matrix $`b`$-boson Green’s function in the $`(𝐫,\tau )`$-representation,

$$\widehat{D}^0(x)^1=_\tau \tau _3\left(\frac{\widehat{𝐩}^2}{2M}+2\nu 2\mu \right).$$
(127)
Since the action $`\overline{S}`$ in Eq. (120) has a bi-linear form with respect to $`\widehat{\mathrm{\Psi }}`$ and $`\widehat{\mathrm{\Psi }}^{}`$, one can immediately integrate out the fermion degrees of freedom using the formula

$$𝒟\mathrm{\Psi }_\sigma ^{}𝒟\mathrm{\Psi }_\sigma e^{_{ij}\mathrm{\Psi }_i^{}M_{ij}\mathrm{\Psi }_j}=\mathrm{det}[\widehat{M}],$$
(128)
where $`i`$ and $`j`$ represent the variables ($`𝐫,\tau ,\sigma `$) involved in $`\mathrm{\Psi }_\sigma (𝐫,\tau )`$. After the functional integration, the action is reduced to ($`\stackrel{~}{S}`$)

$`\stackrel{~}{S}`$ $`=`$ $`\mathrm{Tr}\mathrm{ln}[\widehat{G}^1]+{\displaystyle _0^\beta }𝑑\tau {\displaystyle 𝑑𝐫\left[\frac{|\mathrm{\Delta }(x)|^2}{U}\frac{1}{2}\widehat{\mathrm{\Phi }}(x)^{}\widehat{D}^0(x)^1\widehat{\mathrm{\Phi }}(x)\right]}`$ (129)
$`+`$ $`\beta {\displaystyle \underset{\text{p}}{}}\xi _\text{p}{\displaystyle \frac{\beta }{2}}{\displaystyle \underset{\text{q}}{}}\xi _{B\text{q}}.`$ (130)
The trace in $`\mathrm{Tr}\mathrm{ln}[\widehat{G}^1]`$ involves the integrals over $`\tau `$ and $`𝐫`$, as well as a trace over the particle-hole space. Eq. (130) has a very physical form because it is described by two kinds of Bose fields $`\mathrm{\Delta }(x)`$ and $`\mathrm{\Phi }(x)`$ directly related to superfluid fluctuations.

In order to extract the fluctuation contribution from Eq. (130), we divide the Bose fields $`\mathrm{\Delta }(x)`$ and $`\mathrm{\Phi }(x)`$ into their mean-field values ($`\mathrm{\Delta }`$, $`\varphi _\mathrm{m}`$) and fluctuations from these values. In the functional integral formalism, the mean-field approximation corresponds to the saddle point solution, given by

$$\frac{\delta \stackrel{~}{S}}{\delta \mathrm{\Delta }^{}(x)}=\frac{\delta \stackrel{~}{S}}{\delta \mathrm{\Phi }^{}(x)}=0.$$
(131)
Since we are working with a uniform system (no trap), we can take a uniform (and real) solution as the saddle point solution. Then Eq. (131) gives the two equations:

$`\{\begin{array}{c}\mathrm{\Delta }=U\widehat{G}_{12}(x)={\displaystyle \frac{U}{\beta }}{\displaystyle \underset{\text{p},\omega _m}{}}\widehat{G}_{12}(𝐩,i\omega _m)={\displaystyle \frac{U}{\beta }}{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{2E_\text{p}}}\mathrm{tanh}{\displaystyle \frac{\beta }{2}}E_\text{p},\hfill \\ (2\nu 2\mu )\varphi _\mathrm{m}=g_\mathrm{r}\widehat{G}_{12}(x)={\displaystyle \frac{g_\mathrm{r}}{\beta }}{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{2E_\text{p}}}\mathrm{tanh}{\displaystyle \frac{\beta }{2}}E_\text{p},\hfill \end{array}`$ (134)
where $`\widehat{G}_{12}(𝐩,i\omega _m)`$ is the off-diagonal component of the fermion Green’s function in Eq. (27). Eq. (134) is equivalent to Eqs. (33)-(36).

Writing $`\mathrm{\Delta }(x)=\mathrm{\Delta }+\delta \mathrm{\Delta }(x)`$ and $`\mathrm{\Phi }(x)=\varphi _\mathrm{m}+\delta \mathrm{\Phi }(x)`$, we can express the partition function $`Z`$ in terms of the fluctuation fields $`\delta \mathrm{\Delta }(x)`$ and $`\delta \mathrm{\Phi }(x)`$ as

$$Z=e^{\stackrel{~}{S}_0}𝒟\delta \mathrm{\Phi }^{}𝒟\delta \mathrm{\Phi }𝒟\delta \mathrm{\Delta }^{}𝒟\delta \mathrm{\Delta }e^{\delta \stackrel{~}{S}}.$$
(135)
Here the action $`\stackrel{~}{S}_0`$ does not involve the fluctuation fields, being given by

$`\stackrel{~}{S}_0`$ $`=`$ $`\beta {\displaystyle \underset{\text{p}}{}}\xi _\text{p}{\displaystyle \frac{\beta }{2}}{\displaystyle \underset{\text{q}}{}}\xi _{B\text{q}}\mathrm{Tr}\mathrm{ln}[\widehat{G}(x)^1]+\beta {\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}^2}{U_{\mathrm{eff}}}}`$ (136)
$`=`$ $`\beta {\displaystyle \underset{\text{p}}{}}\xi _\text{p}{\displaystyle \frac{\beta }{2}}{\displaystyle \underset{\text{q}}{}}\xi _{B\text{q}}{\displaystyle \underset{\text{p},\omega _m}{}}\mathrm{Tr}\mathrm{ln}[\widehat{G}(𝐩,i\omega _m)^1]+\beta {\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}^2}{U_{\mathrm{eff}}}}.`$ (137)
Here and in the following equations, $`G(x)`$ is given by Eq. (121), where $`\stackrel{~}{\mathrm{\Delta }}(x)`$ is replaced by the saddle point solution $`\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }g_\mathrm{r}\varphi _\mathrm{m}`$. In the last expression, $`\widehat{G}(𝐩,i\omega _m)`$ is given by Eq. (27) and the trace is over the particle-hole space. (The trace in the first line involves the integrals over both $`𝐫`$ and $`\tau `$.) In Eq. (135), $`\delta \stackrel{~}{S}`$ describes fluctuations around $`\mathrm{\Delta }`$ and $`\varphi _\mathrm{m}`$,

$`\delta \stackrel{~}{S}={\displaystyle _0^\beta }𝑑\tau {\displaystyle 𝑑𝐫\left[\frac{|\delta \mathrm{\Delta }(x)|^2}{U}\frac{1}{2}\delta \widehat{\mathrm{\Phi }}^{}(x)\widehat{D}^0(x)^1\delta \widehat{\mathrm{\Phi }}(x)\right]}+{\displaystyle \underset{n=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^n}{n}}\mathrm{Tr}[(\widehat{G}(x)\delta \widehat{\stackrel{~}{\mathrm{\Delta }}}(x))^n],`$ (138)
where $`\delta \widehat{\stackrel{~}{\mathrm{\Delta }}}`$ has the same form as Eq. (126), in which $`\mathrm{\Delta }(x)`$ and $`\mathrm{\Phi }(x)`$ are now replaced by the corresponding fluctuation variables $`\delta \mathrm{\Delta }(x)`$ and $`\delta \mathrm{\Phi }(x)`$, respectively. We note that since we expand $`\mathrm{\Delta }(x)`$ and $`\mathrm{\Phi }(x)`$ around the saddle point solutions ($`\mathrm{\Delta }`$, $`\varphi _\mathrm{m}`$), the first-order terms in $`\delta \mathrm{\Delta }(x)`$ and $`\delta \mathrm{\Phi }(x)`$ vanish.

We now employ the Gaussian fluctuation approximation to Eq. (138). This approximation corresponds to only retaining quadratic terms with respect to the fluctuation fields, giving

$$\delta \stackrel{~}{S}_0^\beta 𝑑\tau 𝑑𝐫\left[\frac{|\delta \mathrm{\Delta }(x)|^2}{U}\frac{1}{2}\delta \mathrm{\Phi }^{}(x)\widehat{D}^0(x)^1\delta \mathrm{\Phi }(x)\right]+\frac{1}{2}\mathrm{Tr}\left[\widehat{G}(x)\delta \widehat{\stackrel{~}{\mathrm{\Delta }}}(x)\widehat{G}(x)\delta \widehat{\stackrel{~}{\mathrm{\Delta }}}(x)\right].$$
(139)
Eq. (139) involves $`^2\stackrel{~}{S}/\stackrel{~}{\mathrm{\Delta }}^2`$ while neglects higher order terms $`^n\stackrel{~}{S}/\stackrel{~}{\mathrm{\Delta }}^n`$ ($`n3`$). This point is crucial when we calculate the equation of the total number of particles from the identity $`N=\mathrm{\Omega }/\mu `$: The thermodynamic potential $`\mathrm{\Omega }=T\mathrm{ln}Z`$ has the form

$`\mathrm{\Omega }`$ $`=`$ $`T\stackrel{~}{S}_0T\mathrm{ln}{\displaystyle 𝒟\delta \mathrm{\Phi }^{}𝒟\delta \mathrm{\Phi }𝒟\delta \mathrm{\Delta }^{}𝒟\delta \mathrm{\Delta }e^{\delta \stackrel{~}{S}}}`$ (140)
$``$ $`T\stackrel{~}{S}_0+\mathrm{\Lambda }(\mathrm{\Delta },\varphi _\mathrm{m},\mu ).`$ (141)
Since the last term in Eq. (141) is $`O(^2\stackrel{~}{S}/\stackrel{~}{\mathrm{\Delta }}^2`$), the derivative of this term with respect to $`\stackrel{~}{\mathrm{\Delta }}`$ gives a higher-order correction under the Gaussian fluctuation approximation. In addition, the $`\stackrel{~}{\mathrm{\Delta }}`$-derivative of $`\stackrel{~}{S}_0`$ vanishes due to the saddle point condition given in Eq. (131). As a result, although the saddle point solution ($`\mathrm{\Delta },\varphi _\mathrm{m})`$ depends on the chemical potential $`\mu `$, one does not need to keep the contribution in Eq. (84) in evaluating $`N=\mathrm{\Omega }/\mu `$, as we mentioned in Section III.

In order to carry out the functional integrals in $`\mathrm{\Lambda }`$ in Eq. (141), we divide the fluctuation fields into the real part and imaginary part as $`\delta \mathrm{\Delta }(x)=\delta \mathrm{\Delta }_R(x)i\delta \mathrm{\Delta }_I(x)`$ and $`\delta \mathrm{\Phi }(x)=\delta \mathrm{\Phi }_R(x)i\delta \mathrm{\Phi }_I(x)`$. Then we can write $`\delta \widehat{\stackrel{~}{\mathrm{\Delta }}}(x)=A(x)\tau _1+B(x)\tau _2`$, where $`A(x)`$ and $`B(x)`$ are real quantities. When we transform $`\delta \stackrel{~}{S}`$ in Eq. (139) from the $`x(=(𝐫,\tau ))`$-representation to $`q(=(𝐪,\nu _n))`$-representation, we find

$`\delta \stackrel{~}{S}`$ $`=`$ $`{\displaystyle \frac{2}{\beta U}}{\displaystyle \underset{q}{}^{}}\widehat{X}_q^{}\left[1+{\displaystyle \frac{U}{2}}\widehat{\mathrm{\Pi }}^0(q)\right]\widehat{X}_q{\displaystyle \frac{g_\mathrm{r}}{\beta }}{\displaystyle \underset{q}{}^{}}\left[\widehat{X}_q^{}\widehat{\mathrm{\Pi }}^0(q)\widehat{Y}_q+\widehat{Y}_q^{}\mathrm{\Pi }^0(q)\widehat{X}_q\right]`$ (142)
$`+`$ $`{\displaystyle \frac{g_\mathrm{r}^2}{\beta }}{\displaystyle \underset{q}{}^{}}\widehat{Y}_q^{}\widehat{\mathrm{\Pi }}^0\widehat{Y}_q{\displaystyle \frac{2}{\beta }}{\displaystyle \underset{q}{}^{}}\widehat{Y}_q^{}\widehat{W}\widehat{D}^0(q)^1\widehat{W}^1\widehat{Y}_q.`$ (143)
Here $`\widehat{X}_q^{}=(\delta \mathrm{\Delta }_R(q),\delta \mathrm{\Delta }_I(q))`$, $`\widehat{Y}_q^{}=(\delta \mathrm{\Phi }_R(q),\delta \mathrm{\Phi }_I(q))`$. (Note that $`\delta \mathrm{\Delta }_{R,I}^{}(q)=\delta \mathrm{\Delta }_{R,I}(q)`$ and $`\delta \mathrm{\Phi }_{R,I}^{}(q)=\delta \mathrm{\Phi }_{R,I}(q)`$.) The prime in the $`q`$-summation means that one does not take $`q`$ when one takes $`q`$, in order to avoid double counting. $`\widehat{\mathrm{\Pi }}^0(q)`$, $`\widehat{W}`$ and $`\widehat{D}^0(q)`$ are, respectively, defined in Eqs. (63), (75) and (65). Executing the functional integral with respect to $`\widehat{X}_q`$, we obtain

$`\delta \stackrel{~}{S}`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{q}{}^{}}\mathrm{Tr}\mathrm{ln}\left[1+{\displaystyle \frac{U}{2}}\widehat{\mathrm{\Pi }}^0(q)\right]{\displaystyle \frac{g_\mathrm{r}^2}{\beta }}{\displaystyle \underset{q}{}^{}}\widehat{Y}_q^{}\widehat{\mathrm{\Pi }}^0(q)\left[1+{\displaystyle \frac{U}{2}}\widehat{\mathrm{\Pi }}^0(q)\right]^1\widehat{\mathrm{\Pi }}^0(q)\widehat{Y}_q`$ (144)
$`+`$ $`{\displaystyle \frac{g_\mathrm{r}^2}{\beta }}{\displaystyle \underset{q}{}^{}}\widehat{Y}_q^{}\widehat{\mathrm{\Pi }}^0\widehat{Y}_q{\displaystyle \frac{2}{\beta }}{\displaystyle \underset{q}{}^{}}\widehat{Y}_q^{}\widehat{W}D^0(q)^1\widehat{W}^1\widehat{Y}_q`$ (145)
$`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{q}{}}\mathrm{Tr}\left[1+{\displaystyle \frac{U}{2}}\widehat{\mathrm{\Pi }}^0(q)\right]{\displaystyle \frac{2}{\beta }}{\displaystyle \underset{q}{}^{}}\widehat{Y}_q^{}\left[\widehat{W}\widehat{D}^0(q)^1\widehat{W}^1{\displaystyle \frac{g_\mathrm{r}^2}{2}}\widehat{\mathrm{\Pi }}^0(q)\left[1+{\displaystyle \frac{U}{2}}\widehat{\mathrm{\Pi }}^0(q)\right]^1\right]\widehat{Y}_q.`$ (146)
In obtaining Eq. (LABEL:eq.2.54), we have used the formula

$$𝒟\widehat{X}_q^{}𝒟\widehat{X}_qe^{_q[\widehat{X}_q^{}\widehat{M}_q\widehat{X}_q+B_q^{}\widehat{X}_q+\widehat{X}_q^{}C_q]}=\frac{1}{\mathrm{det}[\widehat{M}]}e^{_qB_q^{}M_q^1C_q}.$$
(148)
Carrying out the functional integral with respect to $`\widehat{Y}_q`$ using this formula, the second term in Eq. (141) is given by

$`\mathrm{\Lambda }(\mathrm{\Delta },\varphi _\mathrm{m},\mu )`$ $`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{q}{}}\mathrm{Tr}\mathrm{ln}\left[1+{\displaystyle \frac{U}{2}}\widehat{\mathrm{\Pi }}^0(q)\right]`$ (149)
$`+`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{q}{}}\mathrm{Tr}\mathrm{ln}\left[\widehat{W}\widehat{D}^0(q)^1\widehat{W}^1+{\displaystyle \frac{g_\mathrm{r}^2}{2}}\widehat{\mathrm{\Pi }}^0(q)\left[1+{\displaystyle \frac{U}{2}}\widehat{\mathrm{\Pi }}^0(q)\right]^1\right]`$ (150)
$`=`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{q}{}}\mathrm{Tr}\mathrm{ln}[\widehat{D}^0(q)^1]+{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{q}{}}\mathrm{Tr}\mathrm{ln}\left[1+{\displaystyle \frac{1}{4}}[Ug_\mathrm{r}^2\widehat{D}^0(q)^1]\widehat{\mathrm{\Xi }}^0(q)\right].`$ (151)
The matrix correlation function $`\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)`$ is defined in Eq. (72).

Substituting Eqs. (137) and (151) into Eq. (141), the thermodynamic potential $`\mathrm{\Omega }`$ is given by

$`\mathrm{\Omega }`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}^2}{U_{\mathrm{eff}}}}+{\displaystyle \underset{\text{p}}{}}\left[\xi _\text{p}{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\text{p},\omega _m}{}}\mathrm{Tr}\mathrm{ln}[\widehat{G}(𝐩,i\omega _m)^1]\right]`$ (152)
$`+`$ $`\left[{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{q}{}}\mathrm{Tr}\mathrm{ln}[\widehat{D}^0(q)^1]{\displaystyle \frac{1}{2}}{\displaystyle \underset{\text{q}}{}}\xi _{B\text{q}}\right]`$ (153)
$`+`$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{q}{}}\mathrm{Tr}\mathrm{ln}\left[1+{\displaystyle \frac{1}{4}}(Ug_\mathrm{r}^2\widehat{D}^0(q)^1)\widehat{\mathrm{\Xi }}^0(q)\right].`$ (154)
In turn, the total number of particles $`N=\mathrm{\Omega }/\mu `$ is given by

$`N`$ $`=`$ $`2\varphi _\mathrm{m}^2+{\displaystyle \underset{\text{p}}{}}\left\{1+{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\omega _m}{}}\mathrm{Tr}\left[\widehat{G}(p){\displaystyle \frac{\widehat{G}(p)^1}{\mu }}\right]\right\}{\displaystyle \underset{\text{q}}{}}\left\{1+{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\nu _n}{}}{\displaystyle \frac{}{\mu }}\mathrm{Tr}\left[\widehat{D}^0(q){\displaystyle \frac{\widehat{D}^0(q)^1}{\mu }}\right]\right\}`$ (155)
$``$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}{\displaystyle \frac{}{\mu }}\mathrm{Tr}\mathrm{ln}\left[1+{\displaystyle \frac{1}{4}}(Ug_\mathrm{r}^2\widehat{D}^0(q)^1)\widehat{\mathrm{\Xi }}^0(q)\right]`$ (156)
$`=`$ $`2\varphi _\mathrm{m}^2+{\displaystyle \underset{\text{p}}{}}\left[1{\displaystyle \frac{\xi _\text{p}}{E_\text{p}}}\mathrm{tanh}{\displaystyle \frac{\beta }{2}}E_\text{p}\right]+2{\displaystyle \underset{\text{q}}{}}n_\mathrm{B}(\xi _{B\text{q}})`$ (157)
$``$ $`{\displaystyle \frac{1}{2\beta }}{\displaystyle \underset{\text{q},\nu _n}{}}{\displaystyle \frac{}{\mu }}\mathrm{Tr}\mathrm{ln}\left\{1+{\displaystyle \frac{1}{4}}\left[Ug_\mathrm{r}^2\widehat{D}^0(𝐪,i\nu _n)\right]\widehat{\mathrm{\Xi }}^0(𝐪,i\nu _n)\right\}.`$ (158)
This result is identical to that given in Eq. (86). The reproduction of the results in Eqs. (37) and (86) \[see Eqs. (134) and (158)\] explicitly proves that the strong-coupling theory presented in Section III in terms of summing up diagrams corresponds to keeping Gaussian fluctuations around the (mean-field) order parameters $`\mathrm{\Delta }`$ and $`\varphi _\mathrm{m}`$.

## V Goldstone modes in the BCS-BEC crossover Region

### A Correlation function and $`b`$-boson Green’s function in the HF-RPA

The Goldstone mode in fermion superfluidity (Anderson-Bogoliubov) is a collective phase oscillation (phason) of Cooper-pairs, and thus it appears as a pole in the phase correlation function $`\mathrm{\Pi }_{22}(𝐪,i\nu _n\omega +i\delta )`$. In the present coupled fermion-boson model, we also expect a Bogoliubov phonon mode associated with the BEC of $`b`$-molecules characterized by the Bose order parameter $`\varphi _\mathrm{m}=b_{\text{q}=0}`$. This mode appears in the excitation spectrum of the $`b`$-boson Green’s function $`D(𝐪,\tau )=T_\tau \{b_\text{q}(\tau )b_\text{q}^{}(0)\}`$. However, since the Cooper-pair amplitude $`\mathrm{\Delta }`$ and the $`b`$-boson order parameter $`\varphi _\mathrm{m}`$ are coupled with each other via the Feshbach resonance \[see Eq. (8)\], these two Goldstone modes are strongly hybridized.

In calculating the Goldstone mode, we have to be careful use a consistent approximation for the self-energy and vertex correction. In this regard, apart from the chemical potential, we recall that the gap equation in Eq. (37) is obtained from the Hartree-Fock Green’s function in Eq. (27). Thus we should employ the HF-RPA formalism for the correlation functions $`\widehat{\mathrm{\Pi }}=\{\mathrm{\Pi }_{ij}\}`$ and $`b`$-boson Green’s function. Using the HF-RPA guarantees the gapless behavior of the Goldstone mode.

In calculating the correlation function $`\widehat{\mathrm{\Pi }}`$, it is convenient to first consider $`\widehat{\mathrm{\Xi }}=2\widehat{W}^1\widehat{\mathrm{\Pi }}\widehat{W}`$ \[see Eq. (72)\]. Within the HF-RPA, $`\widehat{\mathrm{\Xi }}`$ is described by the sum of the diagrams in Fig. 5(a). The summation gives

$`\widehat{\mathrm{\Xi }}`$ $`=`$ $`\widehat{\mathrm{\Xi }}_U\left[1{\displaystyle \frac{1}{4}}g_\mathrm{r}^2\widehat{D}^0\widehat{\mathrm{\Xi }}_U\right]^1=\widehat{\mathrm{\Xi }}^0\left[1+{\displaystyle \frac{1}{4}}\left[Ug_\mathrm{r}^2\widehat{D}^0\right]\widehat{\mathrm{\Xi }}^0\right]^1,`$ (159)
where $`\widehat{\mathrm{\Xi }}_U`$ involves the effect of non-resonant interaction $`U`$ using the HF-RPA, as shown in Fig. 5(c) and given by Eq. (76). The correlation function $`\widehat{\mathrm{\Pi }}`$ is then obtained from the inverse relation $`\widehat{\mathrm{\Pi }}=\frac{1}{2}\widehat{W}\widehat{\mathrm{\Xi }}\widehat{W}^1`$, namely

$$\widehat{\mathrm{\Pi }}=\widehat{\mathrm{\Pi }}^0\left[1+\frac{1}{2}\left[Ug_\mathrm{r}^2\widehat{W}\widehat{D}^0\widehat{W}^1\right]\widehat{\mathrm{\Pi }}^0\right]^1.$$
(160)
More explicitly, the amplitude and phase correlation functions $`\mathrm{\Pi }_{11}`$ and $`\mathrm{\Pi }_{22}`$ are given by

$`\{\begin{array}{c}\mathrm{\Pi }_{11}(𝐪,i\nu _n)={\displaystyle \frac{\stackrel{~}{\mathrm{\Pi }}_{11}}{1+{\displaystyle \frac{V_1}{2}}\stackrel{~}{\mathrm{\Pi }}_{11}+{\displaystyle \frac{V_2^2}{4}}{\displaystyle \frac{\mathrm{\Pi }_{11}^0\mathrm{\Pi }_{22}^0\mathrm{\Pi }_{12}^0\mathrm{\Pi }_{21}^0}{1+\frac{V_1}{2}\mathrm{\Pi }_{22}^0}}V_2{\displaystyle \frac{\mathrm{\Pi }_{12}^0}{1+\frac{V_1}{2}\mathrm{\Pi }_{22}^0}}}},\hfill \\ \mathrm{\Pi }_{22}(𝐪,i\nu _n)={\displaystyle \frac{\stackrel{~}{\mathrm{\Pi }}_{22}}{1+{\displaystyle \frac{V_1}{2}}\stackrel{~}{\mathrm{\Pi }}_{22}+{\displaystyle \frac{V_2^2}{4}}{\displaystyle \frac{\mathrm{\Pi }_{11}^0\mathrm{\Pi }_{22}^0\mathrm{\Pi }_{12}^0\mathrm{\Pi }_{21}^0}{1+\frac{V_1}{2}\mathrm{\Pi }_{11}^0}}V_2{\displaystyle \frac{\mathrm{\Pi }_{12}^0}{1+\frac{V_1}{2}\mathrm{\Pi }_{11}^0}}}},\hfill \end{array}`$ (163)
where

$`\{\begin{array}{c}\stackrel{~}{\mathrm{\Pi }}_{11}\mathrm{\Pi }_{11}^0+\mathrm{\Pi }_{12}^0{\displaystyle \frac{\frac{V_1}{2}}{1+\frac{V_1}{2}\mathrm{\Pi }_{22}^0}}\mathrm{\Pi }_{21}^0,\hfill \\ \stackrel{~}{\mathrm{\Pi }}_{22}\mathrm{\Pi }_{22}^0+\mathrm{\Pi }_{21}^0{\displaystyle \frac{\frac{V_1}{2}}{1+\frac{V_1}{2}\mathrm{\Pi }_{11}^0}}\mathrm{\Pi }_{12}^0,\hfill \end{array}`$ (166)
$`\{\begin{array}{c}V_1(𝐪,i\nu _n)Re[U_{\mathrm{eff}}(𝐪,i\nu _n)]=U+g_\mathrm{r}^2{\displaystyle \frac{\xi _{B\text{q}}}{\nu _n^2+\xi _{B\text{q}}^2}},\hfill \\ V_2(𝐪,i\nu _n)Im[U_{\mathrm{eff}}(𝐪,i\nu _n)]=g_\mathrm{r}^2{\displaystyle \frac{\nu _n}{\nu _n^2+\xi _{B\text{q}}^2}}.\hfill \end{array}`$ (169)
Here the lowest-order non-interacting correlation functions $`\mathrm{\Pi }_{ij}^0`$ are defined in Eqs. (43)-(53). $`V_1`$ describes amplitude-amplitude and phase-phase interactions, while $`V_2`$ is an amplitude-phase coupling mediated by $`b`$-molecules. Since both $`V_2`$ and $`\mathrm{\Pi }_{12}^0`$ vanish at $`𝐪=i\nu _n=0`$, the denominator of $`\mathrm{\Pi }_{22}`$ in Eq. (163) reduces to $`1+\frac{U_{\mathrm{eff}}}{2}\mathrm{\Pi }_{22}^0(0,0)`$. This vanishes when the gap equation in Eq. (54) is satisfied below $`T_\mathrm{c}`$. This proves that the collective phase oscillation is always gapless at $`𝐪=0`$, a requirement of the Anderson-Bogoliubov (Goldstone) mode.

The $`b`$-boson Green’s function consistent with treating $`\widehat{\mathrm{\Pi }}`$ in the HF-RPA is shown in terms of diagrams in Fig. 5(b). The result is just the same as the renormalized $`b`$-boson Green’s function $`\widehat{D}`$ given in Eq. (97). As discussed at the end of Section III, the excitation spectrum of $`\widehat{D}`$ is gapless at $`𝐪=0`$ in the superfluid phase below $`T_\mathrm{c}`$, which is again consistent with the expected gapless Bogoliubov phonon mode.

### B Goldstone mode

As discussed in the preceding subsection, the excitation spectrum of the renormalized $`b`$-boson is determined by

$$\mathrm{det}[\widehat{D}(𝐪,i\nu _n\omega +i\delta )^1]=0.$$
(170)
On the other hand, the collective phase oscillation and amplitude oscillation are obtained from the poles of $`\widehat{\mathrm{\Pi }}`$ as given by Eq. (160)

$`0`$ $`=`$ $`\mathrm{det}\left[1+{\displaystyle \frac{1}{2}}\left[Ug_\mathrm{r}^2\widehat{W}\widehat{D}^0\widehat{W}^1\right]\widehat{\mathrm{\Pi }}^0\right]_{i\nu _n\omega _+=\omega +i\delta }`$ (171)
$`=`$ $`\mathrm{det}[\widehat{D}^0(𝐪,\omega _+)]\mathrm{det}\left[1+{\displaystyle \frac{U}{2}}\widehat{\mathrm{\Pi }}^0(𝐪,\omega _+)\right]\mathrm{det}\left[\widehat{D}^0(𝐪,\omega _+)^1{\displaystyle \frac{g_\mathrm{r}^2}{4}}\widehat{\mathrm{\Xi }}_U(𝐪,\omega _+)\right]`$ (172)
$`=`$ $`{\displaystyle \frac{\mathrm{det}\left[1+\frac{U}{2}\widehat{\mathrm{\Pi }}^0(𝐪,\omega _+)\right]}{\xi _{B\text{q}}^2\omega _+^2}}\mathrm{det}[\widehat{D}(𝐪,\omega _+)^1].`$ (173)
Comparing Eqs. (170) with (173), we see that the correlation functions $`\mathrm{\Pi }_{ij}`$ ($`i,j=1,2`$) and the $`b`$-boson Green’s function $`\widehat{D}`$ have the identical poles (unless $`\mathrm{det}[1+\frac{U}{2}\widehat{\mathrm{\Pi }}^0(𝐪,\omega _+)]`$ has zeros). This equivalence is due to hybridizing effects from the coupling between the amplitude fluctuations, phase fluctuations and the $`b`$-bosons with the amplitude-phase coupling $`\mathrm{\Pi }_{12}^0`$ and the Feshbach resonance coupling $`g_\mathrm{r}`$. Thus, in principle, we can consider either one of the correlation functions $`\mathrm{\Pi }_{ij}`$ or the $`b`$-boson Green’s function $`\widehat{D}`$ when we want to study the Goldstone mode. We work with the phase correlation function $`\mathrm{\Pi }_{22}`$ in Eq. (163), in which case the dispersion relation of the collective mode is given by

$$1+\frac{V_1}{2}\stackrel{~}{\mathrm{\Pi }}_{22}+\frac{V_2^2}{4}\frac{\mathrm{\Pi }_{11}^0\mathrm{\Pi }_{22}^0\mathrm{\Pi }_{12}^0\mathrm{\Pi }_{21}^0}{1+\frac{V_1}{2}\mathrm{\Pi }_{11}^0}V_2\frac{\mathrm{\Pi }_{12}^0}{1+\frac{V_1}{2}\mathrm{\Pi }_{11}^0}=0.$$
(174)
Landau-damping associated with thermally excited fermion quasi-particles leads to the imaginary part of collective oscillations and, apart from $`T=0`$, we have to look for a complex solution to Eq. (174). This requires a complicated analysis. In this paper, for simplicity, we only consider the real part of this equation,

$$Re\left[1+\frac{V_1}{2}\stackrel{~}{\mathrm{\Pi }}_{22}+\frac{V_2^2}{4}\frac{\mathrm{\Pi }_{11}^0\mathrm{\Pi }_{22}^0\mathrm{\Pi }_{12}^0\mathrm{\Pi }_{21}^0}{1+\frac{V_1}{2}\mathrm{\Pi }_{11}^0}V_2\frac{\mathrm{\Pi }_{12}^0}{1+\frac{V_1}{2}\mathrm{\Pi }_{11}^0}\right]=0.$$
(175)
From this equation, we can obtain real frequencies as approximate solutions. In order to check the validity of this prescription, we also solve another approximate equation for the mode energy obtained from the renormalized $`b`$-boson Green’s function,

$$Re[\widehat{D}_{11}^0(𝐪,i\nu \omega _+)^1]=0.$$
(176)
As for the damping of the collective mode, we investigate this effect by examining the width of the peak of the collective mode in the spectrum of the correlation functions given by $`Im[\mathrm{\Pi }_{11}(𝐪,\omega _+)]`$ and $`Im[\mathrm{\Pi }_{22}(𝐪,\omega _+)]`$, as well as the $`b`$-boson excitation spectrum $`Im[\widehat{D}_{11}(𝐪,i\nu _n\omega _+)]`$. The structure function

$$S_{jj}(𝐪,\omega )\frac{1}{\pi }[n_\mathrm{B}(\omega )+1]Im[\widehat{\mathrm{\Pi }}_{jj}(𝐪,i\nu \omega _+)](j=1,2).$$
(177)
is more convenient than $`Im[\widehat{\mathrm{\Pi }}_{jj}]`$ in studying collective behavior, because a diffusive mode spectrum which arises appears as a central peak at $`\omega =0`$ in $`S_{jj}`$, which can be easily distinguished from collective modes appearing at a finite frequency. In contrast, this diffusive mode shows up as a peak at a finite frequency $`\omega _d`$ in $`Im[\widehat{\mathrm{\Pi }}_{jj}]`$, which exhibits a structure of the kind $`\omega /(\omega ^2+\omega _d^2)`$.

### C Goldstone mode at zero temperature

At zero temperature, since Landau-damping is absent for modes below the excitation gap $`2\stackrel{~}{\mathrm{\Delta }}`$, we can deal directly with Eq. (174) in determining the Goldstone mode for frequencies $`\omega `$ below the excitation gap. In the long wavelength or phonon limit, we take $`\omega =v_\varphi q`$ and expand Eq. (174) up to the quadratic order in terms of $`\omega `$ and $`q`$. After some calculation, the velocity of the Goldstone mode $`v_\varphi `$ is given by

$`v_\varphi ={\displaystyle \frac{1}{\sqrt{2m}}}{\displaystyle \frac{\sqrt{B+{\displaystyle \frac{1}{U_{\mathrm{eff}}^2}}{\displaystyle \frac{g_\mathrm{r}^2}{(2\nu 2\mu )^2}}}}{\sqrt{A+\eta ^2{\displaystyle \frac{\frac{U_{\mathrm{eff}}}{2}}{1+\frac{U_{\mathrm{eff}}}{2}\mathrm{\Pi }_{11}^0(0,0)}}+{\displaystyle \frac{\frac{2}{U_{\mathrm{eff}}}\frac{g_\mathrm{r}^2}{(2\nu 2\mu )^2}}{1+\frac{U_{\mathrm{eff}}}{2}\mathrm{\Pi }_{11}^0(0,0)}}\left[{\displaystyle \frac{1+\frac{U}{2}\mathrm{\Pi }_{11}^0(0,0)}{[2\nu 2\mu ]U_{\mathrm{eff}}}}\eta \right]}}}.`$ (178)
The factors $`A`$ and $`B`$ are obtained from the expansion of $`\mathrm{\Pi }_{22}^0`$ as $`\mathrm{\Pi }_{22}^0(𝐪,i\nu _n)=\mathrm{\Pi }_{22}^0(0,0)+A\nu _n^2+Bq^2/2m`$, with the explicit expressions

$`\{\begin{array}{c}A={\displaystyle \frac{1}{4}}{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{1}{E_\text{p}^3}},\hfill \\ B={\displaystyle \underset{\text{p}}{}}\left[\epsilon _\text{p}{\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}^2}{2E_\text{p}^5}}+{\displaystyle \frac{\xi _\text{p}}{4E_\text{p}^3}}\right].\hfill \end{array}`$ (181)
Finally, the factor $`\eta `$ is related to the amplitude-phase coupling $`\mathrm{\Pi }_{12}^0`$, namely

$$\eta \frac{\mathrm{\Pi }_{12}^0(𝐪,i\nu _n)}{\nu _n}|_{T=\nu _n=0}=\frac{1}{2}\underset{\text{p}}{}\frac{\xi _\text{p}}{E_\text{p}^3}.$$
(182)
The second term in the denominator in Eq. (178) describes the effect of amplitude-phase coupling (second order in $`\eta `$), while the second term in the numerator and the third term in the denominator involve the effect of the Feshbach resonance coupling (second order in $`g_\mathrm{r}`$).

#### 1 BCS regime: $`2\nu 2\epsilon _\mathrm{F}`$

In the BCS limit ($`2\nu 2\mu `$), the terms involving the factor $`1/(2\nu 2\mu )`$ can be neglected in Eq. (178). In addition, since the region near the Fermi surface dominates just as in ordinary weak-coupling BCS theory, we may take

$$\underset{\text{p}}{}g(\xi _\text{p})N(\mu (0))_{\mathrm{}}^{\mathrm{}}𝑑\xi g(\xi ),$$
(183)
where $`N(\mu (T=0))`$ is the fermion density of states (DOS) at the Fermi surface. In this approximation, $`\eta `$ coming from the amplitude-phase correlation function $`\mathrm{\Pi }_{12}^0`$ vanishes, and Eq. (178) reduces to

$$v_\varphi =\frac{1}{\sqrt{2m}}\sqrt{\frac{B}{A}}=\frac{1}{\sqrt{3}}\overline{v}_\mathrm{F},$$
(184)
where $`\overline{v}_\mathrm{F}\sqrt{2m\mu (0)}`$. In evaluating $`B`$, we have approximated $`\epsilon _\text{p}`$ appearing in Eq. (181) by the Fermi energy $`\mu (0)`$. Eq. (184) is the well-known velocity of the Anderson-Bogoliubov phonon in weak-coupling BCS superfluidity.

The same result as in Eq. (184) is given by the pole of the phase correlation function $`\mathrm{\Pi }_{22}`$ in the BCS limit. $`\mathrm{\Pi }_{22}`$ in this limit is given by \[taking $`\mathrm{\Pi }_{12}^0=0`$, $`V_1=U`$ and $`V_2=0`$ in Eq. (163)\]

$$\mathrm{\Pi }_{22}(𝐪,\omega _+)=\frac{\mathrm{\Pi }_{22}^0(𝐪,\omega _+)}{1+\frac{U}{2}\mathrm{\Pi }_{22}^0(𝐪,\omega _+)}=\frac{\frac{2}{U}\mathrm{\Pi }_{22}^0(𝐪,\omega _+)}{\mathrm{\Pi }_{22}^0(𝐪,\omega _+)\mathrm{\Pi }_{22}^0(0,0)}.$$
(185)
Using $`\mathrm{\Pi }_{22}^0(𝐪,\omega _+)=\mathrm{\Pi }_{22}^0(0,0)A\omega ^2+Bq^2/2m`$, we easily find that the poles of Eq. (185) are phonons with the velocity given by Eq. (184). This shows that the Goldstone mode in the BCS limit is a pure collective phase oscillation of the Cooper-pair order parameter $`\mathrm{\Delta }`$.

The amplitude-phase coupling $`\mathrm{\Pi }_{12}^0`$ in Eq. (20) vanishes in the BCS limit when we use Eq. (183) and retain up to $`O(q/p_\mathrm{F})`$, where $`p_\mathrm{F}`$ is the Fermi momentum. In this case, the amplitude mode is decoupled from Eq. (174). This mode can then be obtained from the amplitude correlation function working within the HF-RPA,

$$\mathrm{\Pi }_{11}(𝐪,\omega _+)=\frac{\mathrm{\Pi }_{11}^0(𝐪,\omega _+)}{1+\frac{U}{2}\mathrm{\Pi }_{11}^0(𝐪,\omega _+)}=\frac{\frac{2}{U}\mathrm{\Pi }_{11}^0(𝐪,\omega _+)}{\mathrm{\Pi }_{11}^0(𝐪,\omega _+)\mathrm{\Pi }_{22}^0(0,0)}.$$
(186)
In particular, at $`𝐪=0`$, Eq. (186) reduces to

$`\mathrm{\Pi }_{11}(0,\omega _+)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Pi }_{11}^0(0,\omega _+)}{UN(\mu (0))}}{\displaystyle \frac{1}{\mathrm{tan}^1\frac{\omega }{\sqrt{4\mathrm{\Delta }^2\omega ^2}}}}{\displaystyle \frac{\omega }{\sqrt{4\mathrm{\Delta }^2\omega ^2}}}(\omega 2\mathrm{\Delta })`$ (187)
$``$ $`{\displaystyle \frac{2\mathrm{\Pi }_{11}^0(0,2\mathrm{\Delta })}{\pi UN(\mu (0))}}{\displaystyle \frac{\mathrm{\Delta }}{\sqrt{4\mathrm{\Delta }^2\omega ^2}}}(\omega 2\mathrm{\Delta }).`$ (188)
Thus, although $`\omega =2\mathrm{\Delta }`$ is seen to be a branch cut rather than a pole, a strong peak is expected at the excitation gap $`\omega =2\mathrm{\Delta }`$ in the spectrum of the amplitude correlation function $`Im[\widehat{\mathrm{\Pi }}_{11}]`$ and hence $`S_{11}(𝐪,\omega )`$.

#### 2 BEC regime: $`\nu <0`$

In the BEC regime ($`\nu <0`$), since the chemical potential $`\mu `$ approaches $`\nu `$ as $`\nu `$ decreases, the effective interaction $`U_{\mathrm{eff}}`$ defined in Eq. (11) is dominated by the Feshbach contribution, $`g_\mathrm{r}^2/(2\nu 2\mu )U`$. On the other hand, the correlation functions $`\mathrm{\Pi }_{ij}^0`$ become less dominant. Since $`|\mu ||\nu |\stackrel{~}{\mathrm{\Delta }}`$, the energy gap is also less important in the excitation spectrum of fermion quasi-particles, so that we can approximate $`E_\text{p}=\sqrt{(\epsilon _\text{p}\mu )^2+\stackrel{~}{\mathrm{\Delta }}^2}\epsilon _\text{p}+|\nu |`$. In this limiting case, Eq. (178) reduces to

$`v_\varphi `$ $`=`$ $`{\displaystyle \frac{g_\mathrm{r}\stackrel{~}{\mathrm{\Delta }}}{\sqrt{8m}}}\left[{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{1}{(\epsilon _\text{p}+|\nu |)^3}}\right]^{1/2}`$ (189)
$`=`$ $`{\displaystyle \frac{g_\mathrm{r}^{}\stackrel{~}{\mathrm{\Delta }}}{\sqrt{8m}}}\left[{\displaystyle \frac{3\pi }{32|\nu \epsilon _\mathrm{F}|^{3/2}}}\right]^{1/2},`$ (190)
where we have rescaled the Feshbach coupling as $`g_\mathrm{r}\sqrt{N}g_\mathrm{r}^{}`$ in the last expression. In evaluating the $`𝐩`$-summation in Eq. (190), we have taken into account the correct energy dependence of the DOS ($`\sqrt{\epsilon _\text{p}}`$) in contrast to the approximation in Eq. (183). Since $`\varphi _\mathrm{m}`$ is the dominant contribution to $`\stackrel{~}{\mathrm{\Delta }}`$ in the BEC regime, Eq. (190) is proportional to $`\sqrt{N_\mathrm{B}^\mathrm{c}}`$ (where $`N_\mathrm{B}^\mathrm{c}=\varphi _\mathrm{m}^2`$ is the number of condensed $`b`$-molecules). This dependence on $`N_\mathrm{B}^\mathrm{c}`$ is characteristic of the Bogoliubov phonon mode in a Bose-condensed gas. Thus, Eq. (190) may be regarded as the velocity of the Bogoliubov phonon associated with a condensate of $`b`$-molecules. We also note that the expression in Eq. (190) gives as $`1/|\nu |^{3/4}`$ and this approaches zero when $`\nu \mathrm{}`$. In this limit, the superfluidity is described by a BEC of a free gas with $`N/2`$ bosons having the particle-like excitation spectrum $`\omega =q^2/2M`$, with no linear (or phonon) component.

The Bogoliubov phonon associated with a Bose condensate of paired fermions has also been discussed in strong-coupling superconductivity, as well as for excitons in optically-excited semiconductors. Eq. (178) reproduces the Bogoliubov phonon velocity in strong-coupling superconductivity. In this case, the Feshbach coupling $`g_\mathrm{r}`$ and the $`b`$-boson are absent, while the non-resonant interaction $`U`$ is taken to be strong. Then Eq. (178) reduces to

$`v_\varphi `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2m}}}B^{1/2}\left[A+{\displaystyle \frac{U}{2}}{\displaystyle \frac{\eta ^2}{1+\frac{U}{2}\mathrm{\Pi }_{11}^0(0,0)}}\right]^{1/2}`$ (191)
$``$ $`{\displaystyle \frac{1}{\sqrt{2m}}}\left[{\displaystyle \frac{B(\mathrm{\Delta }=0)}{\eta ^2(\mathrm{\Delta }=0)}}{\displaystyle \underset{\text{p}}{}}{\displaystyle \frac{\mathrm{\Delta }^2}{(\epsilon _\text{p}+|\mu |^2)^3}}\right]^{1/2}`$ (192)
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{\sqrt{8m|\mu |}}},`$ (193)
where we have taken the strong-coupling limit ($`U\mathrm{}`$), with $`\mu `$ large and negative. The last expression can be shown to be equivalent to the Bogoliubov phonon velocity in strong-coupling superconductivity, as discussed in Ref.. Using the strong-coupling expressions, $`\mathrm{\Delta }=(16/3\pi )^{1/2}\epsilon _\mathrm{F}/\sqrt{p_\mathrm{F}a_s}`$ and $`\mu =1/2ma_s^2`$, where $`a_s`$ is an $`s`$-wave scattering length of the pairing interaction between fermions and $`p_\mathrm{F}`$ is the Fermi momentum, Eq. (193) reduces to a more familiar form,

$$v_\varphi =\left[\frac{n_\mathrm{B}U_\mathrm{B}}{M}\right]^{1/2}.$$
(194)
Here $`n_\mathrm{B}N/2V`$ is the number density of Cooper-pair bosons (where $`V`$ is the volume of a system) and $`U_\mathrm{B}4\pi a_\mathrm{B}/M`$ ($`M=2m`$) is an effective interaction between Cooper-pairs with an $`s`$-wave scattering length $`a_\mathrm{B}2a_s`$.

## VI BCS-BEC crossover in the superfluid phase

In this Section, we present numerical results for the various thermodynamic parameters and correlation functions discussed in earlier Sections, as one passes through the BCS-BEC crossover region. We give results as a function of the temperature in the superfluid phase ($`T<T_\mathrm{c}`$) for different values of the $`b`$-molecule threshold $`2\nu `$. We take the Fermi energy $`\epsilon _\mathrm{F}`$, Fermi momentum $`p_\mathrm{F}`$ and Fermi velocity $`v_\mathrm{F}`$ as the units of energy, momentum and velocity, respectively, where $`\epsilon _\mathrm{F}`$, $`p_\mathrm{F}`$, and $`v_\mathrm{F}`$ are all in the absence of coupling to $`b`$-bosons. As the unit for the number for particles, we take the total number of Fermi atoms $`N`$. Since the interactions $`U`$ and $`g_\mathrm{r}`$ always appear with $`N`$ as $`UN`$ and $`g_\mathrm{r}^2N`$, we rescale them as $`UNU`$ and $`g_\mathrm{r}\sqrt{N}g_\mathrm{r}`$. As for the energy cutoff which is necessary in the gap equation in Eq. (37) and the correlation functions $`\mathrm{\Pi }_{ij}^0`$ in Eqs. (43)-(53), we employ the Gaussian cutoff $`e^{(\epsilon _\text{p}/\omega _c)^2}`$ with $`\omega _c=2\epsilon _\mathrm{F}`$, as in our previous work for $`T_\mathrm{c}`$. Since our “strong-coupling” theory is based on perturbative expansions with respect to $`U`$ and $`g_\mathrm{r}`$, these coupling terms are assumed to be weak perturbations. For this reason, we take $`g_\mathrm{r}=0.6\epsilon _\mathrm{F}`$ and $`U=0.3\epsilon _\mathrm{F}`$ in all our numerical calculations. In Ref., we have discussed the BCS-BEC crossover at $`T_\mathrm{c}`$ for the case of a very broad Feshbach resonance ($`g_\mathrm{r}\epsilon _\mathrm{F}`$).

### A Temperature dependence of the order parameter and chemical potential

The self-consistent solutions ($`\stackrel{~}{\mathrm{\Delta }},\mu `$) of the coupled equations (37) and (86) are shown in Figs. 6 and 7. In Fig. 6(a), we find that the temperature dependence of the order parameter $`\stackrel{~}{\mathrm{\Delta }}\mathrm{\Delta }g_\mathrm{r}\varphi _\mathrm{m}`$ agrees well with the BCS theory (‘BCS’ in the figure). Since $`\stackrel{~}{\mathrm{\Delta }}`$ determines the energy gap of the fermion quasi-particle spectrum as $`E_\text{p}=\sqrt{\xi _\text{p}^2+\stackrel{~}{\mathrm{\Delta }}^2}`$, the superfluid character at $`\nu /\epsilon _\mathrm{F}=2`$ is found to be the BCS-type. However, Fig. 6(a) also shows a sizable difference between $`\stackrel{~}{\mathrm{\Delta }}`$ and the Cooper-pair amplitude $`\mathrm{\Delta }=_\text{p}c_\text{p}c_\text{p}`$. This is because the BEC order parameter $`\varphi _\mathrm{m}=b_{\text{q}=0}`$ is always induced below $`T_\mathrm{c}`$ due to the Feshbach coupling effect in Eq. (8), which contributes to $`\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }g_\mathrm{r}\varphi _\mathrm{m}`$. As illustration, when we substitute $`\nu /\epsilon _\mathrm{F}=2`$, $`g_\mathrm{r}/\epsilon _\mathrm{F}=0.6`$, $`U/\epsilon _\mathrm{F}=0.3`$ and $`\mu /\epsilon _\mathrm{F}0.84`$ into Eq. (8), the contribution of the condensed $`b`$-boson to $`\stackrel{~}{\mathrm{\Delta }}`$ is found to be $`g_\mathrm{r}\varphi _\mathrm{m}=0.517\mathrm{\Delta }`$ at $`T=0`$. This means about one third of the order parameter $`\stackrel{~}{\mathrm{\Delta }}`$ is associated with $`\varphi _\mathrm{m}`$ even in the BCS regime at $`\nu /\epsilon _\mathrm{F}=2`$.

As the threshold energy $`2\nu `$ is lowered, the Cooper-pair component $`\mathrm{\Delta }`$ becomes less dominant while the condensed $`b`$-boson component $`\varphi _\mathrm{m}`$ increases (Fig. 6(a) $``$ (d)). At the same time, the temperature dependence of $`\stackrel{~}{\mathrm{\Delta }}`$ deviates from the weak-coupling BCS theory. At $`\nu /\epsilon _\mathrm{F}=2`$ in Fig. 6(d), $`\stackrel{~}{\mathrm{\Delta }}`$ is dominated by $`\varphi _\mathrm{m}`$ (although the Cooper-pair component also exists due to the Feshbach coupling effect). In this case, $`\varphi _\mathrm{m}`$ is well described by the BEC order parameter of $`N/2`$ atoms in a free Bose gas, given by

$$\varphi _\mathrm{m}^{\mathrm{BEC}}=\sqrt{\frac{N}{2}}\sqrt{1\left(\frac{T}{T_\mathrm{c}}\right)^{3/2}}.$$
(195)
We denote $`\varphi _\mathrm{m}^{\mathrm{BEC}}`$ in Fig. 6 as ‘BEC.’

In the intermediate region of the BCS-BEC crossover, $`\stackrel{~}{\mathrm{\Delta }}`$ (and also $`\mathrm{\Delta }`$ and $`\varphi _\mathrm{m}`$) becomes double-valued near $`T_\mathrm{c}`$. For example, in Fig. 6(b), the two self-consistent values for $`\stackrel{~}{\mathrm{\Delta }}`$ at $`T_\mathrm{c}=0.213T_\mathrm{F}`$ ($`T_\mathrm{c}^\mathrm{L}`$, where L stands for lower transition temperature) are zero and $`0.14\epsilon _\mathrm{F}`$. In this case, when one decreases the temperature, $`\stackrel{~}{\mathrm{\Delta }}`$ jumps abruptly at $`T_\mathrm{c}^\mathrm{L}`$. This is because the phase transition is suppressed by superfluid particle-particle fluctuations, once the phase transition occurs, the opening up of the fermion quasi-particle excitation gap ($`2\stackrel{~}{\mathrm{\Delta }}`$) strongly suppresses these fluctuation effects, which accelerates the increase of $`\stackrel{~}{\mathrm{\Delta }}`$. On the other hand, when we raise the temperature from below, since the superfluid fluctuations are suppressed by the excitation gap, we can exceed $`T_\mathrm{c}^\mathrm{L}`$ staying in the superfluid phase up to a higher temperature $`(T_\mathrm{c}^\mathrm{H})`$. In Fig. 6(b), we see that $`T_\mathrm{c}^\mathrm{H}=0.215T_\mathrm{F}`$, slightly higher than $`T_\mathrm{c}^\mathrm{L}`$. At $`T_\mathrm{c}^\mathrm{H}`$, $`\stackrel{~}{\mathrm{\Delta }}`$ vanishes discontinuously.

Fig. 7 shows that the chemical potential $`\mu `$ is decreased as the system approaches the BEC regime ($`\nu <0`$). The temperature dependence of $`\mu `$ is found to be weak below $`T_\mathrm{c}`$ (except just below $`T_\mathrm{c}`$ in the crossover regime) compared with the normal phase. As discussed in Section III, the chemical potential is temperature-independent below $`T_\mathrm{c}`$ in both the BCS limit ($`\mu =\epsilon _\mathrm{F}`$) and the BEC limit ($`\mu =\nu `$). Fig. 7 shows that this feature also holds in the intermediate region of the BCS-BEC crossover except just below $`T_\mathrm{c}`$, at least for the model parameters we have chosen.

Now that we have calculated self-consistently the values of the composite order parameter $`\stackrel{~}{\mathrm{\Delta }}`$ as function of both $`T`$ and $`2\nu `$, we can use the results to discuss the spectrum of the BCS-Bogoliubov single-particle Green’s function. The excitation spectrum of the BCS-Bogoliubov quasi-particles is given by

$$\rho _\mathrm{F}(\omega )\frac{1}{\pi }\underset{\text{p}}{}Im[\widehat{G}_{11}(𝐩,i\omega _m\omega +i\delta )],$$
(196)
where $`\widehat{G}_{11}(𝐩,\omega +i\delta )`$ is the (11)-component of Eq. (27). When the chemical potential $`\mu `$ is positive, the excitation spectrum in Eq. (196) has a finite gap $`\stackrel{~}{\mathrm{\Delta }}`$ as in BCS theory,

$`\rho _\mathrm{F}(\omega ,\mu 0)={\displaystyle \frac{3N}{8\epsilon _\mathrm{F}^{3/2}}}`$ $`[`$ $`\left({\displaystyle \frac{\omega }{\sqrt{\omega ^2\stackrel{~}{\mathrm{\Delta }}^2}}}+1\right)\sqrt{\mu +\sqrt{\omega ^2\stackrel{~}{\mathrm{\Delta }}^2}}`$ (197)
$`+`$ $`({\displaystyle \frac{\omega }{\sqrt{\omega ^2\stackrel{~}{\mathrm{\Delta }}^2}}}1)\sqrt{\mu \sqrt{\omega ^2\stackrel{~}{\mathrm{\Delta }}^2}}\mathrm{\Theta }(\sqrt{\stackrel{~}{\mathrm{\Delta }}^2+\mu ^2}\omega )]\mathrm{\Theta }(\omega \stackrel{~}{\mathrm{\Delta }}),`$ (198)
where $`\mathrm{\Theta }(x)`$ is the step function. On the other hand, when $`\mu `$ is negative in the BEC regime, one finds,

$`\rho _\mathrm{F}(\omega ,\mu <0)={\displaystyle \frac{3N}{8\epsilon _\mathrm{F}^{3/2}}}\left({\displaystyle \frac{\omega }{\sqrt{\omega ^2\stackrel{~}{\mathrm{\Delta }}^2}}}+1\right)\sqrt{\mu +\sqrt{\omega ^2+\stackrel{~}{\mathrm{\Delta }}^2}}\mathrm{\Theta }(\omega \sqrt{\stackrel{~}{\mathrm{\Delta }}^2+\mu ^2}).`$ (199)
In this case, the Fermi quasiparticle excitation gap is found to be $`\sqrt{\stackrel{~}{\mathrm{\Delta }}^2+\mu ^2}`$, rather than $`\stackrel{~}{\mathrm{\Delta }}`$. This reflects the fact that the threshold energy of free fermion excitations is given by $`|\mu |`$ when $`\mu <0`$. In the BEC limit ($`\mu \nu \epsilon _\mathrm{F}`$), we have $`|\mu |\stackrel{~}{\mathrm{\Delta }}`$ and the excitation gap reduces to $`|\mu |`$.

Fig. 8(a) shows the excitation spectrum of BCS-Bogoliubov quasi-particles in the BCS-BEC crossover. In the BCS regime ($`\nu \epsilon _\mathrm{F}`$), we find a peak at the excitation edge $`\omega =\stackrel{~}{\mathrm{\Delta }}`$ in the spectrum. This is the well-known coherence peak discussed in the superconductivity literature, and the quasi-particle spectrum is found to be the BCS-type in this regime. This coherence is absent when $`\nu 0`$, where the excitation gap gradually changes from $`\stackrel{~}{\mathrm{\Delta }}`$ to $`|\mu |`$ as the threshold energy $`2\nu `$ is lowered. Since $`\mu \nu `$ in the BEC regime, the energy gap ($`\sqrt{\mu ^2+\stackrel{~}{\mathrm{\Delta }}^2}`$) becomes larger for lower values of $`\nu `$ in the BEC regime.

Fig. 8(b) shows the momentum distribution function of Fermi atoms at $`T=0`$, which is given by

$$v_𝐩^2c_{\text{p}\sigma }^{}c_{\text{p}\sigma }=\frac{1}{2}\left(1\frac{\xi _\text{p}}{E_\text{p}}\right).$$
(200)
Since the energy gap $`\stackrel{~}{\mathrm{\Delta }}`$ in $`E_\text{p}=\sqrt{\xi _\text{p}^2+\stackrel{~}{\mathrm{\Delta }}^2}`$ is larger for smaller values of the threshold energy $`2\nu `$, the steep decrease of $`v_𝐩^2`$ around $`\mu (T=0)`$ at $`\nu /\epsilon _\mathrm{F}=2`$ gradually disappears as $`2\nu `$ is lowered. In addition, the magnitude of $`v_𝐩^2`$ decreases as one approaches the BEC regime due to the decrease of the chemical potential $`\mu (T=0)`$, which reflects the fact that that most Fermi atoms form $`b`$-bosons in the BEC regime. The quantity $`v_𝐩`$ also enters into the Bogoliubov-transformation to BCS quasiparticles

$$\gamma _\text{p}^{}=u_\text{p}c_\text{p}^{}+v_\text{p}c_\text{p},$$
(201)
where $`\gamma _\text{p}^{}`$ is a creation operator of a BCS-Bogoliubov quasi-particle and $`u_\text{p}^21v_\text{p}^2`$.

### B Velocity of the Goldstone phonon mode

Fig. 9 shows the velocity of the Goldstone mode $`v_\varphi `$ at $`T=0`$ as obtained from Eq. (178). In the BCS regime ($`\nu >\epsilon _\mathrm{F}`$), the mode velocity $`v_\varphi `$ agrees with the well-known Anderson-Bogoliubov phonon velocity $`v_\varphi =\overline{v}_\mathrm{F}/\sqrt{3}`$ in Eq. (184). \[At $`\nu /\epsilon _\mathrm{F}=2`$, we obtain $`\overline{v}_\mathrm{F}=\sqrt{2m\mu (0)}=0.92v_\mathrm{F}`$ for $`\mu (0)0.84\epsilon _\mathrm{F}`$, which gives $`v_\varphi =0.53v_\mathrm{F}`$.\] As the threshold energy $`2\nu `$ is lowered, $`v_\varphi `$ decreases sharply and approaches the Bogoliubov phonon mode given by Eq. (190). Fig. 9 indicates that $`v_\varphi `$ is strongly dependent on the threshold energy $`2\nu `$ in a uniform Fermi gas.

Fig. 10 shows the velocity of the Goldstone mode at finite temperatures, obtained from Eq. (175). In Fig. 10(a), $`v_\varphi `$ is found to be well described by the Anderson-Bogoliubov mode in the weak-coupling BCS theory (‘BCS’ in the figure) in the whole temperature region, as expected for the value $`\nu =2\epsilon _\mathrm{F}`$. On the other hand, $`v_\varphi `$ becomes less than the BCS result for the Anderson-Bogoliubov mode as the threshold energy $`2\nu `$ is decreased. This is shown in Fig. 10(b). Since the order parameter $`\stackrel{~}{\mathrm{\Delta }}`$ vanishes discontinuously due to the fluctuation effect discussed in the previous subsection, $`v_\varphi `$ shows a finite jump at $`T_\mathrm{c}`$ in Fig. 10(b).

### C Dispersion relation of the Goldstone mode

Fig. 11 shows the dispersion of the Goldstone mode at $`T=0.5T_\mathrm{c}`$. In the BCS regime (Fig. 11(a)), the gapless dispersion is convex and is confined below the excitation gap at $`2\stackrel{~}{\mathrm{\Delta }}`$. This convex dispersion relation gradually changes to a concave one as one goes from Figs. 11(a) to 11(d). In the BEC limit ($`\nu <0`$), the dispersion is particle-like ($`q^2`$), characteristic of free $`b`$-bosons. Indeed the dispersion in Fig. 11(d) is well described by $`\omega =q^2/2M`$, except for the region of linear dispersion $`\omega =v_\varphi q`$ around $`𝐪=0`$.

In the BCS regime, Eq. (176) also has a high-energy solution at $`\omega /\epsilon _\mathrm{F}2.4`$, as shown in Fig. 12. This solution is also obtained at $`T_\mathrm{c}`$, as shown in Fig. 13(a). Since the energy of this solution is close to the threshold energy of the excitation spectrum of free $`b`$-bosons as $`\xi _{B\text{q}=0}=2\nu 2\mu =2.32\epsilon _\mathrm{F}`$ ($`\mu 0.84\epsilon _\mathrm{F}`$), this high energy solution is interpreted as an excitation of a $`b`$-boson.

The high-energy mode solution at $`T_\mathrm{c}`$ shown in Fig. 13(a) may be useful as a guide to see how the Goldstone mode changes in the BCS-BEC crossover. Fig. 13(a) shows that the energy of this mode gradually decreases as the threshold energy $`2\nu `$ is lowered. At the same time, other solutions appear around $`\omega =𝐪=0`$, as shown for $`\nu /\epsilon _\mathrm{F}=1.29`$ in panel (b). At $`\nu /\epsilon _\mathrm{F}=1.14`$ and below, two finite energy and one gapless excitations combine to produce a single phonon. Then the Goldstone mode dispersion below $`T_\mathrm{c}`$ starts to change from convex to concave. In Fig. 11(b), although the dispersion below $`2\stackrel{~}{\mathrm{\Delta }}`$ is still convex, characteristic of the Anderson-Bogoliubov mode in BCS superfluidity, the dispersion above $`2\stackrel{~}{\mathrm{\Delta }}`$ obtained from Eq. (176) approaches the concave dispersion obtained at $`T_\mathrm{c}`$. In Figs. 11(c) and (d), the region where linear dispersion occurs ($`\omega =v_\varphi q`$) is very narrow. The overall behavior is the same as the dispersion relation obtained at $`T_\mathrm{c}`$, which is particle-like, $`q^2/2M`$.

Finally we briefly comment on the approximate Eqs. (175) and (176). Since the excitation gap $`2\stackrel{~}{\mathrm{\Delta }}`$ strongly suppresses the fermion quasi-particle excitations far below $`T_\mathrm{c}`$, the intraband terms of $`\mathrm{\Pi }_{ij}^0`$ and also the Landau damping below $`2\stackrel{~}{\mathrm{\Delta }}`$ are less important at $`T=0.5T_\mathrm{c}`$. As a result, Eqs. (175) and (176) both give good approximations to the solution of Eq. (173) at $`T=0.5T_\mathrm{c}`$, and give almost the same results below $`2\stackrel{~}{\mathrm{\Delta }}`$, as shown in Figs. 11(a) and (b). On the other hand, the interband terms given by the second lines in Eqs. (43)-(53) give rise to Landau damping when $`\omega 2\stackrel{~}{\mathrm{\Delta }}`$. Thus, above $`2\stackrel{~}{\mathrm{\Delta }}`$, Eqs. (175) and (176) may have different solutions, because they do not include the imaginary part in the same way. Indeed, the dashed line above $`2\stackrel{~}{\mathrm{\Delta }}`$ shown in Fig. 11(b), which is a solution of Eq. (176), is not obtained from Eq. (175). In such a case, a more careful analysis is necessary to obtain the correct dispersion relation. However, since the fermion quasi-particles are absent in the BEC limit ($`\nu <0`$), the Landau damping becomes weak in the BEC regime even near $`T_\mathrm{c}`$. Thus Eqs. (175) and (176) again give almost identical results in the whole energy region shown in Figs. 11(c) and (d).

### D Spectral weight and damping of the Goldstone mode

In Sections VI.B and VI.C, we considered the Goldstone mode neglecting the Landau damping. In this section, we evaluate the damping from the width of the collective mode in the structure function $`S_{jj}(𝐪,\omega _+)`$ in Eq. (177) as well as in the $`b`$-boson spectral density

$$\rho _\mathrm{B}(𝐪,\omega )\frac{1}{\pi }Im[\widehat{D}_{11}(𝐪,\omega _+)].$$
(202)
In Figs. 14-17, we show the $`b`$-boson excitation spectrum, as well as the phase ($`S_{22}`$) and amplitude ($`S_{11}`$) structure functions for $`T`$, below $`T_\mathrm{c}`$. In Fig. 14, we find that the Anderson-Bogoliubov Goldstone mode does not appear as a visible peak in the spectrum at $`T/T_\mathrm{c}=0.99`$. (The inset in Fig. 14 shows that the mode energy is $`\omega /\epsilon _\mathrm{F}=0.01`$ for $`q=0.02p_\mathrm{F}`$.) This means that the Anderson-Bogoliubov mode is over-damped near $`T_\mathrm{c}`$ because of strong Landau damping from the thermally-excited fermion quasi-particles. We note that $`S_{22}(𝐪,\omega )`$ and $`S_{11}(𝐪,\omega )`$ both exhibit strong central peaks at $`\omega =0`$ for $`T/T_\mathrm{c}=0.99`$, which indicates the presence of a large number of thermally-excited fermion quasi-particles expected at this temperature.

At low temperatures, when the Landau damping becomes weaker, the Anderson-Bogoliubov mode appears as a visible peak in $`\rho _\mathrm{B}(𝐪,\omega )`$ and $`S_{22}(𝐪,\omega )`$ as shown in Fig. 14. The peak width becomes narrower at lower temperatures, reflecting the weaker Landau damping by fermions. (At $`T=0`$, it becomes a sharp $`\delta `$-function peak.) The peak position ($`\omega /\epsilon _\mathrm{F}0.021`$) at $`T/T_\mathrm{c}=0.5`$ agrees well with the dispersion relation shown in Fig. 11(a).

Since the Anderson-Bogoliubov mode is a collective phase oscillation of the Cooper-pair order parameter $`\mathrm{\Delta }`$, the appearance of this mode in $`\rho _\mathrm{B}(𝐪,\omega )`$ indicates the presence of the coupling between the $`b`$-bosons and the phase fluctuations of $`\mathrm{\Delta }`$. On the other hand, no peak structure is observed except for the central peak at $`\omega =0`$ in the amplitude structure function $`S_{11}`$ shown in Fig. 14. Only a slight structure appears at $`\omega /\epsilon _\mathrm{F}0.02`$ at $`T/T_\mathrm{c}=0.7`$, as shown in part (a) of Fig. 18. (The peak at $`\omega /\epsilon _\mathrm{F}0.07`$ in Fig. 18 is the amplitude mode at the edge of interband excitations $`\omega =2\stackrel{~}{\mathrm{\Delta }}`$.) This is because the amplitude-phase coupling $`\mathrm{\Pi }_{12}^0`$ is very weak in the BCS regime, so that the collective phase oscillation does not strongly couple into the amplitude fluctuations.

As shown in Fig. 12, Eq. (176) also has a high energy solution in the BCS regime at $`\nu /\epsilon _\mathrm{F}=2`$. This solution is the strong resonance in $`\rho _\mathrm{B}(𝐪,\omega )`$, as shown in part (b) of Fig. 18.

In Fig. 15, we find a broad peak in $`S_{22}`$ at $`T/T_\mathrm{c}=0.99`$. As shown in the inset in Fig. 15, the peak position is different from the one expected from Eqs. (175) and (176). However, since this resonance shows gapless behavior, it clearly must be the Goldstone mode. (The difference shown in the inset is due to the Landau damping effect.) Indeed, the peak energy at $`T/T_\mathrm{c}=0.5`$ agrees well with the dispersion in Fig. 11(b). As expected, Landau damping from fermions becomes weak as one approaches the BEC regime. When $`\nu 0`$, we can observe the sharp peak structure even near $`T_\mathrm{c}`$ in Figs. 16 and 17, and the peak position always agrees well with the dispersion relation given in Fig. 1. The dispersion relation of the Goldstone mode approaches the temperature-independent particle-like one in the BEC regime, and hence the temperature dependence of the peak energy in Fig. 17 is weak.

We also find in Figs. 15-17 that the Goldstone mode appears in the amplitude structure factor $`S_{11}(𝐪,\omega )`$, which indicates that the amplitude-phase coupling becomes stronger as $`2\nu `$ is decreased. Although the amplitude fluctuations described by $`\mathrm{\Pi }_{11}`$ are not important in the BCS regime, we cannot neglect these fluctuations in the BCS-BEC crossover regime.

## VII Density-density correlation function

An important problem is how to experimentally observe the Goldstone mode discussed in Section V and VI. In this section, we show that the density-density correlation function $`\mathrm{\Pi }_{33}`$ defined in Eq. (41) exhibits this mode as a pole. This correlation function can be experimentally probed by many techniques, including two-photon Bragg scattering.

Since the Goldstone mode is associated with the collective phase fluctuations of the Cooper-pairs, a coupling between density fluctuations and phase fluctuations is necessary in order for the Goldstone mode to appear in the spectrum exhibited by $`\mathrm{\Pi }_{33}(𝐪,\omega _+)`$. In this regard, we note that a phase-density coupling exists in fermion superfluidity because of the presence of the Josephson effect. This coupling is described by $`\mathrm{\Pi }_{23}^0`$ (and $`\mathrm{\Pi }_{32}^0`$) defined in Eq. (51). In contrast to the amplitude-phase correlation function $`\mathrm{\Pi }_{12}^0`$, which is very weak in the BCS regime, $`\mathrm{\Pi }_{23}^0`$ is finite even if one works with the approximation in Eq. (183). Besides $`\mathrm{\Pi }_{23}^0`$, density fluctuations also couple with superfluid fluctuations through amplitude-density coupling $`\mathrm{\Pi }_{13}^0`$ (and $`\mathrm{\Pi }_{31}^0`$) defined in Eq. (53). This coupling is not important in the BCS regime, because $`\mathrm{\Pi }_{13}^0`$ vanishes when one uses the approximation in Eq. (183) and neglects term of order $`O((q/p_\mathrm{F})^2)`$). However, as in the case of the amplitude-phase correlation function $`\mathrm{\Pi }_{12}^0`$, this coupling effect becomes stronger as one goes into the BCS-BEC crossover regime. Then the density-density correlation function $`\mathrm{\Pi }_{33}`$ is found to couple with superfluid fluctuations (amplitude and phase) via $`\mathrm{\Pi }_{23}^0`$ and $`\mathrm{\Pi }_{13}^0`$, as shown diagrammatically in Fig. 19.

The density-density correlation function $`\mathrm{\Pi }_{33}`$ given by the HF-RPA is obtained by extending the method discussed in Section V. When we introduce a $`3\times 3`$-matrix correlation function $`\widehat{\mathrm{\Pi }}=\{\mathrm{\Pi }_{ij}\}(i,j=1,2,3)`$, which involves the density-density correlation function as the $`(33)`$-component, we obtain an equation similar to Eq. (160), namely

$$\widehat{\mathrm{\Pi }}=\widehat{\mathrm{\Pi }}^0\left[1+\frac{1}{2}\overline{U}\widehat{\mathrm{\Pi }}^0\right]^1.$$
(203)
Here the interaction $`3\times 3`$-matrix $`\overline{U}`$ is defined by

$`\overline{U}\left(\begin{array}{ccc}\overline{U}_{11}& \overline{U}_{12}& 0\\ \overline{U}_{21}& \overline{U}_{22}& 0\\ 0& 0& 0\end{array}\right),`$ (207)
where $`\overline{U}_{ij}\{Ug_\mathrm{r}^2\widehat{W}\widehat{D}^0\widehat{W}^1\}_{ij}`$. If we included a density-density interaction term to our coupled fermion-boson model in Eq. (2), it would be included as the $`(33)`$-component in Eq. (207). The density-density correlation function is obtained from the $`(33)`$-component of Eq. (203),

$$\mathrm{\Pi }_{33}=\mathrm{\Pi }_{33}^0+\frac{\overline{\mathrm{\Pi }}_{33}}{1+\frac{V_1}{2}[\mathrm{\Pi }_{11}^0+\mathrm{\Pi }_{22}^0]\frac{V_1}{2}[\mathrm{\Pi }_{12}^0\mathrm{\Pi }_{21}^0]+[\frac{V_1^2}{4}+\frac{V_2^2}{4}][\mathrm{\Pi }_{11}^0\mathrm{\Pi }_{22}^0\mathrm{\Pi }_{12}^0\mathrm{\Pi }_{21}^0]},$$
(208)
where the numerator is

$`\overline{\mathrm{\Pi }}_{33}`$ $``$ $`[{\displaystyle \frac{V_1^2}{4}}+{\displaystyle \frac{V_2^2}{4}}][\mathrm{\Pi }_{31}^0\mathrm{\Pi }_{12}^0\mathrm{\Pi }_{23}^0\mathrm{\Pi }_{31}^0\mathrm{\Pi }_{22}^0\mathrm{\Pi }_{13}^0+\mathrm{\Pi }_{32}^0\mathrm{\Pi }_{21}^0\mathrm{\Pi }_{13}^0\mathrm{\Pi }_{32}^0\mathrm{\Pi }_{11}^0\mathrm{\Pi }_{23}^0]`$ (209)
$``$ $`{\displaystyle \frac{V_1}{2}}[\mathrm{\Pi }_{31}^0\mathrm{\Pi }_{13}^0+\mathrm{\Pi }_{32}^0\mathrm{\Pi }_{23}^0]{\displaystyle \frac{V_2}{2}}[\mathrm{\Pi }_{31}^0\mathrm{\Pi }_{23}^0\mathrm{\Pi }_{32}^0\mathrm{\Pi }_{13}^0].`$ (210)
Each term in Eq. (210) involves correlation functions between density fluctuations and superfluid fluctuations, such as $`\mathrm{\Pi }_{23}^0`$ and $`\mathrm{\Pi }_{13}^0`$.

Fig. 20 shows the dynamic structure function related to the spectral density of density-density correlation function,

$$S_{33}(𝐪,\omega )\frac{1}{\pi }[n_\mathrm{B}(\omega )+1]Im[\mathrm{\Pi }_{33}(𝐪,i\nu _n\omega _+)].$$
(211)
At low temperatures in Figs. 20(a) and (b) (BCS regime), and at all the temperatures in Figs. 20(c) and (d) (BEC regime), we can clearly see the Goldstone mode as a resonant peak in $`S_{33}(𝐪,\omega )`$. The peak position is the same as that in the spectra shown in Figs. 14-17. Since the density-density correlation function $`\mathrm{\Pi }_{33}`$ can be measured more easily than the phase correlation function $`\mathrm{\Pi }_{22}`$, and may be the most useful way of observing the Goldstone mode.

## VIII Summary and Discussion

In this paper, we have investigated the BCS-BEC crossover in the superfluid phase of a uniform gas of Fermi atoms with a Feshbach resonance. Going past the simple weak-coupling mean-field theory, we included the strong-coupling effect originating from the pairing interaction associated with the Feshbach resonance. We have extended our previous work at and above $`T_\mathrm{c}`$ to the superfluid region below $`T_\mathrm{c}`$. We showed that the superfluid order parameter continuously changes from the Cooper-pair amplitude $`\mathrm{\Delta }=U_\text{p}c_\text{p}c_\text{p}`$ in the BCS regime to the square root of the number of condensed $`b`$-molecules $`\varphi _\mathrm{m}=b_{\text{q}=0}`$ associated with the Feshbach resonance in the BEC regime, as one lowers the threshold energy $`2\nu `$ of the Feshbach resonance. In the intermediate regime in the BCS-BEC crossover, superfluidity is described by the composite order parameter $`\stackrel{~}{\mathrm{\Delta }}\mathrm{\Delta }g_\mathrm{r}\varphi _\mathrm{m}`$.

The Goldstone mode is one of the most fundamental phenomena in an ordered system with spontaneous breakdown of a continuous symmetry. In this paper, we investigated how the Anderson-Bogoliubov mode, which is the Goldstone mode in BCS superfluidity, changes to the Bogoliubov phonon mode in the BCS-BEC crossover. The velocity of the Goldstone mode $`v_\varphi `$ strongly depends on the threshold energy $`2\nu `$, and decreases as one approaches the BEC regime. The Anderson-Bogoliubov phonon mode may be a useful way of monitoring the BCS-BEC crossover. Since it is difficult to strongly modify the strength of the pairing interaction in metallic superconductors, the tunable pairing interaction associated with the Feshbach resonance in Fermi atomic gases gives one a unique tool to clarify the physics in the BCS-BEC crossover region.

We also investigated the damping of the Goldstone mode. The fermion Landau damping of the Goldstone mode due to coupling to fermions becomes weak far below $`T_\mathrm{c}`$, reflecting the fact that the thermal excitation of fermion quasi-particles is negligible. This damping effect is always weak in the BEC regime because the system is then dominated by $`b`$-molecules, with the suppression of the Fermi quasiparticle spectrum. Thus, except in a small region near $`T_\mathrm{c}`$ in the BCS regime, the Goldstone mode appears as a strong resonance both in the spectrum of the phase correlation function and in the excitation spectrum of $`b`$-bosons. As a way to experimentally observe this Goldstone mode, we noted that since the amplitude and phase fluctuations couple with the density fluctuations, one can observe this collective mode in the spectrum of the density-density correlation function $`\mathrm{\Pi }_{33}(𝐪,\omega _+)`$. In cigar-shaped trapped Fermi gases, the direct observation of a density fluctuation pulse might be possible, in analogy with the observation of a Bogoliubov phonon in a superfluid Bose gas with a very weak axial trapping potential. We also note that two photon Bragg scattering experiments provide a convenient way of studying density fluctuations in a trapped atomic gas. In such experiments, $`Im\mathrm{\Pi }_{33}(𝐪,\omega )`$ is measured directly, rather than $`S_{33}(𝐪,\omega )`$.

It is important to remember that our treatment of the BCS-BEC crossover leaves out several important contributions. In our calculation of various response functions, we always ignored the fermion self-energy arising from the coupling to the non-Bose-condensed bosons. In the crossover region, these two-particle states are strongly damped. The numerical results given in Section VI and VII show how Landau damping due to fermions decreases as we approach the BEC regime (small or negative values of $`\nu /\epsilon _\mathrm{F}`$). This is simply because more and more the fermions are forming bound states. However, we expect new forms of damping to arise in the BEC region, namely the Landau and Baliaev damping associated with the interaction between Bogoliubov excitations.

In this paper, we have not considered the effect of a trap. The atomic density profile in the BCS-BEC crossover in a trap was recently investigated within the LDA by the authors at and above $`T_\mathrm{c}`$, using the coupled fermion-boson model, and at $`T=0`$ by Perali et al., using the strong-coupling BCS model. A trap potential also leads to various discrete low energy collective modes associated with the confined geometry, and it will be interesting to study how these collective oscillations behave in the BCS-BEC crossover. This will be the subject of a future paper, but we conclude with a few brief remarks on this. In the BCS limit of our model, one has the extensive theoretical work on the Cooper-pair condensate modes in trapped Fermi gases. To be specific, we consider the quadrupole mode $`\omega _Q=\sqrt{2}\omega _0`$ (where $`\omega _0`$ is the trap frequency of a spherical trap). This is an analogue of the Goldstone phonon mode in a uniform system. If $`\omega _Q`$ is greater than the effective threshold $`2\mathrm{\Delta }_{\mathrm{eff}}`$ for breaking up Cooper-pairs, this mode is damped and has small spectral weight. However, as one approaches the BCS-BEC crossover, the spectral weight of the fermion quasi-particles decreases, shifting to the quadrupole mode $`\omega _Q`$. In the BEC limit, the mode $`\omega _Q=\sqrt{2}\omega _0`$ is the well known quadrupole oscillation of a trapped Bose gas. In this paper, we have seen that the Goldstone phonon mode in a uniform gas also persists in the BCS-BEC crossover, but the phonon velocity changes. In contrast, the quadrupole mode frequency does not depend explicitly on the interaction strength, having the same frequency $`\omega _Q=\sqrt{2}\omega _0`$ in both the BCS and BEC limits. What does change is its spectral weight and damping, which could be used as an experimental signature.

###### Acknowledgements.

Y. O. would like to thank Professors S. Takada and M. Ueda for useful discussions, and the Japanese Government for financial support while a visiting professor at the University of Toronto. A. G. acknowledges financial support from the NSERC of Canada.