Dynamic structure factor of a strongly correlated Fermi superfluid within a density functional theory approach

The realization of ultracold Fermi gases of ⁶Li and ⁴⁰K atoms near Feshbach resonances provides a new paradigm for studying strongly correlated many-body systems [1]. At low temperature, these systems display an intriguing crossover from Bardeen–Cooper–Schrieffer (BCS) superfluids to Bose–Einstein condensates (BEC) [2, 3]. At a special point in between the two regimes, where the s-wave scattering length diverges, the gas exhibits universal properties, which might also exist in other strongly interacting Fermi superfluids [4, 5], such as high-temperature superconductors or nuclear matter in neutron stars. This is called unitary Fermi gas and corresponds a novel type of superfluid with neither dominant bosonic nor fermionic character. This new superfluid has already been intensively investigated [3], leading to several milestone observations.

Theoretical challenges in describing the BCS-BEC crossover arise from its strongly correlated nature: there is no small interaction parameter to control the accuracy of theories [6]. To date, significant progress has been made in developing better quantum Monte Carlo (QMC) simulations [7–15] and strong-coupling theories [16–24], leading to the quantitative establishment of a number of properties in conjunction with the rapid experimental advances. These include the equation of state [6, 25–29], frequency of collective oscillations [30, 31], pairing gap [10, 11, 32, 33], and superfluid transition temperature [9, 29]. However, some fundamental dynamical properties, such as the single-particle spectral function measured by radio-frequency (rf) spectroscopy [34–37] and the dynamic structure factor probed by Bragg spectroscopy [38–42], are not well understood yet.

As an important fingerprint of quantum gases in some certain states, dynamic structure factor contains rich information of properties of a many-body system [43]. By tuning the transferred momentum or energy from a low value to a high one, we can observe the low-lying collective phonon excitations, Cooper-pair (i.e., molecular) excitations and single-particle atomic excitations, respectively. In particular, at finite temperature the dynamic structure factor can help to judge whether the system is in the superfluid or normal state from the emergence of the phonon excitations. Also, it is reasonable to anticipate that the dynamic structure factor may play a role to solve the debate on the existence of pseudogap pairing or pre-pairing states [34, 35]. Experimentally, the dynamic structure factor can be measured via two-photon Bragg scattering technique [39], at both low and finite temperatures [41]. Theoretically, since no exact solution exists for strongly interacting Fermi gases and the numerically exact QMC approach is less efficient for simulating dynamical quantities, one has to resort to some approximated approaches, which are useful in certain limiting cases [40] (see table 1). For example, at high temperature, as the fugacity is a small parameter, a quantum cluster expansion has been proven to be an efficient method [44, 45], and has been used to calculate the dynamic structure factor [46, 47]. In the limits of both large momentum and high frequency, asymptotically exact Tan relations have been derived to describe the high-frequency tails [48, 49]. On the other hand, in the limit of long wavelength or small momentum, the phenomenological two-fluid hydrodynamic theory may provide a useful description [50].

Table 1.  A list of the existing theories for the dynamic structure of strongly interacting fermions, including the virial expansion [46, 47], Tan relation [48, 49], two-fluid hydrodynamics [50], diagrammatic strong-coupling approach [52, 53], and BdG-RPA [38, 55–57]. The applicable conditions for the transferred momentum q and temperature T, under which each theory is quantitatively useful, are indicated. Here ${k}_{F}={(3{\pi }^{2}n)}^{1/3}$, ${\varepsilon }_{F}={k}_{F}^{2}/{(2m)=(3{\pi }^{2}n)}^{2/3}/(2m),$ and T_F are the Fermi momentum, energy, and temperature, respectively. T_c is the superfluid transition temperature.

Theories	q (applicable)	T (applicable)
Virial expansion	arbitrary	$T\gt {T}_{F}$
Tan relation	$q\gg {k}_{F}$, $\omega \gg {\varepsilon }_{F}$	$T\gg {T}_{F}$
Two-fluid hydrodynamics	$q\ll {k}_{F}$,	$T\lt {T}_{c}$
Diagrammatic approach	arbitrary	arbitrary
BdG-RPA	$q\gg {k}_{F}$	$T\ll {T}_{c}$
SLDA-RPA (this work)	$q\lesssim {k}_{F}$	$T\ll {T}_{c}$

A general theoretical framework of the dynamic structure factor, valid at arbitrary temperature and momentum, can be developed by using the diagrammatic technique [51, 52] or functional path integral approach [53], in parallel with the existing strong-coupling theories of interacting Fermi gases [6]. The expressions for the density and spin responses of strongly interacting Fermi gases have been obtained [53]. However, their numerical calculations turn out to be extremely difficult, except in the limit of zero transferred momentum [52]. A microscopic approach utilizing correlated basis function theory, with static structure factor from QMC as input, has also been used to calculate the dynamic structure factor [54]. A more commonly used approach is the random-phase approximation (RPA) on top of the mean-field Bogoliubov–de Gennes (BdG) theory [38, 55–57]. By comparing the BdG-RPA predictions with the experimental data for the dynamic structure factor of strongly interacting fermions and with the QMC results for the static structure factor [56], it has been surprisingly shown by two of the present authors that the BdG-RPA theory works quantitatively well at sufficiently large transferred momentum (i.e., $q\sim 5{k}_{F}$, where k_F is the Fermi momentum). At small transferred momentum, i.e., $q\lesssim {k}_{F}$, apparently, the BdG-RPA only provides a qualitative description of the dynamic structure factor at the BCS-BEC crossover, since both the sound velocity (associated with the phonon excitations) and pairing gap (associated with the single-particle fermionic excitations) are strongly over-estimated within the BdG framework [3].

In this work, we aim to develop a quantitative theory for the dynamic structure factor of strongly interacting fermions at low transferred momentum and at low temperature, which is amenable for numerical calculations. For this purpose, we adopt a superfluid local density approximation (SLDA) approach [58–60], within the framework of density functional theory (DFT) [62–64], as recently suggested by Bulgac and his co-workers. The SLDA theory assumes an energy density functional (i.e., a function of the density function) to describe a unitary Fermi superfluid and uses the QMC results for the chemical potential and order parameter as two important inputs. It can be well regarded as a better quasi-particle description than the mean-field BdG theory. It has been shown that at low-energy the SLDA theory provides useful results for the equation of state [60] and real-time dynamics [65, 66] of a strongly interacting Fermi superfluid.

Here we apply the random phase approximation on top of the SLDA theory. The use of SLDA in place of the standard BdG equations improves the predictions for the dynamic structure factor in the BCS-BEC crossover near unitarity. The static structure factor at small momentum transfer is in excellent agreement with the results of the latest QMC [14, 15, 55]. A more stringent test can be obtained in the near future by comparing our predictions with the experimental data [42], without any adjustable parameters.

Our paper is organized as follows. In the next section (section 2), we introduce the SLDA theory. In section 3, we review the main idea of RPA. The expression for the dynamic structure factor is derived in section 4. In sections 5 and 6, we present our main results of dynamic structure factor in the unitary limit and the crossover regime, respectively. Finally, section 7 is devoted to conclusions and outlooks. For convenience, we set ${\hslash }={k}_{B}=1$ in the following discussions.

