Dynamic density structure factor of a unitary Fermi gas at finite temperature

We present a theoretical investigation of the dynamic density structure factor of a strongly interacting Fermi gas near a Feshbach resonance at finite temperature. The study is based on a gauge invariant linear response theory. The theory is consistent with a diagrammatic approach for the equilibrium state taking into account the pair fluctuation effects and respects some important restrictions like the $f$-sum rule. Our numerical results show that the dynamic density structure factor at large incoming momentum and at half recoil frequency has a qualitatively similar behavior as the order parameter, which can signify the appearance of the condensate. This qualitatively agrees with the recent Bragg spectroscopy experiment results. We also present the results at small incoming momentum.


Introduction
Ultra-cold Fermi gas has been the focus of a lot of research investigations due to its highly controllable attractive interaction [1][2][3][4][5][6]. Comparing to other condensed matter systems, it is very difficult to establish the appearance of condensate or phase coherence in cold fermi gases, because of the lack of transport measurements. A decade ago, the superfluidity of Fermi gases has been proved experimentally by the observation of the vortex lattices [7]. However, this method still requires a fast sweep of the attractive interaction between the fermions to the deep Bose-Einstein condensation (BEC) limit, in order to see the density depletion in the vortex core. Recently, an alternative approach based on measuring the collective (Nambu-Goldstone) mode through the Bragg spectroscopy has been carried out experimentally to establish the phase coherence of unitary Fermi gases [8].
Due to the inhomogeneity of the clouds of Fermi gases, the integrated Bragg spectroscopy is difficult to extract useful information to indicate the existence of the condensate. This issue has been overcome by a method invented in [8], which makes the local observation of Bragg spectroscopy possible.
A lot of information of the many-body system can be inferred from the dynamic structure factor, such as the excitations related to the pair breaking and quasiparticle scattering. Importantly, the Nambu-Goldstone mode due to the spontaneous breaking of U(1) symmetry in the superfluid/superconducting phase can also be deduced from it, which can help to judge the onset of the ordered phase. A theoretical account of Bragg spectroscopy requires a consistent treatment of the density-density response function or the dynamical density structure factor (susceptibility) [9,10]. The Noether current associated with the continuous U(1) symmetry must be conserved no matter whether it is broken or not, which can be guaranteed in a gauge invariant theory. Hence one must take into account of the contributions from condensate, pair fluctuation and collective modes in a systematic way. A naive tree-level calculation of the density response function cannot even produce the correct locations of poles associated with the collective modes, where the tree-level means that calculation is the zero-th order of the electromagnetic coupling constant. Due to the strong correlations occurred in the unitary limit, the BCS-BEC crossover theory of unitary Fermi gases has to employ certain approximations to compute various thermodynamic properties [11,12]. It is a challenging problem to maintain the gauge invariance of the linear response theory when extensions of these approximations are included. The 'consistency' of the theory can be thought of as an important constraint since no exact solution is known for strongly correlated Fermi gas. This can guarantee that the approximation of the linear response theory is in the same way as how the self-energy is approximated for the strongly correlated Fermi gas.
Important progresses have been made in the theoretical framework of the structure factor. Early development of quantum Monte Carlo (QMC) approach was focused on low temperature simulation [13], and later applied to inhomogeneous unitary Fermi gas at finite temperature [14]. Diagrammatic technique was applied in [10,15,16], and the random-phase approach (RPA) was developed in [17]. The latter was later generalized to the RPA on top of the superfluid local density approximation (SLDA-RPA) approach.
In this paper, we adopt an G 0 G pair fluctuation theory to compute the density response function of unitary Fermi gases. We show that a gauge invariant linear response theory in which the approximation exactly matches that used in the thermodynamic calculation can be systematically constructed. In this way, the current conservation and the corresponding sum rules are satisfied exactly. Moreover, this theory has the advantage of easy implementation in the numerical calculations.
This paper is organized as follows. We first brief introduce G 0 G pair fluctuation theory both below and above T c in section 2. Then we show how to construct a gauge invariant linear response theory in section 3. In section 4, we present the numerical results of our theory and compare with experiments.

