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

We theoretically investigate the dynamic structure factor of a strongly interacting Fermi gas at the crossover from Bardeen-Cooper-Schrieffer superfluids to Bose-Einstein condensates, by developing an improved random phase approximation within the framework of a density functional theory - the so-called superfluid local density approximation. Compared with the previous random-phase-approximation studies based on the standard Bogoliubov-de Gennes equations, the use of the density functional theory greatly improves the accuracy of the equation of state at the crossover, and leads to a better description of both collective Bogoliubov-Anderson-Goldstone phonon mode and single-particle fermionic excitations at small transferred momentum. Near unitarity, where the s-wave scattering length diverges, we show that the single-particle excitations start to significantly contribute to the spectrum of dynamic structure factor once the frequency is above a threshold of the energy gap at $2\Delta$. The sharp rise in the spectrum at this threshold can be utilized to measure the pairing gap $\Delta$. Together with the sound velocity determined from the phonon branch, the dynamic structure factor provides us some key information of the crossover Fermi superfluid. Our predictions could be examined in experiments with $^{6}$Li or $^{40}$K atoms using Bragg spectroscopy.

As an important fingerprint of quantum gases in some certain states, dynamic structure factor contains rich information of properties of a many-body system [41]. 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 [32,33]. Experimentally, the dynamic structure factor can be measured via two-photon Bragg scattering technique [37], at both low and finite temperatures [39]. 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  [44,45], Tan relation [46,47], two-fluid hydrodynamics [48], diagrammatic strong-coupling approach [50,51], and BdG-RPA [36,[52][53][54]. The applicable conditions for the transferred momentum q and temperature T , under which each theory is quantitatively useful, are indicated. Here kF = (3π 2 n) 1/3 , εF = k 2 F /(2m) = (3π 2 n) 2/3 /(2m), and TF are the Fermi momentum, energy, and temperature, respectively. Tc is the superfluid transition temperature.
in certain limiting cases [38] (see Table I). For example, at high temperature, as the fugacity is a small parameter, a quantum cluster expansion has been proven to be an efficient method [42,43], and has been used to calculate the dynamic structure factor [44,45]. In the limits of both large momentum and high frequency, asymptotically exact Tan relations have been derived to describe the high-frequency tails [46,47]. On the other hand, in the limit of long wavelength or small momentum, the phenomenological two-fluid hydrodynamic theory may provide a useful description [48].
A general theoretical framework of the dynamic structure factor, valid at arbitrary temperature and momentum, can be developed by using the diagrammatic technique [49,50] or functional path integral approach [51], 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 [51]. However, their numerical calculations turn out to be extremely difficult, except in the limit of zero transferred momentum [50]. A more commonly used approach is the random-phase approximation (RPA) on top of the mean-field Bogoliubov-de Gennes (BdG) theory [36,[52][53][54]. 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 [53], 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 ∼ 5k F , where k F is the Fermi momentum). At small transferred momentum, i.e., q 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 [55][56][57], within the framework of density functional theory [59][60][61], 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 [57] and real-time dynamics [62,63] 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,52]. A more stringent test can be obtained in the near future by comparing our predictions with the experimental data [40], without any adjustable parameters.
Our paper is organized as follows. In the next section (Sec. II), we introduce the SLDA theory. In Sec. III, we review the main idea of RPA. The expression for the dynamic structure factor is derived in Sec. IV. In Sec. V and Sec. VI, we present our main results of dynamic structure factor in the unitary limit and the crossover regime, respectively. Finally, Sec. VII is devoted to conclusions and outlooks. For convenience, we set = k B = 1 in the following discussions.