The DFT developed by Hohenberg and Kohn [62], together with the local density approximation (LDA) by Kohn and Sham [63], is a powerful tool to understand the properties of many-electron systems. The DFT was initially used for electrons in the normal, non-superconducting state. It is based on the assumptions that there is a unique mapping between the external potential and the total wave function of the system (or the normal density), and that the exact energy of the system can be written as a density functional. A limitation of the DFT is that the exact form of the density functional is often not known. Therefore, approximated phenomenological functionals are introduced, which should be optimized for a specific system. Typically, those functionals rely on the Kohn-Sham orbitals [64] and thus can not effectively deal with superfluidity. The generalization of the DFT to superfluid cold-atom systems—referred to as SLDA as we mentioned earlier—was recently introduced by Bulgac and Yu [58–60]. This SLDA originates from a similar DFT previously used in the context of nuclear physics [59, 66].

A nice feature of ultracold fermions is that, in the unitary limit the form of the energy density functional is restricted by dimensional arguments. Another advantage is the availability of ab initio QMC results and accurate experimental data for both homogeneous and inhomogeneous systems, which can be used to fix the parameters of the density functional, as we shall see below.

For a superfluid atomic Fermi gas, two atoms with mass m in different spin state can form a Cooper pair. As a result, the system possesses an anomalous Cooper-pair density $\nu ({\bf{r}},t)$, in addition to the number density $n({\bf{r}},t)$. The energy density functional ${ \mathcal E }[\tau ({\bf{r}},t),n({\bf{r}},t),\nu ({\bf{r}},t)]$ of the system must include the kinetic density $\tau ({\bf{r}},t)$, number density $n({\bf{r}},t)$, and also the anomalous density $\nu ({\bf{r}},t)$ [59, 60]:

$$\begin{eqnarray}&&{ \mathcal E }[\tau ,n,\nu ]=\alpha \displaystyle \frac{1}{2m}\tau +\beta \displaystyle \frac{3{(3{\pi }^{2})}^{2/3}}{10m}{n}^{5/3}+\gamma \displaystyle \frac{1}{{{mn}}^{1/3}}| \nu {| }^{2},\end{eqnarray} \tag{ 1 }$$

where the kinetic density τ, number density n and anomalous density ν are given by,

$$\begin{eqnarray}&&\tau =2\displaystyle \sum _{{\bf{k}}}| {\rm{\nabla }}{v}_{{\bf{k}}}{| }^{2},n=2\displaystyle \sum _{{\bf{k}}}| {v}_{{\bf{k}}}{| }^{2},\nu =\displaystyle \sum _{{\bf{k}}}{u}_{{\bf{k}}}{v}_{{\bf{k}}}^{* },\end{eqnarray} \tag{ 2 }$$

and ${u}_{{\bf{k}}}({\bf{r}},t)$ and ${v}_{{\bf{k}}}({\bf{r}},t)$ are the Bogoliubov quasiparticle wavefunctions with ${\bf{k}}$ labeling the quasiparticle states. Three dimensionless constants, the effective mass parameter α, Hartree parameter β and pairing parameter γ, are introduced. These parameters are determined by requiring that the SLDA reproduces exactly the zero temperature chemical potential, pairing gap and energy per particle that are obtained by either QMC simulations or accurate experimental measurements for a uniform system [60, 66].

In the unitary limit at zero temperature, the simple form of the energy density functional equation (1) is inspired by the dimensional analysis: the first and third terms are the unique combination required by the renormalizablity of the theory [58]; while the second term is the only possible form allowed by the scale invariance at unitarity [4]. The above energy density functional has been successfully used by Bulgac and his co-workers to understand the thermodynamics [60] and dynamics [65, 66] of a unitary Fermi gas at zero temperature. It is reasonable to assume that the energy density functional equation (1) can be applied also away from unitarity, but close to it, and at non-zero temperature, but significantly below T_c.

As both the kinetic and anomalous densities diverge due to the use of a pairwise contact interaction, a regularization procedure is needed for the pairing gap and for the energy density [58]. After regularization, the energy density functional with regularized kinetic density ${\tau }_{c}({\bf{r}},t)$ and anomalous density ${\nu }_{c}({\bf{r}},t)$ takes the following form [60],

$$\begin{eqnarray}&&{ \mathcal E }=\alpha \displaystyle \frac{1}{2m}{\tau }_{c}+\beta \displaystyle \frac{3{(3{\pi }^{2})}^{2/3}}{10m}{n}^{5/3}+{g}_{\mathrm{eff}}| {\nu }_{c}{| }^{2},\end{eqnarray} \tag{ 3 }$$

where the effective coupling constant ${g}_{\mathrm{eff}}$ is given by

$$\begin{eqnarray}&&\displaystyle \frac{1}{{g}_{\mathrm{eff}}}=\displaystyle \frac{{{mn}}^{1/3}}{\gamma }-\displaystyle \sum _{| {\bf{k}}| \lt {\boldsymbol{\Lambda }}}\displaystyle \frac{m}{\alpha {{\bf{k}}}^{2}}.\end{eqnarray} \tag{ 4 }$$

We note that, ${g}_{\mathrm{eff}}$ scales to zero once the cut-off momentum Λ runs to infinity. The order parameter ${\rm{\Delta }}({\bf{r}},t)$ is related to the anomalous density ν by

$$\begin{eqnarray}&&{\rm{\Delta }}({\bf{r}},t)=-{g}_{\mathrm{eff}}{\nu }_{c}({\bf{r}},t).\end{eqnarray} \tag{ 5 }$$

The stationary SLDA equations for the quasiparticle wave functions are obtained by the standard functional minimization with respect to the variations ${u}_{{\bf{k}}}$ and ${v}_{{\bf{k}}}$. One obtains

$$\begin{eqnarray}\left[\begin{array}{cc}{{ \mathcal H }}_{s}-\mu & {\rm{\Delta }}\\ {{\rm{\Delta }}}^{* } & -{{ \mathcal H }}_{s}+\mu \end{array}\right]\left[\begin{array}{c}{u}_{{\bf{k}}}\\ {v}_{{\bf{k}}}\end{array}\right]={E}_{{\bf{k}}}\left[\begin{array}{c}{u}_{{\bf{k}}}\\ {v}_{{\bf{k}}}\end{array}\right],\end{eqnarray} \tag{ 6 }$$

with a single quasiparticle Hamiltonian

$$\begin{eqnarray}&&{{ \mathcal H }}_{s}=-\alpha \displaystyle \frac{{{\rm{\nabla }}}^{2}}{2m}+\beta \displaystyle \frac{(3{\pi }^{2}n){}^{2/3}}{2m}-\displaystyle \frac{| {\rm{\Delta }}{| }^{2}}{3\gamma {{mn}}^{2/3}},\end{eqnarray} \tag{ 7 }$$

and the chemical potential μ.