A brief introduction to the pair fluctuation theory
The construction of the gauge invariant linear response theory is based on the t-matrix approach with ladder diagrams made by one bare and one fully dressed Green's functions, which is known as the 'G 0 G' theory. This approach is inspired by the early work of Kadanoff and Martin [18]. A detailed review of this theory can be found in [4]. This asymmetric choice of the ladder series is more compatible with the BCS-leggett ground state [18,19]. Before explaining the main idea of the theory, we first fix some notations as follows.
Assuming that m is the particle mass and μ is the chemical potential, the general one-particle Green's function is ) , n is an integer, β=1/T and T K n k å = å å . It is a central task of various theories for interacting Fermi gases to find the most accurate expression of the self-energy Σ(K ).
For the strongly attractive Fermi gas, the pair fluctuations can be treated as virtual non-condensed pairs in equilibrium with the condensate of Cooper pairs. To describe the effects of non-condensed pairs, we introduce the t-matrix t pg (Q) which can be thought of as an amputated propagator obtained by truncating the infinite series of equations of motion of multi-particle Green's function [20,21]. Below T c in the superfluid phase, the self-energy can be decomposed into two parts. Aside from the usual BCS self-energy Here the pair propagator is given by the summation of infinite ladders made of bare and full Green's functions, and g is the contact coupling constant, which is related to the s-wave scattering length a s through This renormalization relation cancels precisely the ultraviolet divergence in χ(K ). The Bose-Einstein condensation condition, or the Thouless criterion, of the non-condensed pairs in this case can be expressed as the vanishing of the 'pair chemical potential', which is equivalent to the divergence of the pair propagator at zero momentum t g G K GK . It is easy to verify that it is simply the BCS gap equation. Therefore, the G 0 G theory reduces to BCS mean field theory at T=0. There is an undetermined pseudogap self-energy in the t-matrix which determines the pseudogap selfenergy itself. Therefore, the full G 0 G t-matrix theory requires a self-consistent solution to Σ pg from a set of coupled integral equations, which is too complicated in practical calculations. One can employ an approximation to simplify the final result. The pair condensation condition t 0 0 pg 1 = -( ) implies that the main contribution to Σ pg comes from the vicinity of Q=0. Thus, one can simplify the convolution to a multiplication In this way, Σ pg takes the same form as the BCS self-energy, which greatly simplifies the numerics and also provides an explicit expression for the pseudogap Δ pg . Hence the total self-energy is given by Then the energy gap can be decomposed into the superconducting (sc) and the pseudogap (pg) parts 2 Above T c , there is no condensate and Δ sc =0. We expect a non-zero pair chemical potential appearing in the following expansion Accordingly the gap and pseudogap equations are replaced by So far we have summarized the most important ingredients of G 0 G pair fluctuation theory. This theory has been applied to compute various thermodynamic quantities in the BCS-BEC crossover of ultra-cold Fermi gases, and agrees with the experiments well [22]. A gauge invariant linear response theory based on it must include the approximations in the same way as how the pair fluctuation effect is included in the self-energy.