II. SUPERFLUID LOCAL DENSITY APPROXIMATION
The density functional theory (DFT) developed by Hohenberg and Kohn [59], together with the local density approximation (LDA) by Kohn and Sham [60], 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 [61] 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 [55][56][57]. This SLDA originates from a similar DFT previously used in the context of nuclear physics [56,63].
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 ν(r, t), in addition to the number density n(r, t).
The energy density functional E[τ (r, t), n(r, t), ν(r, t)] of the system must include the kinetic density τ (r, t), number density n(r, t), and also the anomalous density ν(r, t) [56,57]: where the kinetic density τ , number density n and anomalous density ν are given by, and u k (r, t) and v k (r, t) are the Bogoliubov quasiparticle wavefunctions with 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 [57,63].
In the unitary limit at zero temperature, the simple form of the energy density functional Eq. (1) is inspired by the dimensional analysis: the first and third terms are the unique combination required by the renormalizablity of the theory [55]; 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 [57] and dynamics [62,63] of a unitary Fermi gas at zero temperature. It is reasonable to assume that the energy density functional Eq. (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 [55]. After regularization, the energy density functional with regularized kinetic density τ c (r, t) and anomalous density ν c (r, t) takes the following form [57], where the effective coupling constant g eff is given by We note that, g eff scales to zero once the cut-off momentum Λ runs to infinity. The order parameter ∆(r, t) is related to the anomalous density ν by The stationary SLDA equations for the quasiparticle wave functions are obtained by the standard functional minimization with respect to the variations u k and v k . One obtains with a single quasiparticle Hamiltonian and the chemical potential µ. By requiring that a homogeneous Fermi gas of the number density n = N/V = k 3 F /(3π 2 ) has an energy per particle E/N = (3/5)ξ E ε F , a chemical potential µ = ξ µ ε F , and a pairing order parameter ∆ = ηε F at zero temperature, one can determine the value of dimensionless parameters α, β and γ in Eq. (3) through the following equations, which are independent on the cut-off momentum (i.e., Λ → ∞): and where In these three constraint equations, ξ µ , η and ξ E are the three inputs, whose value can be reliably determined by using QMC simulations [11][12][13] or from the experimental measurements [26,27,30]. The parameter α can be determined using the single particle dispersion, near the unitary limit, typically the parameter α is very close to 1 [57,58], indicating that the effective mass only differs slightly from the bare atomic mass m. For simplicity, throughout the work we take α = 1 and use the density equation Eq. (8) and the gap equation Eq. (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 ξ µ = ξ E = ξ due to the scale invariance [4], the latest auxiliary field QMC provides ξ ≃ 0.372 [12], which is quite close to the experimental value ξ = 0.376(5) [27]. As to the parameter η, its accurate value is to be determined yet. An earlier rf-spectroscopy experiment reports η ≃ 0.44 [31] and the latest QMC result is η = 0.504 [10,11]. In this work, for a unitary Fermi gas we choose the experimental result µ = 0.376ε F for the chemical potential and the QMC prediction ∆ = 0.5ε F for the pairing gap. This leads to β = −0.430 and 1/γ = −0.0767. It is worth noting that, when α = 1, our SLDA result reduces that of the standard BdG theory, if we set β and 1/γ 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,21] as the inputs, since these theoretical results have already been shown to provide a satisfactory explanation for the experimentally measured chemical potential [26].

III. RANDOM PHASE APPROXIMATION
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 [36,52,[64][65][66][67]. 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 ↑ (r, t) and n ↓ (r, t), and Cooper-pairs density ν(r, t) or its complex conjugate ν * (r, t). In the following, for convenience, we denote these four densities n ↑ , n ↓ , ν and ν * as n 1 , n 2 , n 3 and n 4 , respectively.
In the SLDA energy density functional, it is easy to see that the interaction contribution to the functional is given by, The resulting fluctuating potential is simply j E I ij δn j , where [67] E I ij = and δn i=1,2,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 i ext , the total effective perturbative potential takes the form, Using this effective perturbation, the density fluctuations δn i can be written down straightforwardly, according to the standard linear response theory, where χ 0 is the bare response function of the quasiparticle reference system described by the SLDA equation (6), which is easy to calculate (see Appendix A). By combining Eqs. (13) and (14), we arrive at, where χ is the RPA response function, Once the bare response function χ 0 and the second order derivative E I ij are known, we obtain directly χ. The density response function χ D is a summation of χ ij in the density channel: χ D = χ 11 + χ 12 + χ 21 + χ 22 = 2(χ 11 + χ 12 ). The dynamic structure factor is connected to the imaginary part of the density response function, with 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 [64] and collective oscillations [65] 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 [36,52]. Some properties of the density response of unitary Fermi gas for the SLDA has also been studied in [68]. In the following, we examine the improved RPA based on the SLDA theory.

IV. DYNAMIC STRUCTURE FACTOR IN SLDA THEORY
The calculation of the second-order derivative matrix E I is straightforward. It reads, where I n and I ν are two dimensionless variables, We note the existence of the crossing term I ν , due to the (implicit) coupling between the number density and the anomalous Cooper-pair density in the interaction energy density functional Eq. (11). In the unitary limit, in comparison to the BdG-RPA theory, we note also that the matrix element in the number density channel, I n ε F /n, changes from a vanishingly small number (i.e., at the order of g eff ) to a finite value. The response function of the quasiparticle reference system χ 0 can be constructed by solving the stationary SLDA equation (6). It is a 4 by 4 matrix. However, as we shown in Appendix A, only six of all 16 matrix elements are independent: The detailed expressions of the elements χ 0 11 , χ 0 12 , χ 0 13 , χ 0 14 , χ 0 34 and χ 0 43 are we show in Appendix A. By solving the RPA equation (16), we obtain all the matrix elements χ ij of the RPA response function χ. The resulting density response function is given by, It is well known that the anomalous density correlated functions, χ 0 34 and χ 0 43 , are divergent, because of the use of the contact interatomic interaction [65]. Thus, we introduce the regularized functions χ 0 34 = χ 0 34 − 1/g eff and χ 0 43 = χ 0 43 − 1/g eff , with which the density response function now takes the form, Here, To obtain the expression of 1 − χ 0 E I /g 2 eff , it should be noted that g eff is a vanishingly small quantity. Therefore, it is useful to arrange different terms in terms of the powers of g eff . For instance, for the matrix elements of E n , I n ε F /n has the order of [g eff ] 0 , while I ν g eff has the order of [g eff ] 1 . For the determinant 1 − χ 0 E I , there are no terms at the order of O(g eff ) or O(1), as anticipated. The order of most terms is O([g eff ] 2 ). By collecting those terms, we find that, where B ν1 = χ 0 12 χ 0 13 + χ 0 12 χ 0 14 + χ 0 In the unitary limit, if we set both I n and I ν to zero, 1 − χ 0 E I /g 2 eff is just χ 0 34 χ 0 43 − χ 0 12 2 , and then we recover the BdG-RPA expression for the density response function [36,53,54]. We use Eqs. (21) and (22) to obtain the density response function χ D and then calculate the dynamic structure factor S(q, ω) via the fluctuation-dissipation theorem Eq. (17). To take the analytic continuation numerically, i.e., iν n → ω + iδ, where δ = 0 + , we use a small broadening parameter δ = 10 −3 ε F , unless specified elsewhere.

V. DYNAMIC STRUCTURE FACTOR OF A UNITARY FERMI SUPERFLUID
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 ≤ k F by comparing the resulting static structure factor Eq. (23) with the latest QMC data [15]. Fig. 1 reports a contour plot of S(q, ω) in the momentum range from q = 0 to q = 2k F . Two types of contributions are clearly visible: one is the collective Bogoliubov-Anderson phonon excitations within the energy gap ω < E g = 2∆ Here, to better represent the distribution of a delta function at small q, a broadening width δ = 10 −4 εF has been used. The dynamic structure factor is measured in units of N/εF . [36], 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. A close examination of the phonon excitations is shown in Fig. 2 for a very small transferred momentum q = 0.01k 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 ω = 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 ≃ 0.354v F , which coincides, within the accuracy of our numerical calculations, with the value obtained using the macroscopic definition of the sound speed, c s = (n/m)∂µ/∂n = ξ µ /3v 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,   [15] and the BdG-RPA prediction (red dashed line). Our SLDA-RPA theory is expected to be quantitatively reliable at q ≤ kF , as highlighted by the yellow area.
the BdG-RPA theory predicts a sound speed c s ≃ 0.444v F , which is about 30% larger than the above mentioned SLDA-RPA result. At larger transferred momentum, i.e., q 0.5k F , the single-particle excitations start to make a notable contribution to the dynamic structure factor above the threshold ω = 2∆ = ε F , as shown in Fig. 3. The sharp rise of the singleparticle contribution at ω = 2∆ 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 [40]. 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.
A test of the accuracy of the theory can be obtained by looking at the static structure factor for which for which QMC results are available [15,52]. The comparison of our SLDA-RPA predictions with the latest diffusion Monte Carlo data [15] is shown in Fig. 4, together with the predictions of BdG-RPA. The excellent agreement between SLDA-RPA and QMC at q ≤ 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 [53]. In Fig. 5, we show the dynamic structure factor at the momentum q = 4k F . At such a large momentum, one can still separately resolve the bosonic Cooper-pair excitations (i.e., a molecular peak structure at ω = q 2 /4m = 8ε F ) and fermionic single-particle excitations (i.e., the broader distribution at ω = q 2 /2m = 16ε 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 ω ≫ ε F . We note that, experimentally, there is a finite energy resolution in the measurement of the dynamic structure factor [53]. 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 > k F . It should also be noted that an independent check of the SLDA-RPA theory is provided by the f -sum rule [49] dωωS(q, ω) = N q 2 2m , which should be satisfied. We have numerically checked that our SLDA-RPA calculations obey this sum-rule within 1% relative accuracy.

VI. DYNAMIC STRUCTURE FACTOR AT THE BCS-BEC CROSSOVER
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 Eq. (1) -obtained under the scale invariance assumption -is supposed to work well slightly away from the unitary limit. Fig. 6 reports the dynamic structure factor at the BCS-BEC crossover at two different transferred momenta q = 0.5k 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 ω B ∼ c s q is close to the twoparticle scattering threshold 2∆, 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.5k F , the broader single-particle distribution can be barely seen on the BEC side with 1/(k F a) > 0.
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.

VII. CONCLUSIONS
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 density functional theory approach [57,63]. The theory is expected to be quantitatively reliable at low transferred momentum (i.e., q < k F ) and at low temperature (i.e., T ≪ T c ), where the predicted static structure factor agrees excellently well with the result of the latest ab-initio diffusion quantum Monte Carlo [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 Fig. 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 [40].  [44,45], BdG-RPA theory [36,52,53] and diagrammatic approach [50,51]. The applicable parameter space of our SLDA-RPA theory is enclosed by the red dashed line at small transferred momentum q ≤ kF and at low temperature T ≪ Tc. The two-photon Bragg scattering experiment has so far been carried out at q ∼ 0.5kF [40] and q ≥ 3kF [37,39]. The dashed borders of the domains should not be considered as sharp boundaries, but just as an illustrative guide.
where the abbreviation χ 0 ij = n i n j 0 is used. The derivation of these matrix elements is cumbersome. We show here, as an example, the derivation of χ 0 ↑↑ ≡ χ 0 11 . According to the Wick theorem, and following the BCS theory, which assume that only propagators -like Ψ † ↑ Ψ ↑ , Ψ † ↓ Ψ ↓ , Ψ ↓ Ψ ↑ and Ψ † ↑ Ψ † ↓ -have a non-zero value, the imaginary-time Green's function χ 0 11 (r, r ′ , τ ) ≡ − T τ [n 1 (r, τ )n 1 (r ′ , 0)] can be written as where τ is the imaginary time and we assume τ > 0. By using the Bogoliubov transformations Ψ ↑ = j u j↑ (r)c j↑ e −iE j↑ t + v * j↓ (r)c † j↓ e iE j↓ t , Using the expressions for u k and u k+p , at zero temperature we obtain, Through a similar process, we can derive the other 15 matrix elements of χ 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 k → −k − q. In the following, we list the other five expressions for χ 0 12 , χ 0 13 , χ 0 14 , χ 0 34 and χ 0 43 at zero temperature: We note that, χ 0 34 and χ 0 43 should be regularized in order to remove the ultraviolet divergence.