By requiring that a homogeneous Fermi gas of the number density $n=N/V={k}_{F}^{3}/(3{\pi }^{2})$ has an energy per particle $E/N=(3/5){\xi }_{E}{\varepsilon }_{F}$, a chemical potential $\mu ={\xi }_{\mu }{\varepsilon }_{F}$, and a pairing order parameter ${\rm{\Delta }}=\eta {\varepsilon }_{F}$ at zero temperature, one can determine the value of dimensionless parameters α, β and γ in equation (3) through the following equations, which are independent on the cut-off momentum (i.e., ${\rm{\Lambda }}\to \infty$):

$$\begin{eqnarray}&&n=\displaystyle \sum _{{\bf{k}}}\left(1-\displaystyle \frac{{\xi }_{{\bf{k}}}}{{E}_{{\bf{k}}}}\right),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{mn}}^{1/3}}{\gamma }=\displaystyle \sum _{{\bf{k}}}\left(\displaystyle \frac{m}{\alpha {{\bf{k}}}^{2}}-\displaystyle \frac{1}{2{E}_{{\bf{k}}}}\right),\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{3}{5}{E}_{F}n({\xi }_{E}-\beta )=\displaystyle \sum _{{\bf{k}}}\left[\alpha \displaystyle \frac{{{\bf{k}}}^{2}}{2m}\left(1-\displaystyle \frac{{\xi }_{{\bf{k}}}}{{E}_{{\bf{k}}}}\right)-\displaystyle \frac{{\rm{\Delta }}}{2{E}_{{\bf{k}}}}\right],\end{eqnarray} \tag{ 10 }$$

where ${\xi }_{{\bf{k}}}=\alpha {{\bf{k}}}^{2}/{(2m)+[\beta -(3{\pi }^{2})}^{2/3}{\eta }^{2}/(6\gamma )-{\xi }_{\mu }]{\varepsilon }_{F}$ and ${E}_{{\bf{k}}}=\sqrt{{\xi }_{{\bf{k}}}^{2}+{{\rm{\Delta }}}^{2}}$. In these three constraint equations, ${\xi }_{\mu }$, η and ${\xi }_{E}$ are the three inputs, whose value can be reliably determined by using QMC simulations [11–13] or from the experimental measurements [28, 29, 32]. The parameter α can be determined using the single particle dispersion, near the unitary limit, typically the parameter α is very close to 1 [60, 61], indicating that the effective mass only differs slightly from the bare atomic mass m. For simplicity, throughout the work we take $\alpha =1$ and use the density equation equation (8) and the gap equation equation (9) to determine the parameters β and γ. As we shall see, this simple choice also ensures that the f-sum rule of the dynamic structure factor is strictly satisfied.

For a unitary Fermi superfluid, where ${\xi }_{\mu }={\xi }_{E}=\xi$ due to the scale invariance [4], the latest auxiliary field QMC provides $\xi \simeq 0.372$ [12], which is quite close to the experimental value $\xi =0.376(5)$ [29]. Other field theoretical approaches provide very similar predictions. These include the dimensional - expansion ($\xi =0.377\pm 0.014$ at the next-to-leading order [20] and $\xi =0.360\pm 0.020$ at the next-to-next-to-leading order [23]), self-consistent Luttinger-Ward theory ($\xi =0.36$) [21], and non-self-consistent Gaussian pair fluctuation theory ($\xi =0.401$) [19, 22]. As to the parameter η, its accurate value is to be determined yet. An earlier rf-spectroscopy experiment reports $\eta \simeq 0.44$ [33] and the latest QMC result is $\eta =0.504$ [10, 11]. The field theoretical predictions (i.e., $\eta =0.46$ from both the Luttinger-Ward theory [21] and Gaussian pair-fluctuation theory [19, 22]) also give similar results. In this work, for a unitary Fermi gas we choose the experimental result $\mu =0.376{\varepsilon }_{F}$ for the chemical potential and the QMC prediction ${\rm{\Delta }}=0.5{\varepsilon }_{F}$ for the pairing gap. This leads to $\beta =-0.430$ and $1/\gamma =-0.0767$. It is worth noting that, when $\alpha =1$, our SLDA result reduces that of the standard BdG theory, if we set β and $1/\gamma$ to zero.

Away from the unitary limit, the knowledge on the pairing gap is not complete. We use the predictions of a Gaussian pair fluctuation theory [19, 22] as the inputs, since these theoretical results have already been shown to provide a satisfactory explanation for the experimentally measured chemical potential [28].

If a superfluid Fermi gas is perturbed by a small external potential, usually the number density and anomalous density will fluctuate. Due to the interatomic interactions, the fluctuating densities will feedback and induce an additional perturbation potential. One way to include these fluctuation effects is to use the linear response theory within the RPA [38, 55, 67–70]. The essential idea of RPA is that the induced fluctuation potential is assumed to be a self-generated mean-field potential experienced by quasiparticles, due to the local changes in the number densities ${n}_{\uparrow }({\bf{r}},t)$ and ${n}_{\downarrow }({\bf{r}},t)$, and Cooper-pairs density $\nu ({\bf{r}},t)$ or its complex conjugate ${\nu }^{* }({\bf{r}},t)$. In the following, for convenience, we denote these four densities n_↑, n_↓, ν and ${\nu }^{* }$ as n₁, n₂, n₃ and n₄, respectively.

In the SLDA energy density functional, it is easy to see that the interaction contribution to the functional is given by,

$$\begin{eqnarray}&&{{ \mathcal E }}_{\mathrm{int}}=\beta \displaystyle \frac{3{(3{\pi }^{2})}^{2/3}}{10m}{n}^{5/3}({\bf{r}},t)+\displaystyle \frac{| {\rm{\Delta }}({\bf{r}},t){| }^{2}}{{g}_{\mathrm{eff}}}.\end{eqnarray} \tag{ 11 }$$

The resulting fluctuating potential is simply ${\sum }_{j}{E}_{{ij}}^{I}\delta {n}_{j}$, where [70]

$$\begin{eqnarray}&&{E}_{{ij}}^{I}={\left(\displaystyle \frac{{\partial }^{2}{{ \mathcal E }}_{\mathrm{int}}}{\partial {n}_{i}\partial {n}_{j}}\right)}_{0}\end{eqnarray} \tag{ 12 }$$

and $\delta {n}_{i=\mathrm{1,2},\mathrm{3,4}}$ are the density fluctuations around equilibrium, which are to be determined. The suffix 0 indicates that the derivatives are calculated at equilibrium. Therefore, together with the external potential ${V}_{\mathrm{ext}}^{i}$, the total effective perturbative potential takes the form,

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}^{i}\equiv {V}_{\mathrm{ext}}^{i}+\sum _{j}{E}_{{ij}}^{I}\delta {n}_{j}.\end{eqnarray} \tag{ 13 }$$

Using this effective perturbation, the density fluctuations $\delta {n}_{i}$ can be written down straightforwardly, according to the standard linear response theory,

$$\begin{eqnarray}&&\delta {n}_{i}=\sum _{j}{\chi }_{{ij}}^{0}{V}_{\mathrm{eff}}^{i},\end{eqnarray} \tag{ 14 }$$