Gauge invariant linear response theory with pair fluctuation effects
The Hamiltonian of the ultra-cold Fermi gas has a U(1) symmetry, of which the Noether current J μ satisfies the conservation law J 0 ¶ = m m . The density-density response function function can be obtained from a linear response theory subjected to an effective external electromagnetic (EM) potential A μ where μ=0, 1, 2, 3 is vector index of the pseudo Minkowski space with metric )is the four-momentum of external field. The EM response functions are given by )is the full interaction vertex to be determined later. The densitydensity response function is the '00' component of the tensor K mn . The density susceptibility and dynamical structure factors are Note when ω is large or T is low, there is almost no difference between χ″ and S since the factor coth T 2 w is almost 1.
In a gauge invariant linear response theory, the perturbed current must also be conserved, i.e.
, in the whole BCS-BEC crossover regime at general temperature. This can be guaranteed by the Ward-Takahashi identity (WI) of the EM response function q K Q 0 = m mn ( ) . These WIs satisfied by the response functions can be further inferred from the WIs satisfied by the interaction vertices and equation (13). The first WI can be easily verified, and the second WI indicates that the correction to the EM vertex must be in the same approximation as the self-energy effect is included in the Green's function. A gauge invariant linear response theory for BCS mean-field theory can be developed by the integral equation formalism [23], matrix linear response theory incorporating with consistent fluctuations of order parameter [24][25][26] or functional path integral approach [27]. When stronger-than-BCS attractive interaction is considered, the gauge invariance has to be maintained in a non-trivial way due to the pair fluctuation effects. There have been some formal discussions on this subject [27][28][29][30][31]. A general principle to find a gauge invariant interaction vertex is the same as proving WI in quantum field theory: inserting the bare EM vertex to the selfenergy diagram in all possible ways. For the pseudogap self-energy associated with the non-condensed fermion pairs, these insertions give rise to the Maki-Thompson(MT) diagram related to pseudogap, denoted by MT pg , and the two different Aslamazov-Larkin (AL 1,2 ) diagrams. For the BCS self-energy associated with condensed fermion pairs, these insertions give rise to MT diagram related to superfluid, denoted by MT sc , and also a contribution related to collective modes, denoted by Coll sc of which we haven't found a diagram representation [30]. In summary, the full vertex is now given by P Q P P Q P P Q P P Q P P Q P P Q P P Q P A diagrammatic representation of the full vertex is shown in figure 1. Here the MT sc diagram is obtained by inserting the EM vertex to the bare propagator of the superfluid self-energy Similarly, the MT pg diagram is obtained by inserting the EM vertex to the bare propagator of the pseudogap selfenergy The two AL diagrams can be obtained by inserting EM vertex to the pair propagator

P Q P t K t K Q G K P G K L G L Q L Q L G L P Q P t K t K Q G K P G K L G L Q
Note that due to the asymmetric choice of t-matrix, we have two different AL terms in G 0 G theory, while there is only on AL term in the Nozieres Schmitt-Rink(NSR) crossover theory [2]. The last one is the collective-mode diagram  = + and the response functions Q ij are given in appendix A. This diagram only appears in the superfluid phase (below T c ) where the symmetry is spontaneously broken. In the strict mean field theory, this diagram can be obtained by inserting the EM vertex into the gap equation at all possible places [31].
The full EM vertex given by equation (17) satisfies the Ward identity (16) q K Q K G K Q G K , 1 1 A proof by combining equations (B4), (B1) and (B6) can be found in appendix B. The full density-density response function is now expressed by

