Universal contact of strongly interacting fermions at finite temperatures

The recently discovered universal thermodynamic behaviour of dilute, strongly interacting Fermi gases also implies a universal structure in the many-body pair-correlation function at short distances, as quantified by the contact ${\cal I}$. This quantity is an excellent indicator of the presence of strong correlations in these systems, which provide a highly accessible physical model for other strongly correlated quantum fluids. Here we theoretically calculate the temperature dependence of this universal contact for a Fermi gas in free space and in a harmonic trap. At high temperatures above the Fermi degeneracy temperature, $T\gtrsim T_{F}$, we obtain a reliable non-perturbative quantum virial expansion up to third order. At low temperatures we compare different approximate strong coupling theories. These make different predictions, which need to be tested either by future experiments or advanced quantum Monte Carlo simulations. We conjecture that in the universal unitarity limit, the contact or correlation decreases monotonically with increasing temperature, unless the temperature is significantly lower than the critical temperature, $T\ll T_{c}\sim0.2T_{F}$. We also discuss briefly how to measure the universal contact either in homogeneous or harmonically trapped Fermi gases.


I. INTRODUCTION
Understanding strongly interacting fermions is one of the most challenging problems in present-day physics [1]. While the behavior of interacting fermions can be understood in some well-defined regions of parameter space, it remains elusive in the correlated strongly interacting regime [2,3]. A better understanding of strongly interacting fermions has wide-ranging implications for systems such as quark matter in neutron stars and hightemperature superconductors. An important generic idea in this field is fermionic universality [4,5]: all strongly interacting, dilute Fermi gases should behave identically, depending only on a scaling factor equal to the average particle separation, but not on the details of the interaction.
The recent realization of broad collisional (Feshbach) resonances with an interaction potential range r 0 small relative to the inter-particle spacing in ultracold atomic Fermi gases provides a highly controlled environment for studying the general problem of strongly interacting fermions [2,3]. By applying an external magnetic field across the resonance, the interparticle interaction can be accurately tuned from weak to infinitely strong [6]. This has led to the observation of crossover from Bardeen-Cooper-Schrieffer (BCS) superfluids to Bose-Einstein condensations (BEC) [7][8][9].
The most strongly interacting regime lies at the resonance, where the s-wave scattering length a s diverges and the two-body scattering amplitude reaches its maximum value allowed by quantum mechanics, i.e., it becomes unitarity limited [4]. A number of properties in the uni- * hhu@swin.edu.au † xiajiliu@swin.edu.au ‡ pdrummond@swin.edu.au tarity limit have been characterized. In particular the universal thermodynamic behavior [5,[10][11][12][13][14][15] has been experimentally observed, and agrees with approximate analytic theories and Monte Carlo simulations. However, it is still nontrivial to calculate these universal properties quantitatively, as there is no small interaction parameter in a perturbation theory expansion. In 2008, Tan gave new insight into this difficult problem by deriving a set of relations which link the shortdistance, large momentum correlations to the bulk thermodynamic properties of a fermion system with shortrange interparticle interactions [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31]. These relations generically apply to any dilute Fermi gas with an interparticle spacing much larger than r 0 . All of these Tan relations are exact in the limit of a vanishingly small range of the interaction potential, so that k F r 0 → 0, where k F is the Fermi momentum. They are connected by a single coefficient I, referred to as the integrated contact intensity or "contact". For instance, the pair correlation function is predicted to diverge as I/(16πr 2 ) at short distance r(≫ r 0 ) → 0 [16] and the momentum distribution will fall off as I/k 4 at large momentum k(≪ 1/r 0 ) [17]. The original definition of the contact given by Tan is, where ρ σ (k) is the momentum distribution in one spin component. The contact I is an extensive quantity and has the unit of N k F , where N is the total number of atoms of the system. For a homogeneous Fermi gas, it is also convenient to use a contact intensity, C = I/V , where V is the volume of the system. The immediate significance of the Tan relations can be seen most clearly from Tan's adiabatic sweep theorem [17], which states that the rate of adiabatic change of the total energy with respect to the inverse scattering length is proportional to the contact. While it is surprising that the short-range, high-energy correlations can be related so seamlessly to the low-energy equation of state, the Tan relations are closely related to the Feynman-Hellman theorem [32]. This method has also been used to calculate correlations in one-dimensional Bose gases at finite temperature from thermodynamic properties [33], with experimental verification [34].
In the fermionic strongly interacting regime, universal thermodynamics implies a universal contact. The determination of the contact therefore provides a very important means of characterizing the many-body properties and phases of strongly interacting fermions, complementary to existing measurements of the equation of state. As it is related to the derivative of the observed energy, the contact is a more sensitive test of theoretical predictions than direct thermodynamic measurements. We note that some of the Tan relations have been confirmed experimentally at JILA in the USA [35] and at Swinburne University of Technology in Australia [36]. It is compared with a zero-temperature perturbation prediction (solid line) [36] and a finite-temperature second-order virial expansion calculation at 0.5TF (dashed line) [25]. The measurements include Bragg spectroscopy (empty squares) [36], momentum distributions (solid or empty circles) [35], rf spectroscopy (stars) [35], and molecular fraction (triangles) [38].
Experimentally, the contact is obtained by a number of methods, including measuring: • the molecular fraction [37,38], • the momentum distribution tail [35], • the radio-frequency (rf) spectroscopy signal at large frequencies [35,39], or • the spin-antiparallel static structure factor at large momenta, using two-photon Bragg spectroscopy [25,36]. Figure 1 summarizes all the measured contact data to date for a trapped interacting Fermi gas (symbols). The experimental results are compared with our earlier predictions of a zero-temperature result from an approximate Gaussian pair fluctuation theory (solid line) [36] and from a finite-temperature calculation using the quantum virial expansion (dashed line) [25]. The experimental data of course are not taken exactly at zero temperature. Accurate determination of temperature is difficult due to the absence of an efficient thermometry. A typical experimental estimate gives T /T F ≃ 0.10 − 0.20. The data in Fig. 1 agree reasonably well with the zerotemperature curve, although near the unitarity limit they seem to lie systematically lower, possibly due to the nonzero temperature in these experiments. A similar theoretical calculation at zero temperature has been reported by Werner and collaborators [38] by interpolating the known perturbation results in the BEC and BCS limits.
In this paper, following the pioneering studies by Yu, Palestini and co-workers, we systematically investigate the temperature dependence of the contact for the most interesting case of a unitary Fermi gas. Yu et al. estimated the finite-temperature behavior of the universal contact [40], based on an assumption of phonon excitations at low temperatures and a second-order quantum virial expansion at high temperatures; while Palestini et al. presented ab-initio results using a non-self-consistent T -matrix approximation [41]. Both studies predicted a pronounced maximum in the contact at finite temperatures of about 0.6T F and 0.2T F , respectively. The present investigation, however, calls into question this prediction of a temperature-dependent peak in the universal contact. Our main results are summarized below, following a two-fold motivation.
Firstly, we give a high-temperature quantum virial expansion study of the contact up to third order. This provides us a general idea of how the virial expansion converges with decreasing temperature, and therefore enables us to estimate the temperature window for its applicability. We find that the virial expansion of the contact is quantitatively accurate down to the Fermi degenerate temperature T F and 0.5T F for a homogeneous and trapped Fermi gas at unitarity, respectively. We show that, in exact analogy with the virial expansion of thermodynamic potential [42], the contact at unitarity can also be expanded in terms of universal, temperature independent coefficients, c n = ∂∆b n /∂(λ/a s ), referred to as the contact virial coefficients. Here, ∆b n is the n-th virial coefficient and λ is the thermal wavelength. This is a direct consequence of fermionic universality. We predict that for a uniform unitary Fermi gas, the second and third contact coefficients are given by, respectively, c 2,∞ = 1/π ≃ 0.318 and c 3,∞ ≃ −0.141. The subscript "∞" stands for the unitarity limit with a s = ±∞.
Secondly, we are interested in the low-temperature be-havior of the universal contact. In the strongly interacting regime, as there is no controllable small interaction parameter, a comparative study using different strong coupling theories is necessary [43][44][45]. We find that the theoretical descriptions of the finite-temperature contact from different strong coupling theories show considerable discrepancies near the critical temperature T c (∼ 0.2T F ) for the onset of superfluidity. In particular, the enhancement (peak structure) of the contact near T c , found earlier in a non-self-consistent T -matrix calculation [41], is completely absent in other strong coupling theories. Therefore, we conjecture that the universal contact decreases monotonically with increasing temperature, unless the temperature is much smaller than the critical temperature. This discrepancy remains to be resolved by future experiments and advanced quantum Monte Carlo simulations [46][47][48], exploring the critical temperature regime close to resonance. Our paper is organized as follows: in Sec. II we present the quantum virial expansion for the contact and calculate the universal contact coefficient at unitarity, for both trapped and homogeneous Fermi gases. We then discuss how to determine the contact from strong-coupling theories in Sec. III, using different means via the tail of momentum distribution, Tan's adiabatic sweep relation and pressure relation. Theoretical predictions for the uniform universal contact from different strong-coupling theories are compared with each other and with accurate virial expansion results in the high-temperature regime. In the next section (Sec. IV), we describe how to calculate the trapped universal contact from the homogeneous data and report the results of different strong-coupling theories and virial expansion. In Sec. V, we discuss briefly the measurement of the temperature dependence of the universal contact in a homogeneous Fermi gas. A summary and outlook is given in Sec. VI.

II. QUANTUM VIRIAL EXPANSION OF THE CONTACT
The quantum virial expansion has proved to be an efficient method for studying the high-temperature properties of ultracold atomic Fermi gases [42,[49][50][51][52][53][54]. This method utilizes the fact that in the high temperature limit, the chemical potential µ diverges to −∞ and the fugacity z ≡ exp(µ/k B T ) ≪ 1 is a well-defined small expansion parameter. Thus, one may expand any physical quantities of interest in powers of fugacity, no matter how strong the interaction is.
For the grand thermodynamic potential of a twocomponent Fermi gas, the expansion is given by [42,52], where Q n = T r n [exp(−H/k B T )] is the partition function of a cluster containing n particles and b n is the n-th virial expansion coefficient. The trace T r n takes into account all the quantum states of n-particles with a proper symmetry. The virial coefficient b n can be calculated from the cluster partition functions [42,52], for instance, and In practice, it is often convenient to separate out a background, non-interacting thermodynamic potential, n z n + · · · ]. Here, we have used the superscript "1" to indicate the noninteracting systems. The thermodynamic potential of interacting Fermi gases can then be rewritten as, n . For a homogeneous Fermi gas, The essential assumption of the quantum virial expansion is that the above expansion of Ω − Ω (1) converges for temperatures down to T c , the critical temperature for the onset of superfluidity. We note that the fugacity in the expansion may be larger than unity at these low temperatures. However, convergence is still possible given small enough virial coefficients ∆b n . For a homogeneous unitary Fermi gas, the second virial coefficient ∆b 2,∞ = 1/ √ 2 was known 70 years ago [55]. The third virial coefficient, ∆b 3,∞ ≃ −0.35501298, was recently calculated by the present authors [42] and confirmed experimentally in an accurate thermodynamic measurement by Nascimbène and co-workers [13].
The virial expansion of the contact follows directly from an alternative representation of Tan's adiabatic sweep theorem in the grand-canonical ensemble, This is simply because the adiabatic sweep theorem implies the first law of thermodynamics, which can alternatively be written as Therefore, using Eq. (6) we immediately obtain a quantum virial expansion for the contact: where we have defined the dimensionless contact coefficient, c n ≡ ∂∆b n /∂(λ/a s ). For a homogeneous system, we sometime use the contact intensity, C = I/V .
In general, the contact coefficient should be a function of λ/a s and hence is temperature dependent. In the unitarity limit where λ/a s = 0, however, we anticipate a constant, universal contact coefficient, similar to the universal virial coefficient ∆b n,∞ [42,49]. This is a manifestation of fermionic universality, shared by all systems of strongly interacting fermions [4,5].

A. Universal relation between the homogeneous and trapped contact coefficients
In exact analogy with the virial coefficient [42], fermionic universality leads to a very simple relation between the trapped and homogeneous contact coefficients at unitarity. Let us consider the contact of a harmonically trapped Fermi gas with the trapping potential, In the thermodynamic limit of ω → 0, we may use the local density approximation and neglect the discrete energy levels. The whole Fermi system is treated as many cells with a local chemical potential µ(r) = µ − V T (r) and a local fugacity Due to the constant contact coefficients, the spatial integration in the total contact I T =´dr[C(r)] can be easily performed. Here, we use the subscript "T " to indicate explicitly the trapped system. We find that, where the trapped contact coefficient is given by a universal relation, and Q 1,T = 2(k B T ) 3 /( ω) 3 is the single-particle partition function in harmonic traps and in the local density approximation. In the following, using the known solution of two-and three-fermion problems, we calculate the universal second and third contact coefficients, in both homogeneous and trapped configurations.

B. Second universal contact coefficient
The second contact coefficient of a homogeneous interacting Fermi gas can be obtained from the well-known phase-shift expression for the second virial coefficient [49,55], Here, E i b is the energy of the i-th bound state of the two-body attractive interaction and δ 0 (k) is the s-wave scattering phase shift, given by k cot δ 0 ≃ −1/a s +r 0 k 2 /2. To obtain the derivative ∂∆b 2 /∂(λ/a s ), let us choose the BCS side with a s < 0, on which the two-body bound state is absent. As the range of the interaction potential r 0 is vanishingly small so that k B T ≪ 2 /(mr 2 0 ) and k ≪ 1/r 0 , we find that dδ 0 (k)/dk ≃ −a s /[1 + k 2 a 2 s ] and, Near to the resonance, it is readily seen that ∆b 2 ≃ 1/ √ 2 + λ/(πa s ), giving rise to a homogeneous contact coefficient, To calculate the trapped second contact coefficient, we consider the second virial sufficient in an isotropic harmonic trap, which is given by [42,52], Here, ǫ rel,n = (2ν n + 3/2) ω is the n-th relative energy of two fermions with unlike spins and ν n satisfies the secular equation, 2Γ(−ν n )/Γ(−ν n − 1/2) = d/a s , with Γ the gamma function and d = 2 /(mω) the characteristic length scale of the harmonic trap. The non-interacting relative energy is ǫ (1) rel,n = (2n + 3/2) ω (n = 0, 1, 2, ...) and in the unitarity limit, the solution of ν n is known analytically, ν n,∞ = n − 1/2. It is easy to show that, Thus, we find that in the unitarity limit, whereω ≡ ω/(k B T ) ≪ 1 is the reduced trapping frequency. The sum over n can be exactly performed, leading to, The leading term in the above equation is universal, satisfying the universal relation Eq. (12). The second term (∝ω 2 ) is non-universal and is caused by the length scale of the harmonic trap [42]. It represents the finite-size correction to the local density approximation that we have adopted above. This correction however is extremely small ifω → 0, as anticipated. At the Fermi degenerate temperature, T F = (3N ) 1/3 ω/k B , we find that ω 2 ∝ N −2/3 ∼ 10 −4 , for the typical number of atoms N ∼ 10 5 in the current experiments [2].

C. Third universal contact coefficient
The determination of the third contact coefficient is more cumbersome. In the homogeneous case, it is extremely difficult to use the conventional method employing the three-particle scattering matrix to obtain the third virial coefficient. Therefore, it is more practical to determine firstly the trapped contact coefficient c 3,∞,T , and then to use the universal relations at low trap frequency to obtain the homogeneous result, c 3,∞ = 3 √ 3c 3,∞,T . We note that this method relies on exact three-body solutions to the energy eigenvalues in a trap. These are known analytically for all eigenstates, and can be summed numerically to any required accuracy.
An estimate of c 3,∞,T can already be obtained by the known results of ∆b 3,T as a function of the coupling constant 1/k F a s at different temperatures T /T F and at a smallω ∼ 0.1 (see, for example, Fig. 3b in Ref. [42]), simply because, We find that the coefficient c 3,∞,T at resonance is indeed nearly temperature independent and estimate from the slope of ∆b 3,T that, An accurate determination of c 3,∞,T requires a systematic extrapolation to the limit ofω = 0. For this purpose, we calculate numerically the derivative c 3,∞,T (ω) = [∂∆b 3,T /∂(λ/a s )] λ/as=0 as a function ofω. Using the smallω data, a numerical extrapolation toω = 0 gives rise to the trapped third virial contact coefficient, Thus, we obtain immediately from the universal relation, Eq. (12), the homogeneous third contact coefficient, D. Large-T contact: the homogeneous case We are now ready to calculate the universal contact in the high temperature regime. For a homogeneous Fermi system, the single-particle partition function Q 1 = 2V /λ 3 and the dimensionless contact I/(N k F ) is given by, (23) Here N ≡ ρV is the total number of atoms with the homogeneous density ρ.
The fugacity z is determined by the number equation [52],ρ where we have defined a dimensionless densityρ ≡ ρλ 3 /2 = [4/(3 √ π)](T F /T ) 3/2 and the density of a noninteracting Fermi gas as In practice, for a given fugacity, we calculate the dimensionless density using Eq. (24) and hence the reduced temperature T /T F . The dimensionless contact is then obtained from Eq. (23), as a function of T /T F or the inverse fugacity z −1 .
the leading temperature dependence of the contact is given by, as predicted by Yu and co-workers [40]. We note however that the prefactor there is smaller by a factor of 4π 2 , due to a different definition for the contact. For a trapped Fermi gas at unitarity, the dimensionless contact can be written as, (28) The number equation takes the form [52], The trapped virial coefficient is given by ∆b n,∞,T = ∆b n,∞ /n 3/2 . In analogy with the homogeneous case, for a given fugacity we determine the reduced temperature T /T F from the number equation (29) and then calculate the trapped contact using Eq. (28).  Figure 3 presents the virial expansion prediction for the trapped universal contact, expanding up to the second order (dashed line) or third order (solid line). Amazingly, because of the factor of n −3/2 reduction for the n-th contact coefficient in harmonic traps, the convergence of the expansion is much improved. The expansion now seems to be quantitatively reliable down to 0.5T F . The asymptotic behavior of the contact at very high temperatures can be determined by setting We find that, in agreement with our previous result from the large momentum behavior of static structure factor [25]. Thus, the contact in harmonic traps decays at high temperatures much faster than in homogeneous space, due to the reduction of the peak density at the trap center at high temperatures.

III. LOW-T HOMOGENEOUS CONTACT FROM STRONG COUPLING THEORIES
At low temperatures, we have to resort to strong coupling theories to determine the contact. In the absence of a controllable small interaction parameter, however, there is no known a priori theoretical justification for which strong coupling theory is the most appropriate [45]. Here we consider three different theories within the many-body T -matrix approximation: the Nozières and Schmitt-Rink (NSR) Gaussian pair fluctuation theory [7,[56][57][58][59], non-self-consistent G 0 G 0 theory [60][61][62], and self-consistent GG theory [63][64][65].

A. Brief review of T -matrix approximations
The T -matrix approximation is also referred to as ladder approximation, which involves infinite set of Feynman diagrams -the ladder sum in the particle-particle channel. It is now generally accepted that this ladder sum is necessary for accounting for strong pair fluctuations in the strongly interacting regime and is the leading class of all sets of diagrams [66,67]. The need for including infinite sets of Feynman diagrams is self-evident, since we are dealing with a strongly interacting theory in which there is no small coupling constant that would allow a truncation of perturbation theory to finite order.
There are several T -matrix theories, differing in the diagrammatic structure of the T -matrix t(Q), the particleparticle propagator χ (Q), and the self-energy Σ (K). The ones we consider here are respectively given by the following form [66,67], Here and throughout, Q = (q, iν n ), K = (k, iω m ), while is the bare contact interaction renormalized in terms of the s-wave scattering length. We follow the conventional notations and use ǫ k = 2 k 2 /(2m) and K = k B T m k , where q and k are wave vectors, while ν n and ω m are bosonic and fermionic Matsubara frequencies. The Greek subscript α in the expressions for the propagator χ (Q) and the selfenergy Σ (K) may either be set to "0", indicating a noninteracting Green's function G 0 (K), or be absent, indi-cating a full interacting Green's function G (K). The latter can be calculated once the self-energy is determined, using the Dyson equation, By taking different values of α, two different T -matrix approximations are obtained: either self-consistent (no α), or non-selfconsistent (α = 0) . We note that, in the superfluid phase the Green's function has to be a 2 × 2 matrix, accounting for the U (1) symmetry breaking. The simplest choice in the T -matrix approximation is to take a non-interacting Green function G 0 (K) everywhere in χ (Q) and Σ (K). This idea was pioneered by Nozières and Schmitt-Rink (NSR) for an interacting Fermi gas in its normal state [7]. In addition, an approximate Dyson equation for the Green's function, i.e., was used. Interestingly, this non-self-consistent NSR approach combined with a truncated Dyson equation can be reformulated using the functional path-integral language for the grand thermodynamic potential [56], expanded to second order to include Gaussian fluctuations. From this perspective, the NSR theory is the first step in an orderby-order expansion in path-integral fluctuations, which has a systematic path to including higher-order terms. The NSR theory was extended recently to the brokensymmetry superfluid phase by several authors in different manners [57,59,[69][70][71], some of which involve assumptions in order to ease the computational workload. A full extension following the original NSR approach was reported by the present authors in 2006 [57], with the use of a mean-field (2 × 2 matrix) BCS Green's function as "G 0 ". Hereafter we shall refer to this extension as a Gaussian pair fluctuation theory or the GPF approach.
The truncation of the Dyson equation, Eq. (36), can be avoided, as shown by Strinati [60,61] and Combescot [62] and their co-workers. In particular, in the superfluid state the coupled T -matrix equations have been solved [61]. In the following, we shall refer to this non-selfconsistent T -matrix approximation with no truncation in the Dyson equation as the G 0 G 0 theory. However, since this approach is no longer part of a systematic expansion of the path integral, it is an open question as to whether it is more or less accurate than the NSR theory.
More sophisticated strong-coupling theories can also be obtained by using the full (dressed) Green function G (K) in χ (Q) and Σ (K). This modification, referred to as the self-consistent GG theory, was investigated in detail at the BEC-BCS crossover by Haussmann and collaborators [63,65], above and below the superfluid transition temperature. One advantage of the GG approximation is that the theory satisfies the so-called Φ−derivability and is thus conserving, although the conservation is at the single-particle level only. The GG theory has also been intensively studied in the condensed matter physics, particularly for high-temperature materials, under the name of the fluctuation exchange (FLEX) approximation [66,67]. Below threshold the GG approximation is difficult to treat at unitarity, and generally involves additional assumptions in order to obtain a scale-invariant theory.
It is clear that in both the NSR or G 0 G 0 approximations one omits infinite classes of diagrams that are responsible for multiparticle interactions. The fully selfconsistent GG theory attempts to correct for this, by modifying the single-particle Green function in the ladder diagram. However, the more crucial interaction vertices remain unchanged. At this stage, there are no general grounds for deciding which strong-coupling theory is the most appropriate. Therefore, any theoretical predictions produced by different strong coupling theories should be treated conservatively and be examined critically using more accurate ab-initio quantum Monte Carlo simulations or future experiments. We note that, for the equation of state of strongly interacting fermions, a comprehensive comparison between different strong-coupling theories and most recent experimental results was reported by the present authors [45]. In this earlier comparison we also included intermediate (GG 0 ) theories which are partially self-consistent [68]. These are omitted here for simplicity.
In the following, we calculate the universal contact at unitarity by using the GPF or NSR theory and compare our results with these reported by Strinati and co-workers (G 0 G 0 ) [41] and with the predictions provided by Haussmann and co-workers (GG) [65,72,73]. There are a number of ways to calculate the contact from T -matrix theories according to the different Tan relations, as follows: 1. Following Tan's original definition [16], the most direct way seems to be the use of the large wavevector asymptotic behavior of the fermionic momentum distribution, Eq. (1). The theoretical contact from the G 0 G 0 theory was obtained in this manner by Strinati and co-workers [41]. In the JILA experiment, the contact was extracted by averaging k 4 ρ σ (k) over an interval at k C,expt ∼ 2k F [35]. Figure 4 shows k 4 ρ σ (k) for a unitary Fermi gas in its normal state at T = 0.3T F , obtained by solving (iteratively) the coupled T -matrix equations [64], Eqs. (32), (33), (34), and (35). The power-law k −4 tail is clearly identified since k 4 ρ σ (k) approaches a constant above a characteristic wavevector k C ∼ 5k F , which is however larger than k C,expt . The small wiggles in Fig. 4 at about 10k F are due to the numerical inaccuracy. We note that the power-law tail of k −4 in the momentum distribution was first pointed out by Haussmann [63], using the self-consistent GG theory for a normal interacting Fermi gas. In the GG theory, the tail in ρ σ (k) implies that the contact can be directly calculated from the order parameter ∆ and the vertex function Γ(x, τ ), A brief explanation of this expression is given in the Appendix A. We also refer to Sec. IID of Ref.
2. Another way to calculate the contact is to use Tan's pressure relation [17], which results directly from the adiabatic sweep relation. As discussed in Ref. [17], this is evident if we write the total energy in the form, where f E is a dimensionless function. The adiabatic sweep relation leads to, 2 I/(4πm) = −N ǫ F f ′ E /k F , where the derivative is taken with respect to the variable 1/(k F a s ). We then find that, A direct application of the pressure relation in the unitarity limit is prohibited as a −1 s = 0 and therefore the dependence on the contact is lost. However, we can use the pressure relation to calculate the contact around the unitarity and then do a smooth interpolation to the limit of 1/(k F a s ) = 0.
3. Finally, the contact can also be calculated from the entropy by using, which can be obtained directly by taking temperature derivative in Eq. (7). This important relation was first shown by Yu and co-workers [40]. Thus, the temperature dependence of the contact is determined by the variation in the entropy with respect to the coupling constant.  Figure 5. (Color online) Temperature dependence of the contact in the unitarity limit, obtained by using the variation in the entropy Eq. (42), Tan's pressure relation Eq. (38), and the momentum distribution tail Eq. (1). Here we use the Gaussian pair fluctuation theory [57]. Figure 5 reports the universal contact calculated using the above mentioned methods, within the GPF (NSR) theory. There is an apparently discontinuous behaviour at the critical temperature T c , due to the breakdown of the NSR approach just above T c [44,45]. We find that all three theoretical methods yield nearly the same prediction, owing to the consistency in thermodynamic relations. The self-consistent GG theory has the same advantage of thermodynamic consistency.   Table I. Zero temperature contact in the unitarity limit. The contact I/(N kF ) = 6πζ/5 can be calculated from the ground state energy of an interacting Fermi gas near unitarity [17]: , where ξ and ζ are two universal parameters. The contact may also be determined from the tail of spin-antiparallel static structure factor [25], S ↑↓ (q) = I/(4N q). Here, we compare different results obtained from the three strong-coupling theories with the predictions by using ζ, either measured at ENS (ζ = 0.93(5)) [14] or extracted from quantum Monte simulation for ground state energy (ζ = 0.901(2)) [77] or for the static structure factor (or pair correlation function (ζ = 0.90)) [74].
zero temperature results are listed in the Table I, as well as the contact calculated by using the low temperature experimental data of the equation of state [14] and by using the quantum Monte Carlo data for pair correlation function [74] or ground state energy [77]. The ex-perimental data of course is not at zero temperature, as this is not experimentally achievable. In these measurements, low temperature thermometry was difficult. Similar experiments carried out elsewhere have reported lowest achievable temperatures of around T /T F ≃ 0.05. At low temperatures, we observe two distinct predictions for the behavior of the contact with increasing temperature: while the G 0 G 0 theory predicts an increase of the contact and hence a maximum around the critical temperature T c , both the GPF and GG theory suggest that the contact decreases monotonically as the temperature increases, unless T ≪ T c . This qualitative discrepancy deserves a careful analysis.
There are known arguments for the presence of a maximum in the contact at low temperature, as follows: • The enhancement of the contact just above T c in the G 0 G 0 theory was interpreted to be related to the strengthening of local pairing correlations in the absence of long-range order [41].
• As shown by Yu and co-workers [40], phonon excitations at low temperatures will enhance the contact as T 4 , causing a growth in the contact with temperature.
• The assumption of Fermi-liquid behavior for a weakly coupled Fermi gas [40] would also lead to a maximum in the contact at T ∼ T F .
These arguments, however, are not convincing enough to conclude there is a real maximum in the low-temperature contact in the unitarity limit, as we discuss below. The difficulty is that all the strong coupling theories are less accurate near criticality, suffering from strong pair fluctuations. Thus, the enhancement of the contact just above T c in the G 0 G 0 theory might be related to the breakdown of the approximations in this approach. Indeed, we observe a very similar enhancement of the contact in the NSR theory. However, we know that the NSR approach simply breaks down in the temperature region just above criticality.
On the other hand, while phonon excitations are important at low temperatures, our realistic calculation using the GPF theory at low temperatures suggests that the power-law of T 4 due to phonons is exhibited at T < 0.07T F ≪ T c only, as clearly shown in Fig. 7a. At such low temperatures, the contribution from phonons is of relative order 10 −3 . This is so nearly negligible that it might be easily be compensated by other excitations in this strongly interacting superfluid at unitarity.
We have also investigated the Fermi-liquid theory of the contact for a weak-coupling Fermi gas, as reported in Fig. 7b, where we compare the Fermi-liquid result with the NSR prediction for a weakly interacting Fermi gas with k F a = −0.5. The essential idea of the Fermi-liquid theory is that at low temperatures, the entropy density is given by s = m * k F k 2 B T /(3 2 ), where the effective mass m * = m[1+8(7 ln 2−1)/(15π 2 )k 2 F a 2 s ]. Therefore, by using Eq. (42), the derivative of the contact with temperature is given by, This leads to, where the zero temperature homogeneous contact I 0 can be calculated using the Lee and Yang's ground state energy (up to (k F a s ) 2 ) and has the form [17], It is clear from Fig. 7b that the weak-coupling Fermiliquid prediction works at very low temperatures only and strongly over-estimates the contact once T > 0.2T F , although the theory describes qualitatively the correct behavior of an increasing contact up to T ∼ T F . Near the unitarity limit, the use of Eq. (44) and hence its support for an enhancement of the contact is certainly doubtful.
In brief, we feel that in the unitarity limit it is more reasonable to expect a monotonically decreasing contact as the temperature increases. It is interesting to note that, for a weaker coupling constant, 1/(k F a s ) = −1, all the strong-coupling theories give the same qualitative behavior of the contact, as shown in the inset of Fig. 6. In particular, the contact is predicted to decrease with temperature in the superfluid phase. It is natural to assume that the same decrease would happen as well for a unitary Fermi superfluid.
At high temperatures, we find that the NSR and G 0 G 0 theories predict essentially the same result, as the difference between the two theories becomes rather small in the high-temperature regime. However, both predictions lie systematically below the quantum virial expansion result, indicating the inaccuracies of these theories. In contrast, at high temperatures the self-consistent GG result approaches the second-order virial expansion prediction. None of the strong-coupling theories can produce correctly the third-order expansion result, since the many-body T -matrix approach fails to account for the three-particle scattering process.

IV. LOW-T TRAPPED CONTACT FROM STRONG COUPLING THEORIES
Let us now turn to the contact of a unitary Fermi gas in a harmonic trap. Under the local density approximation this is given by where ρ(r) is the density distribution and k F (r) = [3π 2 ρ(r)] −1/3 .
In the unitarity limit, the local contact [C/(ρk F )](r) is a function of the density only, through the reduced temperature T /T F (r). The density may be expressed using ρ = 8/(3 √ π)λ −3 (T F /T ) 3/2 , so that the trapped contact is given by, Here, for simplicity, we assume an isotropic trap. In the local density approximation, it is convenient to use a local inverse fugacity ζ (r) = e −µ(r)/kB T ≡ ζ 0 exp[V T (r) /k B T ] with ζ 0 = e −µ/kB T . The integration over the radius r can then be converted to an integration about the inverse fugacity, by using r = (2k B T /mω 2 )ℓ(ζ) with ℓ(ζ) ≡ ln(ζ/ζ 0 ). In a harmonic trap, the total number of atoms is related to the Fermi temperature by N = (1/3)(k B T F ) 3 /( ω) 3 and the Fermi wave-vector is given by k F = 2mk B T F / 2 . Thus, we find that, where the inputs, the local dimensionless contactĨ ≡ C/(ρk F ) and the local reduced temperatureT ≡ T /T F , are functions of the inverse fugacity ζ, and can be determined using the homogeneous equation of state. The temperature T /T F in the above equation (48) is yet to be related to the inverse fugacity ζ 0 . For this purpose, we rewrite the number equation N =´drρ(r) into the form, For a given ζ 0 , we use Eqs. (49) and (48) to calculate the reduced temperature T /T F and the trapped contact [I/(N k F )] T , respectively, and in turn express the contact as a function of the reduced temperature. The detailed procedure of numerical calculations can be found in the Appendix B.
A. Trapped contact at zero temperature At zero temperature, the trapped contact at unitarity can be found analytically. In this case, the density profile is known exactly, where the peak density n 0 and the Thomas-Fermi R T F are respectively given by, where ξ ≃ 0.41(1) is the universal parameter of a unitary Fermi gas [14]. Because the temperature is zero, the local dimensionless contact [C/(N ρk)] hom is spatially independent, as it can only depend on a scale factor which has already been included in the definition of this quantity.
Here, for clarity we explicitly indicate the homogeneous contact by the suffix "hom". Thus, the trapped zero-temperature contact is given by, Using k F = 2mk B T F / 2 in harmonic traps with k B T F = (3N ) 1/3 ω, we find then, where the two universal parameters ξ and ζ are related to the ground state energy of an interacting Fermi gas near unitarity [17], E/(N ǫ F ) = (3/5)[ξ − ζ/(k F a s ) + · · · ].  Table II. Zero temperature trapped contact in the unitarity limit. We compare different results obtained from the three strong-coupling theories with the experimental measurements at Swinburne [36] and ENS [14].
with the experimental measurements at the Swinburne University of Technology, Melbourne [36] and at ENS, Paris [14]. We note that, the Swinburne data point is determined from spin-antiparallel static structure factor measured by Bragg spectroscopy. The lowest temperature in the Bragg measurement is about 0.10 ± 0.02T F , which could explain in part the slightly lower value obtained for the contact. For the ENS data point, we have used ξ = 0.41(1) and ζ = 0.93 (5), as determined from the experimental data for ground state energy. The temperature of this measurement was below the limits of the thermometry used, but we can estimate this as 0.05 ± 0.05T F from comparisons with related experimental measurements. The symbols indicate the measurements at ENS, Paris (star) [14] and at the Swinburne University of Technology, Melbourne (square) [36]. Figure 8 reports the temperature dependence of the contact of a trapped Fermi gas at unitarity, obtained from different strong-coupling theories. We compare the theoretical results with two experimental data points measured at low temperatures, while at high temperatures, the results are contrasted with the accurate quantum virial expansion prediction.
At low temperatures, T < 0.3T F , there is a sizable discrepancy between the predictions from different strongcoupling theories. Nevertheless, all the predictions are now in qualitative agreement, as the maximum in the homogeneous contact found previously in the G 0 G 0 theory is washed out by the trap averaging [41]. At high temperatures, the difference between different theories become much smaller. This is in line with our observation for the n-th contact coefficient, which receives a factor of n 3/2 reduction in harmonic traps.
The trapped contact of a unitary Fermi gas can be measured readily in experiment, by extending the existing low-temperature measurements of either the momentum distribution [35] or Bragg spectroscopy of a unitary Fermi gas [36] to the high-temperature regime. The precise determination of the temperature in the unitarity limit is a challenging problem [12]. However, this can be overcome by using the known equation of state in the unitarity limit [13,43,45]. We anticipate that the quantitative discrepancy between different strong-coupling theories shown in Fig. 8 will be clarified by future experiments.

V. EXPERIMENTAL MEASUREMENT OF THE HOMOGENEOUS CONTACT AT UNITARITY
It would be very useful to determine the contact of a homogeneous unitary Fermi gas, without the complication caused by the trap averaging. The result would provide a unique opportunity to test quantum many-body theories of strongly interacting Fermi gases that we have discussed above. However, an in-situ measurement is generally required, which is difficult to setup experimentally. The determination of the universal homogeneous contact could be achieved as follows, by: • Measurement of homogeneous equation of state near unitarity. The Salomon group at ENS, Paris has already demonstrated the accurate measurements of uniform equation of state either in the unitarity limit [13] or at zero temperature [14]. The measurement could be extended straightforwardly to the general case with arbitrary coupling constant and temperature, say, for example, the entropy S(T /T F , 1/k F a). The homogeneous contact can then be determined directly using Eq. (42).
• Measurement of the tail of the spatially-resolved RF spectrum at unitarity. Tomographic RF spectroscopy was recently developed at MIT to measure the local pairing gap of a unitary Fermi gas of 6 Li atoms at finite temperatures [78]. Ideally, the local universal contact could be retrieved from the tail of RF-spectrum that obeys a ω −3/2 power-law. However, the analysis of the tail structure is complicated for 6 Li atoms due to a strong final state effect. It would be desirable to develop a spatiallyresolved RF-spectroscopy for 40 K atoms, for which the final state effect may be negligible.
• Measurement of the spatially-resolved Bragg scattering spectrum at unitarity [36,79]. Thus, one obtains the local spin-antiparallel static structure factor at large momentum q, which at unitarity satisfies the asymptotic behavior of S ↑↓ (q) = [I/(N k F )]/(4q/k F ). The local universal contact is then determined. All of these measurements are based on the local density approximation. The local contact is then to be determined as a function of the local chemical potential µ(r) or local inverse fugacity ζ(r) = exp(−µ(r)/k B T ). In Fig.  9, we show the inverse fugacity dependence of the universal contact of a unitary Fermi gas, to be confronted with future experiments.

VI. CONCLUSIONS
To summarize, we have presented a comparative theoretical study of the universal contact for a strongly interacting Fermi gas at unitarity. The contact at high temperature has been accurately determined by using a quantum virial expansion method, while at low temperatures, we have employed different strong-coupling theories to try to estimate the temperature dependence of the contact.
The temperature dependence of the contact near the critical temperature, however, remains unresolved, since one of the strong-coupling theories, the G 0 G 0 theory, predicts an enhancement or a maximum near criticality, which is not observed in the other two strong-coupling theories. We believe that this enhancement could be due to the inaccuracy of the G 0 G 0 approach and conjecture that the universal contact should decrease monotonically with increasing temperature, except possibly at very low temperatures T ≪ T c .
As there is no solid justification for any of these strongcoupling theories, our conjecture should be examined critically by future experiments or advanced quantum Monte Carlo simulations. Therefore, we have proposed several experimental ways to determine the homogeneous contact of a unitary Fermi gas, using measurements of equation of state, tomographic RF-spectroscopy, and spatially-resolved Bragg spectroscopy. All of these measurements are within reach of current experimental techniques.
We note finally that fermionic universality and the detailed behavior of a strongly interacting Fermi gas near the normal-superfluid transition is of great interest. Tan's contact provides an entirely new means to characterize the universal properties and phases of strongly interacting fermions, in addition to the equation of state. Our comprehensive study of the universal contact, built on the most recent theoretical methods, should provide useful insights for future research.
Note added.
-Since the submission of our manuscript, the finite-temperature contact of a trapped unitary Fermi gas was measured at Swinburne [80]. The experimental result agrees fairly well the theoretical predictions in Fig. 8. A quantum Monte Carlo simulation of the finite-temperature contact of a homogeneous unitary Fermi gas was also carried out by Drut, Lähde and Ten [81].
To do the integration, for instance, for the integral in the number equation (B5), we use the following discretized version, where the two types of integrals I (ζ i ) and I ∞ are defined by,