where ${\chi }^{0}$ is the bare response function of the quasiparticle reference system described by the SLDA equation (6), which is easy to calculate (see appendix ). By combining equations (13) and (14), we arrive at,

$$\begin{eqnarray}&&\delta {n}_{i}=\sum _{j}{\chi }_{{ij}}{V}_{\mathrm{ext}}^{i},\end{eqnarray} \tag{ 15 }$$

where χ is the RPA response function,

$$\begin{eqnarray}&&\chi ={\chi }^{0}{[1-{\chi }^{0}{E}^{I}]}^{-1}.\end{eqnarray} \tag{ 16 }$$

Once the bare response function ${\chi }^{0}$ and the second order derivative E_ij^I are known, we obtain directly χ. The density response function ${\chi }_{D}$ is a summation of ${\chi }_{{ij}}$ in the density channel: ${\chi }_{D}={\chi }_{11}+{\chi }_{12}+{\chi }_{21}+{\chi }_{22}=2({\chi }_{11}+{\chi }_{12})$. The dynamic structure factor is connected to the imaginary part of the density response function,

$$\begin{eqnarray}&&S({\bf{q}},\omega )=-\displaystyle \frac{1}{\pi }\displaystyle \frac{\mathrm{Im}{\chi }_{D}({\bf{q}},{\rm{i}}{\nu }_{n}\to \omega +{\rm{i}}{0}^{+})}{1-{{\rm{e}}}^{-\omega /T}},\end{eqnarray} \tag{ 17 }$$

with ${\bf{q}}$ and ω being the transferred momentum and energy, respectively.

The RPA on top of the mean-field BdG theory has previously been used to study the dynamic structure factor [67] and collective oscillations [68] of weakly interacting Fermi superfluids. A dynamical mean-field approach, identical to the RPA but based on kinetic equations, was also developed to investigate dynamic and static structure factors and collective modes of strongly interacting Fermi superfluids [38, 55]. Some properties of the density response of unitary Fermi gas for the SLDA has also been studied in [71]. In the following, we examine the improved RPA based on the SLDA theory.

The calculation of the second-order derivative matrix E^I is straightforward. It reads,

$$\begin{eqnarray}{E}^{I}=\left[\begin{array}{cccc}{{ \mathcal I }}_{n}{\varepsilon }_{F}/n & {{ \mathcal I }}_{n}{\varepsilon }_{F}/n & {{ \mathcal I }}_{\nu }{g}_{e\mathrm{ff}} & {{ \mathcal I }}_{\nu }{g}_{e\mathrm{ff}}\\ {{ \mathcal I }}_{n}{\varepsilon }_{F}/n & {{ \mathcal I }}_{n}{\varepsilon }_{F}/n & {{ \mathcal I }}_{\nu }{g}_{e\mathrm{ff}} & {{ \mathcal I }}_{\nu }{g}_{e\mathrm{ff}}\\ {{ \mathcal I }}_{\nu }{g}_{e\mathrm{ff}} & {{ \mathcal I }}_{\nu }{g}_{e\mathrm{ff}} & 0 & {g}_{\mathrm{eff}}\\ {{ \mathcal I }}_{\nu }{g}_{e\mathrm{ff}} & {{ \mathcal I }}_{\nu }{g}_{e\mathrm{ff}} & {g}_{\mathrm{eff}} & 0\end{array}\right],\end{eqnarray} \tag{ 18 }$$

where ${{ \mathcal I }}_{n}$ and ${{ \mathcal I }}_{\nu }$ are two dimensionless variables,

$$\begin{eqnarray*}{{ \mathcal I }}_{n} & = & \displaystyle \frac{2\beta }{3}+\displaystyle \frac{{(3{\pi }^{2})}^{2/3}}{9\gamma }\displaystyle \frac{{\rm{\Delta }}{}^{2}}{{\varepsilon }_{F}^{2}},\\ {{ \mathcal I }}_{\nu } & = & \displaystyle \frac{{(3{\pi }^{2})}^{2/3}}{6\gamma }\displaystyle \frac{{\rm{\Delta }}}{{\varepsilon }_{F}}.\end{eqnarray*}$$

We note the existence of the crossing term ${{ \mathcal I }}_{\nu }$, due to the (implicit) coupling between the number density and the anomalous Cooper-pair density in the interaction energy density functional equation (11). In the unitary limit, in comparison to the BdG-RPA theory, we note also that the matrix element in the number density channel, ${{ \mathcal I }}_{n}{\varepsilon }_{F}/n$, changes from a vanishingly small number (i.e., at the order of ${g}_{\mathrm{eff}}$) to a finite value. The response function of the quasiparticle reference system ${\chi }^{0}$ can be constructed by solving the stationary SLDA equation (6). It is a 4 by 4 matrix. However, as we shown in appendix, only six of all 16 matrix elements are independent:

$$\begin{eqnarray}{\chi }^{0}=\left[\begin{array}{cccc}{\chi }_{11}^{0} & {\chi }_{12}^{0} & {\chi }_{13}^{0} & {\chi }_{14}^{0}\\ {\chi }_{12}^{0} & {\chi }_{11}^{0} & {\chi }_{13}^{0} & {\chi }_{14}^{0}\\ {\chi }_{14}^{0} & {\chi }_{14}^{0} & -{\chi }_{12}^{0} & {\chi }_{34}^{0}\\ {\chi }_{13}^{0} & {\chi }_{13}^{0} & {\chi }_{43}^{0} & -{\chi }_{12}^{0}\end{array}\right],\end{eqnarray} \tag{ 19 }$$

The detailed expressions of the elements ${\chi }_{11}^{0}$, ${\chi }_{12}^{0}$, ${\chi }_{13}^{0}$, ${\chi }_{14}^{0}$, ${\chi }_{34}^{0}$ and ${\chi }_{43}^{0}$ are we show in appendix. By solving the RPA equation (16), we obtain all the matrix elements ${\chi }_{{ij}}$ of the RPA response function χ. The resulting density response function is given by,

$$\begin{eqnarray}{\chi }_{D}=2\left|\begin{array}{ccc}{\chi }_{12}^{0}+{\chi }_{11}^{0} & -{\chi }_{14}^{0}{g}_{\mathrm{eff}} & -{\chi }_{13}^{0}{g}_{\mathrm{eff}}\\ 2{\chi }_{14}^{0} & 1-{\chi }_{34}^{0}{g}_{\mathrm{eff}} & {\chi }_{12}^{0}{g}_{\mathrm{eff}}\\ 2{\chi }_{13}^{0} & {\chi }_{12}^{0}{g}_{\mathrm{eff}} & 1-{\chi }_{43}^{0}{g}_{\mathrm{eff}}\end{array}\right|/| 1-{\chi }^{0}{E}^{I}| .\end{eqnarray} \tag{ 20 }$$

It is well known that the anomalous density correlated functions, ${\chi }_{34}^{0}$ and ${\chi }_{43}^{0}$, are divergent, because of the use of the contact interatomic interaction [68]. Thus, we introduce the regularized functions ${\widetilde{\chi }}_{34}^{0}={\chi }_{34}^{0}-1/{g}_{\mathrm{eff}}$ and ${\widetilde{\chi }}_{43}^{0}={\chi }_{43}^{0}-1/{g}_{\mathrm{eff}}$, with which the density response function now takes the form,