numerical results and discussions
Now we present the numerical results of our theory and compare them with the experiments. We first calculate the dynamic density structure factor of uniform unitary Fermi gases from low to high temperatures with relatively small incoming momentum. The results are shown in figure 2 where k k 0.5 F = in (a) and k 1.0 F in (b). Here k F is the Fermi momentum of a non-interacting Fermi gas with the same density. In the G 0 G theory, the critical temperature is T T 0.26 c F  at unitary limit, which is about 40% higher than the experiment value T c /T F =0.16 [32]. This difference is due to the inaccurate estimation of the mean field background, which is also known from other thermodynamic quantities. However, we can expect that the dynamical density susceptibility is generally insensitive to this background.
When T T c < , there is clearly a peak associated with the Nambu-Goldstone mode due to the spontaneous breaking of the U(1) symmetry in the superfluid phase, which gives rise to a collective motion of the condensate. When the transferring momentum is small enough, this corresponds to the well known Bogoliubov-Anderson mode. In the numerics, the imaginary part of ω, which comes from the complex continuation i i0 n w w  + + in the response functions, is assumed to be infinitesimally small below T c . Since the inverse of this imaginary part is recognized as the lifetime of the pairs, this assumption is reasonable because the Cooper pairs has a very long lifetime below T c . This makes the peaks in figure 2 much sharper than those of the experiments. There are also two single-particle continua, one corresponds to the quasiparticle scattering process with . For details please refer to the expressions of the response function list in appendix A. When the transfer frequency is low enough, the former branch of the continuum becomes quite apparent in the curves denoted by the T T 0.1 c = , T 0.5 c . This agrees with the recent experiment [33]. In the inset of panel (a), we present a full scale diagram of χ″ to give a comparison between the peak and the continuum. When T>T c , the collective-mode peak disappears due to the restoration of the symmetry, and the continua join together. As we increase the transferring momentum k, the collective-mode peak approaches the pair-breaking continuum (see figure 2(b)). Finally they will merge together if k is large enough.
In figure 3, we now plot the dynamical density susceptibility as a function of ω with large transferring momentum. As we discussed before, all continua and the collective-mode peak merge together when k is large enough. The five different curves correspond to T/T c =0.42, 0.72, 0.96, 1.26 and 1.60 respectively with k/k F is fixed to be 4.1. Below T c , χ″ displays a sharp peak roughly at 1/2 of the so-called recoil frequency k m 2 Since the transferring momentum is fixed, then all peaks emerge at the same position. The sharp peak of χ″ clearly corresponds to the scattering of the pairs from momentum conservation consideration. Since the lifetime of the pair is assumed to be long enough below T c as we discussed previously, the peaks in figure 3 are also much sharper than those of the experiments. Above T c , this sharp peak disappears again due to vanishing of condensed pairs, and only a very broad continuum roughly centered around ω r is left, which corresponds to the preformed pair breaking and quasiparticle scattering. Here we have assumed that the recoil frequency in the response functions has a nonzero imaginary part i 1 w + t above T c after implementing analytical continuation, since the non-condensed pairs have a finite lifetime. The lifetime is chosen approximately as to ensure a suitable width of the peak. This qualitatively agrees with what has been observed in the experiments [8].
In figure 4, we plot the dynamical density susceptibility χ″ as a function of T/T c for fixed large momentum k/k F =5 and fixed frequency ω=ω 0 , where ω 0 is the location of the peaks below T c and ω 0 =ω r /2 above T c .  in (b)). In the inset of panel (a), a full scale diagram of χ″(ω, k) at T T 0.9 c = is presented to show the comparison between the collective-mode peak and the single-particle continuum. Here n k 3 Below T c , the peak location is not exactly at ω r /2. Only when the momentum transfer is large enough, the peak location is close to ω r /2 since the characteristic length is much less than the size of condensed pairs. Hence we choose to plot the value of χ″ at the top of the peak, which makes the curve blow T c not very smooth. Above T c , there is no sharp peak any more, thus we simply plot the value of χ″ at ω=ω r /2. One can clearly see that the temperature dependence of χ″(ω 0 ) is roughly linear both below and above T c , but with quite different slopes. Above T c , χ″(ω 0 ) is almost flat since there are no condensed pairs, while below T c , χ″(ω 0 ) increases much faster with decreasing T. This indicates that the appearing of condensed pairs greatly increases the pair scattering. Therefore, 0 c w ( ) might qualitatively serve as the order parameter of the condensate because the difference between above and below T c mainly comes from the scattering of condensed pairs, which can be clearly found from the inset of panle (a) of figure 2 that the contribution from the peak is significantly larger even at T T 0.9 c = . In the insets, we show the blow up of the details around T c . One can see the flat curve above T c actually increases slightly with the increasing temperature. This trend is the same as the ideal gas.

Conclusion
In summary, we construct a manifestly gauge invariant linear response theory for a strongly correlated Fermi gas, from which the dynamic density structure factor can be studied near the Feshbach resonance. Our numerical results qualitatively agrees with the known experimental results.