$$\begin{eqnarray}&&{\chi }_{D}=2\displaystyle \frac{({{ \mathcal B }}_{n1}-2{{ \mathcal B }}_{n2})}{| 1-{\chi }^{0}{E}^{I}| /{g}_{\mathrm{eff}}^{2}}.\end{eqnarray} \tag{ 21 }$$

Here,

$$\begin{eqnarray*}{{ \mathcal B }}_{n1} & = & ({\chi }_{11}^{0}+{\chi }_{12}^{0})[{\widetilde{\chi }}_{34}^{0}{\widetilde{\chi }}_{43}^{0}-{({\chi }_{12}^{0})}^{2}],\\ {{ \mathcal B }}_{n2} & = & 2{\chi }_{12}^{0}{\chi }_{13}^{0}{\chi }_{14}^{0}+{({\chi }_{13}^{0})}^{2}{\widetilde{\chi }}_{34}^{0}+{({\chi }_{14}^{0})}^{2}{\widetilde{\chi }}_{43}^{0}.\end{eqnarray*}$$

To obtain the expression of $| 1-{\chi }^{0}{E}^{I}| /{g}_{\mathrm{eff}}^{2}$, it should be noted that ${g}_{\mathrm{eff}}$ is a vanishingly small quantity. Therefore, it is useful to arrange different terms in terms of the powers of ${g}_{{\rm{e}}\mathrm{ff}}$. For instance, for the matrix elements of E_n, ${{ \mathcal I }}_{n}{\varepsilon }_{F}/n$ has the order of ${[{g}_{\mathrm{eff}}]}^{0}$, while ${{ \mathcal I }}_{\nu }{g}_{\mathrm{eff}}$ has the order of ${[{g}_{\mathrm{eff}}]}^{1}$. For the determinant $| 1-{\chi }^{0}{E}^{I}|$, there are no terms at the order of ${ \mathcal O }({g}_{\mathrm{eff}})$ or ${ \mathcal O }(1)$, as anticipated. The order of most terms is ${ \mathcal O }({[{g}_{\mathrm{eff}}]}^{2})$. By collecting those terms, we find that,

$$\begin{eqnarray}&&\displaystyle \frac{| 1-{\chi }^{0}{E}^{I}| }{{g}_{\mathrm{eff}}^{2}}=4{{ \mathcal I }}_{\nu }{{ \mathcal B }}_{\nu 1}+2{{ \mathcal I }}_{\nu }^{2}{{ \mathcal B }}_{\nu 2}-2{{ \mathcal I }}_{n}({{ \mathcal B }}_{n1}-2{{ \mathcal B }}_{n2})+{\widetilde{\chi }}_{34}^{0}{\widetilde{\chi }}_{43}^{0}-{({\chi }_{12}^{0})}^{2},\end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray*}{{ \mathcal B }}_{\nu 1} & = & {\chi }_{12}^{0}{\chi }_{13}^{0}+{\chi }_{12}^{0}{\chi }_{14}^{0}+{\chi }_{13}^{0}{\widetilde{\chi }}_{34}^{0}+{\chi }_{14}^{0}{\widetilde{\chi }}_{43}^{0},\\ {{ \mathcal B }}_{\nu 2} & = & ({\chi }_{11}^{0}+{\chi }_{12}^{0})(2{\chi }_{12}^{0}+{\widetilde{\chi }}_{34}^{0}+{\widetilde{\chi }}_{43}^{0})-2{({\chi }_{13}^{0}-{\chi }_{14}^{0})}^{2},\end{eqnarray*}$$

In the unitary limit, if we set both ${{ \mathcal I }}_{n}$ and ${{ \mathcal I }}_{\nu }$ to zero, $| 1-{\chi }^{0}{E}^{I}| /{g}_{\mathrm{eff}}^{2}$ is just ${\widetilde{\chi }}_{34}^{0}{\widetilde{\chi }}_{43}^{0}-{({\chi }_{12}^{0})}^{2}$, and then we recover the BdG-RPA expression for the density response function [38, 56, 57].

We use equations (21) and (22) to obtain the density response function ${\chi }_{D}$ and then calculate the dynamic structure factor $S({\bf{q}},\omega )$ via the fluctuation-dissipation theorem equation (17). To take the analytic continuation numerically, i.e., ${\rm{i}}{\nu }_{n}\to \omega +{\rm{i}}\delta$, where $\delta ={0}^{+}$, we use a small broadening parameter $\delta ={10}^{-3}{\varepsilon }_{F}$, unless specified elsewhere.

In this section, we present the results for the dynamic structure factor of a unitary Fermi gas at zero temperature within SLDA-RPA, and justify our theory at low transferred momentum $q\leqslant {k}_{F}$ by comparing the resulting static structure factor equation (23) with the latest QMC data [15].

Figure 1 reports a contour plot of $S({\bf{q}},\omega )$ in the momentum range from q = 0 to $q=2{k}_{F}$. Two types of contributions are clearly visible: one is the collective Bogoliubov-Anderson phonon excitations within the energy gap $\omega \lt {E}_{g}=2{\rm{\Delta }}$ [38], which exhibit themselves as a sharp δ-peak in the structure factor spectrum. Right above the energy gap, a much broader distribution emerges, which should be attributed to the fermionic single-particle excitations by breaking Cooper pairs.

Zoom In Zoom Out Reset image size

A close examination of the phonon excitations is shown in figure 2 for a very small transferred momentum $q=0.01{k}_{F}$. For comparison, we also plot the result of the standard BdG-RPA prediction by a red dashed line. It is anticipated that the dispersion of the phonon excitations should follow $\omega ={c}_{s}q$, where c_s is the sound velocity. By fitting the position of the phonon peak as a function of q, we numerically extract a value ${c}_{s}\simeq 0.354{v}_{F}$, which coincides, within the accuracy of our numerical calculations, with the value obtained using the macroscopic definition of the sound speed, ${c}_{s}=\sqrt{(n/m)\partial \mu /\partial n}=\sqrt{{\xi }_{\mu }/3}{v}_{F}$. This value is also consistent with the results determined from the experiments and from the ab initio Monte Carlo calculations. The agreement is not surprising, since the SLDA parameters have been chosen to reproduce the known equation of state and hence the sound speed. It is worth noting that a similar phonon peak is also predicted by the BdG-RPA theory (i.e., using the BdG energy density functional). However, the BdG-RPA theory predicts a sound speed ${c}_{s}\simeq 0.444{v}_{F}$, which is about 30% larger than the above mentioned SLDA-RPA result.

Zoom In Zoom Out Reset image size

At larger transferred momentum, i.e., $q\gtrsim 0.5{k}_{F}$, the single-particle excitations start to make a notable contribution to the dynamic structure factor above the threshold $\omega =2{\rm{\Delta }}={\varepsilon }_{F}$, as shown in figure 3. The sharp rise of the single-particle contribution at $\omega =2{\rm{\Delta }}$ is unlikely to be destroyed by the possible residue interactions between Cooper pairs and unpaired fermions, which is not accounted for in our theory. Therefore, it could serve as a useful feature to experimentally determine the pairing gap in the two-photon Bragg scattering experiments [42]. We also note that, compared with our SLDA-RPA results, the BdG-RPA theory predicts a much weaker response of the single-particle excitations at a larger threshold. This difference between the SLDA- and BdG-RPA predictions could be easily resolved experimentally.

Zoom In Zoom Out Reset image size

A test of the accuracy of the theory can be obtained by looking at the static structure factor

$$\begin{eqnarray}&&S({\bf{q}})=\int {\rm{d}}\omega S({\bf{q}},\omega )\end{eqnarray} \tag{ 23 }$$

for which QMC results are available [15, 55]. The comparison of our SLDA-RPA predictions with the latest diffusion Monte Carlo data [15] is shown in figure 4, together with the predictions of BdG-RPA. The excellent agreement between SLDA-RPA and QMC at $q\leqslant {k}_{F}$ is non-trivial and suggests that our theory can be quantitatively reliable at small momentum transfer. Above the Fermi momentum, instead, there are significant deviations. It is worth noticing that the BdG-RPA theory gives results closer to QMC at large momentum transfer, where the physics is dominated by single-particle excitations and where BdG-RPA theory is known to work well [56].

Zoom In Zoom Out Reset image size

In figure 5, we show the dynamic structure factor at the momentum $q=4{k}_{F}$. At such a large momentum, one can still separately resolve the bosonic Cooper-pair excitations (i.e., a molecular peak structure at $\omega ={q}^{2}/4m=8{\varepsilon }_{F}$) and fermionic single-particle excitations (i.e., the broader distribution at $\omega ={q}^{2}/2m=16{\varepsilon }_{F}$). Compared with the BdG-RPA result, our SLDA-RPA theory predicts a much smaller molecular peak. This is understandable, since the SLDA theory is effectively a low-energy theory and hence becomes less efficient at $\omega \gg {\varepsilon }_{F}$. We note that, experimentally, there is a finite energy resolution in the measurement of the dynamic structure factor [56]. The notable difference in the predictions for the molecular peak will be easily smeared out by the finite energy resolution. As a result, the SLDA-RPA approach may predict nearly the same line shape as the BdG-RPA theory. The difference in the line shape is characterized by the relative difference in the static structure factor, which is about 5%. In the sense of predicting the experimental line shape for the dynamic structure factor, we may argue that the SLDA-RPA is semi-quantitatively valid at large transferred momentum $q\gt {k}_{F}$.

Zoom In Zoom Out Reset image size

It should also be noted that an independent check of the SLDA-RPA theory is provided by the f-sum rule [51]

$$\begin{eqnarray}&&\int {\rm{d}}\omega \omega S({\bf{q}},\omega )=\displaystyle \frac{N{{\bf{q}}}^{2}}{2m},\end{eqnarray} \tag{ 24 }$$

which should be satisfied. We have numerically checked that our SLDA-RPA calculations obey this sum-rule within 1% relative accuracy.

In this section, we apply the SLDA-RPA theory to determine the dynamic structure factor at the whole BCS-BEC crossover, by using the zero-temperature chemical potential and pairing gap calculated from a Gaussian pair fluctuation theory [19] as the inputs. The energy density functional equation (1)—obtained under the scale invariance assumption—is supposed to work well slightly away from the unitary limit.

Figure 6 reports the dynamic structure factor at the BCS-BEC crossover at two different transferred momenta $q=0.5{k}_{F}$ (a) and $q={k}_{F}$ (b). On the BCS side, the single-particle contributions become significant, as one may anticipate. Furthermore, at $q={k}_{F}$ and $1/({k}_{F}a)=-0.4$, where the bosonic peak position ${\omega }_{B}\sim {c}_{s}q$ is close to the two-particle scattering threshold $2{\rm{\Delta }}$, there is a strong overlap between the phonon and single-particle contributions, leading to an interesting peak-dip-bump structure. When the system crosses over to the BEC limit with increasing $1/({k}_{F}a)$, the phonon peak moves to the low energy, due to the decreasing sound velocity. The single-particle contributions get suppressed very quickly. In particular, at $q=0.5{k}_{F}$, the broader single-particle distribution can be barely seen on the BEC side with $1/({k}_{F}a)\gt 0$.

Zoom In Zoom Out Reset image size

Apparently, the experimental determination of the phonon peaks can be ideally used to measure the sound velocity across the BCS-BEC crossover. The measurement of the broader single-particle contributions may also be useful to determine the pairing gap on the BCS side.

In summary, we have developed a random phase approximation theory for calculating the dynamic structure factor of a strongly interacting Fermi gas at unitarity and in the BCS-BEC crossover, within the framework of a DFT approach [60, 66]. The theory is expected to be quantitatively reliable at low transferred momentum (i.e., $q\lt {k}_{F}$) and at low temperature (i.e., $T\ll {T}_{c}$), where the predicted static structure factor agrees excellently well with the result of the latest ab initio diffusion QMC [15]. Therefore, our theory is useful to understand the dynamic structure factor in the previously un-explored territory of low transferred momentum, as schematically illustrated in figure 7 by a red rectangle. A stringent test of the applicability of our theory could be obtained by comparing our predictions with the results of on-going experiments [42].

Zoom In Zoom Out Reset image size

We are grateful to Chris Vale, Sandro Stringari, Aurel Bulgac, Michael McNeil Forbes and Lianyi He for fruitful discussions, and Stefano Giorgini and Stefano Gandolfi for sharing their QMC data. PZ is indebted to the BEC Center at Trento for hospitality when this work started. This work was supported by the ARC Discovery Projects: FT130100815 and DP140103231 (HH), DP140100637 , and FT140100003 (XJL). RS acknowledges support from DAE, Government of India. The work is also supported by Provincia Autonoma di Trento (FD). Correspondence should be addressed to PZ at phy.zoupeng@gmail.com.

In this appendix, we discuss how to calculate the response function ${\chi }^{0}$, by solving the stationary SLDA equation. The existence of four different densities means that there will be 16 correlation functions in ${\chi }^{0}$:

$$\begin{eqnarray}{\chi }^{0}\equiv \left[\begin{array}{cccc}\langle {n}_{1}{n}_{1}{\rangle }_{0} & \langle {n}_{1}{n}_{2}{\rangle }_{0} & \langle {n}_{1}{n}_{3}{\rangle }_{0} & \langle {n}_{1}{n}_{4}{\rangle }_{0}\\ \langle {n}_{2}{n}_{1}{\rangle }_{0} & \langle {n}_{2}{n}_{2}{\rangle }_{0} & \langle {n}_{2}{n}_{3}{\rangle }_{0} & \langle {n}_{2}{n}_{4}{\rangle }_{0}\\ \langle {n}_{3}{n}_{1}{\rangle }_{0} & \langle {n}_{3}{n}_{2}{\rangle }_{0} & \langle {n}_{3}{n}_{3}{\rangle }_{0} & \langle {n}_{3}{n}_{4}{\rangle }_{0}\\ \langle {n}_{4}{n}_{1}{\rangle }_{0} & \langle {n}_{4}{n}_{2}{\rangle }_{0} & \langle {n}_{4}{n}_{3}{\rangle }_{0} & \langle {n}_{4}{n}_{4}{\rangle }_{0}\end{array}\right],\end{eqnarray} \tag{ A1 }$$

where the abbreviation ${\chi }_{{ij}}^{0}=\langle {n}_{i}{n}_{j}{\rangle }_{0}$ is used. The derivation of these matrix elements is cumbersome. We show here, as an example, the derivation of ${\chi }_{\uparrow \uparrow }^{0}\equiv {\chi }_{11}^{0}$. According to the Wick theorem, and following the BCS theory, which assume that only propagators—like $\langle {{\rm{\Psi }}}_{\uparrow }^{\dagger }{{\rm{\Psi }}}_{\uparrow }\rangle$, $\langle {{\rm{\Psi }}}_{\downarrow }^{\dagger }{{\rm{\Psi }}}_{\downarrow }\rangle$, $\langle {{\rm{\Psi }}}_{\downarrow }{{\rm{\Psi }}}_{\uparrow }\rangle$ and $\langle {{\rm{\Psi }}}_{\uparrow }^{\dagger }{{\rm{\Psi }}}_{\downarrow }^{\dagger }\rangle$ —have a non-zero value, the imaginary-time Green's function ${\chi }_{11}^{0}({\bf{r}},{\bf{r}}^{\prime} ,\tau )$ $\equiv \,-\langle {T}_{\tau }[{\hat{n}}_{1}({\bf{r}},\tau ){\hat{n}}_{1}({\bf{r}}^{\prime} ,0)]\rangle$ can be written as

$$\begin{eqnarray}&&{\chi }_{11}^{0}=-\langle {{\rm{\Psi }}}_{\uparrow }^{\dagger }({\bf{r}},\tau ){{\rm{\Psi }}}_{\uparrow }({\bf{r}}^{\prime} ,0)\rangle \langle {{\rm{\Psi }}}_{\uparrow }({\bf{r}},\tau ){{\rm{\Psi }}}_{\uparrow }^{\dagger }({\bf{r}}^{\prime} ,0)\rangle ,\end{eqnarray} \tag{ A2 }$$

where τ is the imaginary time and we assume $\tau \gt 0$. By using the Bogoliubov transformations

$$\begin{eqnarray}{{\rm{\Psi }}}_{\uparrow } & = & \displaystyle \sum _{j}[{u}_{j\uparrow }({\bf{r}}){c}_{j\uparrow }{{\rm{e}}}^{-{{\rm{i}}{E}}_{j\uparrow }t}+{v}_{j\downarrow }^{* }({\bf{r}}){c}_{j\downarrow }^{\dagger }{{\rm{e}}}^{{{\rm{i}}{E}}_{j\downarrow }t}],\\ {{\rm{\Psi }}}_{\downarrow }^{\dagger } & = & \displaystyle \sum _{j}[{u}_{j\downarrow }^{* }({\bf{r}}){c}_{j\downarrow }^{\dagger }{{\rm{e}}}^{{{\rm{i}}{E}}_{j\downarrow }t}-{v}_{j\uparrow }({\bf{r}}){c}_{j\uparrow }{{\rm{e}}}^{-{{\rm{i}}{E}}_{j\uparrow }t}],\end{eqnarray} \tag{ A3 }$$

for the field operators ${{\rm{\Psi }}}_{\sigma }$ and ${{\rm{\Psi }}}_{\sigma }^{\dagger }$, one finds

$$\begin{eqnarray}&&{\chi }_{11}^{0}({\bf{r}},{\bf{r}}^{\prime} ,\tau )=-\displaystyle \sum _{i,j}{u}_{i}^{* }({\bf{r}}){u}_{i}({\bf{r}}^{\prime} ){u}_{j}({\bf{r}}){u}_{j}^{* }({\bf{r}}^{\prime} )f({E}_{i})f(-{E}_{j}){{\rm{e}}}^{({E}_{i}-{E}_{j})\tau }.\end{eqnarray} \tag{ A4 }$$

Here we use $\langle {c}_{i}^{\dagger }{c}_{j}\rangle =f({E}_{i}){\delta }_{{ij}}$ and $\langle {c}_{i}{c}_{j}^{\dagger }\rangle =f(-{E}_{i}){\delta }_{{ij}}$, and $f(x)=1/({{\rm{e}}}^{x/T}+1)$ is the Fermi distribution function of quasiparticles. The spin index has been removed owing to the existence of a one-to-one correspondence between the solutions of spin-up and spin-down energy levels. By taking the Fourier transformation in the imaginary time, ${\chi }_{11}^{0}({\bf{r}},{\bf{r}}^{\prime} ,{\rm{i}}{\nu }_{n})={\int }_{0}^{\beta }{\rm{d}}\tau {{\rm{e}}}^{{\rm{i}}{\nu }_{n}\tau }{\chi }_{11}^{0}({\bf{r}},{\bf{r}}^{\prime} ,\tau )$, where ${\nu }_{n}=2n\pi {k}_{B}T$ is the bosonic Matsubara frequency, one obtains,

$$\begin{eqnarray}&&{\chi }_{11}^{0}({\bf{r}},{\bf{r}}^{\prime} ,{\rm{i}}{\nu }_{n})=\displaystyle \sum _{{i},j}{u}_{{\rm{i}}}^{* }({\bf{r}}){u}_{{i}}({\bf{r}}^{\prime} ){u}_{j}({\bf{r}}){u}_{j}^{* }({\bf{r}}^{\prime} )\displaystyle \frac{f({E}_{{i}})-f({E}_{j})}{{\rm{i}}{\nu }_{n}+({E}_{{i}}-{E}_{j})}.\end{eqnarray} \tag{ A5 }$$

For the homogeneous gas, a set of plane wave functions can be used to expand the eigenfunctions u_i in the form ${u}_{i}({\bf{r}})\to {u}_{k}{{\rm{e}}}^{{\rm{i}}{\bf{kr}}}$. By defining the transferring momentum ${\bf{p}}={\bf{k}}^{\prime} -{\bf{k}}$ and the relative coordinate $\delta {\bf{r}}={\bf{r}}-{\bf{r}}^{\prime}$, then

$$\begin{eqnarray}&&{\chi }_{11}^{0}(\delta {\bf{r}},{\rm{i}}{\nu }_{n})=\displaystyle \sum _{{\bf{k}},{\bf{p}}}{| {u}_{{\bf{k}}}| }^{2}{| {u}_{{\bf{k}}+{\bf{p}}}| }^{2}{{\rm{e}}}^{{\rm{i}}{\bf{p}}\delta {\bf{r}}}\displaystyle \frac{f({E}_{{\bf{k}}})-f({E}_{{\bf{k}}+{\bf{p}}})}{{\rm{i}}{\nu }_{n}+({E}_{{\bf{k}}}-{E}_{{\bf{k}}+{\bf{p}}})}.\end{eqnarray} \tag{ A6 }$$

By taking the Fourier transformation of the relative coordinate, ${\chi }_{11}^{0}(q,{\omega }_{n})=\int {\rm{d}}\delta {\bf{r}}{\chi }_{11}^{0}(\delta {\bf{r}},{\rm{i}}{\omega }_{n}){{\rm{e}}}^{-{\rm{i}}{\bf{q}}\delta {\bf{r}}}$, we find that,

$$\begin{eqnarray}&&{\chi }_{11}^{0}({\bf{q}},{\rm{i}}{\nu }_{n})=\displaystyle \sum _{{\bf{k}}}{| {u}_{{\bf{k}}}| }^{2}{| {u}_{{\bf{k}}+{\bf{p}}}| }^{2}\displaystyle \frac{f({E}_{{\bf{k}}})-f({E}_{{\bf{k}}+{\bf{p}}})}{{\rm{i}}{\nu }_{n}+({E}_{{\bf{k}}}-{E}_{{\bf{k}}+{\bf{p}}})}.\end{eqnarray} \tag{ A7 }$$

Using the expressions for ${u}_{{\bf{k}}}$ and ${u}_{{\bf{k}}+{\bf{p}}}$, at zero temperature we obtain,

$$\begin{eqnarray}&&{\chi }_{11}^{0}({\bf{q}},{\rm{i}}{\nu }_{n})=\displaystyle \sum _{{\bf{k}}}\displaystyle \frac{1}{2}\left(1-\displaystyle \frac{{\xi }_{{\bf{k}}}{\xi }_{{\bf{k}}+{\bf{q}}}}{{E}_{{\bf{k}}}{E}_{{\bf{k}}+{\bf{q}}}}\right)\displaystyle \frac{{E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}}}{{({\rm{i}}{\nu }_{n})}^{2}-{({E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}})}^{2}}.\end{eqnarray} \tag{ A8 }$$

Through a similar process, we can derive the other 15 matrix elements of ${\chi }^{0}$. In fact, after checking their expressions, only six of them are independent. The remaining expressions are simply related to each other by, for example, the replacement ${\bf{k}}\to -{\bf{k}}-{\bf{q}}$. In the following, we list the other five expressions for ${\chi }_{12}^{0}$, ${\chi }_{13}^{0}$, ${\chi }_{14}^{0}$, ${\chi }_{34}^{0}$ and ${\chi }_{43}^{0}$ at zero temperature:

$$\begin{eqnarray}&&{\chi }_{12}^{0}=\displaystyle \sum _{{\bf{k}}}\displaystyle \frac{1}{2}\displaystyle \frac{{{\rm{\Delta }}}^{2}}{{E}_{{\bf{k}}}{E}_{{\bf{k}}+{\bf{q}}}}\displaystyle \frac{{E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}}}{{({\rm{i}}{\nu }_{n})}^{2}-{({E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}})}^{2}},\end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray}&&{\chi }_{13}^{0}=\displaystyle \sum _{{\bf{k}}}\displaystyle \frac{{\rm{\Delta }}}{4}\left[\displaystyle \frac{({\xi }_{{\bf{k}}}+{\xi }_{{\bf{k}}+{\bf{q}}})}{{E}_{{\bf{k}}}{E}_{{\bf{k}}+{\bf{q}}}}\displaystyle \frac{{E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}}}{{({\rm{i}}{\nu }_{n})}^{2}-{({E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}})}^{2}}-\left(\displaystyle \frac{1}{{E}_{{\bf{k}}}}+\displaystyle \frac{1}{{E}_{{\bf{k}}+{\bf{q}}}}\right)\displaystyle \frac{{\rm{i}}{\nu }_{n}}{{({\rm{i}}{\nu }_{n})}^{2}-{({E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}})}^{2}}\right],\end{eqnarray} \tag{ A10 }$$

$$\begin{eqnarray}&&{\chi }_{14}^{0}=\displaystyle \sum _{{\bf{k}}}\displaystyle \frac{{\rm{\Delta }}}{4}\left[\displaystyle \frac{({\xi }_{{\bf{k}}}+{\xi }_{{\bf{k}}+{\bf{q}}})}{{E}_{{\bf{k}}}{E}_{{\bf{k}}+{\bf{q}}}}\displaystyle \frac{{E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}}}{{({\rm{i}}{\nu }_{n})}^{2}-{({E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}})}^{2}}+\left(\displaystyle \frac{1}{{E}_{{\bf{k}}}}+\displaystyle \frac{1}{{E}_{{\bf{k}}+{\bf{q}}}}\right)\displaystyle \frac{{\rm{i}}{\nu }_{n}}{{({\rm{i}}{\nu }_{n})}^{2}-{({E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}})}^{2}}\right],\end{eqnarray} \tag{ A11 }$$

$$\begin{eqnarray}&&{\chi }_{34}^{0}=\displaystyle \sum _{{\bf{k}}}\displaystyle \frac{1}{2}\left[\left(1+\displaystyle \frac{{\xi }_{{\bf{k}}}}{{E}_{{\bf{k}}}}\displaystyle \frac{{\xi }_{{\bf{k}}+{\bf{q}}}}{{E}_{{\bf{k}}+{\bf{q}}}}\right)\displaystyle \frac{{E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}}}{{({\rm{i}}{\nu }_{n})}^{2}-{({E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}})}^{2}}+\left(\displaystyle \frac{{\xi }_{{\bf{k}}}}{{E}_{{\bf{k}}}}+\displaystyle \frac{{\xi }_{{\bf{k}}+{\bf{q}}}}{{E}_{{\bf{k}}+{\bf{q}}}}\right)\displaystyle \frac{{\rm{i}}{\nu }_{n}}{{({\rm{i}}{\nu }_{n})}^{2}-{({E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}})}^{2}}\right],\end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray}&&{\chi }_{43}^{0}=\displaystyle \sum _{{\bf{k}}}\displaystyle \frac{1}{2}\left[\left(1+\displaystyle \frac{{\xi }_{{\bf{k}}}}{{E}_{{\bf{k}}}}\displaystyle \frac{{\xi }_{{\bf{k}}+{\bf{q}}}}{{E}_{{\bf{k}}+{\bf{q}}}}\right)\displaystyle \frac{{E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}}}{{({\rm{i}}{\nu }_{n})}^{2}-{({E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}})}^{2}}-\left(\displaystyle \frac{{\xi }_{{\bf{k}}}}{{E}_{{\bf{k}}}}+\displaystyle \frac{{\xi }_{{\bf{k}}+{\bf{q}}}}{{E}_{{\bf{k}}+{\bf{q}}}}\right)\displaystyle \frac{{\rm{i}}{\nu }_{n}}{{({\rm{i}}{\nu }_{n})}^{2}-{({E}_{{\bf{k}}}+{E}_{{\bf{k}}+{\bf{q}}})}^{2}}\right].\end{eqnarray} \tag{ A13 }$$

We note that, ${\chi }_{34}^{0}$ and ${\chi }_{43}^{0}$ should be regularized in order to remove the ultraviolet divergence.