Universal thermodynamics of a strongly interacting Fermi gas: theory versus experiment

Strongly interacting, dilute Fermi gases exhibit a scale-invariant, universal thermodynamic behaviour. This is notoriously difficult to understand theoretically because of the absence of a small interaction parameter. Here we present a systematic comparison of theoretical predictions from different quantum many-body theories with recent experimental data of Nascimbene et. al. (arXiv:0911.0747v1). Our comparisons have no adjustable parameters, either theoretically or experimentally. A simple Gaussian pair fluctuation theory is shown to give the best quantitative agreement in the superfluid state below threshold. In the normal state, we also calculate the equation of state by using a quantum cluster expansion theory and explore in detail its applicability to low temperatures. Using the accurate experimental result for the thermodynamic function $S(T)$, we determine the temperature $T$ of a trapped Fermi gas at unitarity as a function of a non-interacting temperature $T_{i}$ which can be obtained by an adiabatic sweep to the free gas limit. By analyzing the recent experimental data, we find a normal-superfluid transition temperature $(T/T_{F})_{c}=0.19\pm0.02$ or $(T_{i}/T_{F})_{c}=0.16\pm0.02$ in a harmonic trap, where $T_{F}$ is the Fermi temperature for a trapped ideal, non-interacting Fermi gas.


Introduction
The recent discovery of broad Feshbach resonances in two-component atomic Fermi gases has opened a new era in the study of strongly interacting fermions [1,2,3,4,5,6,7]. By tuning an external magnetic field across the Feshbach resonance, the interatomic attractions can be changed precisely from weak to infinitely strong, leading to the observation of a crossover from a Bardeen-Cooper-Schrieffer (BCS) superfluid to a Bose-Einstein condensation (BEC). At resonance, the s-wave scattering length a s diverges (a s = ±∞) and the two-body scattering amplitude reaches the maximum value allowed by quantum mechanics due to unitarity. Many unique properties are anticipated in this strongly interacting limit, including a high superfluid transition temperature and an exotic normal state with a pseudogap.
Most interesting of all is fermionic universality. This means that all strongly interacting, dilute Fermi gases behave identically, regardless of the details of the interaction. Their properties depend only on temperature, together with a scaling factor equal to the average particle separation [8,9,10,11]. This limit promises to bring a new rigour and simplicity to the understanding of strongly correlated Fermi gases. Because of universality, it should be feasible to understand other strongly interacting Fermi superfluids from experiments in the highly controlled environment of an atomic physics laboratory. Possible examples include neutron stars and high-T c superconductors.
Intense experimental investigations have been carried out to understand fermionic universality, in particular, its implication to the thermodynamic properties [12,13,14,15,16,17,18,19,20]. Pioneering observations were carried out at Duke University [12,13,14,15], which made the first attempt to reach the unitarity limit in 6 Li gas in 2002 [12]. The stability of atomic Fermi gases with strongly attractive interactions was observed. The ground state energy was found to be reduced significantly compared to its ideal, non-interacting limit. The reduction factor β (or ξ = 1 + β) has now been determined accurately to within a few percent: β ≃ −0.59 ± 0.01, after substantial experimental effort [20].
In this paper we show that recent highly accurate measurements [19] on strongly interacting 6 Li give the most stringent test to date of fermionic strongly coupling many-body theories. In fact, these experiments determine the whole set of universal thermodynamic functions for a trapped Fermi gas at unitarity. As well as measuring bulk thermodynamic properties, the data can be used to determine the energy and entropy, E(T ) and S(T ), of a trapped gas. Due to the use of larger samples, the accuracy is even better than that achieved at Duke. With unprecedented precision, these new universal functions therefore provide an unbiased test of theoretical predictions. This comparison, without any fitting parameters, indicates that while a BCS-type mean field theory is certainly incorrect, an extension using a simple Gaussian pair fluctuation theory provides overall the best agreement with experiment, except at the critical superfluid transition region.
To give some background to these developments, the first energy and heat capacity measurements as a function of temperature were performed by Kinast et al. [13]. However, due to the lack of reliable thermometry in the strongly interacting regime, an empirical thermometry was used. Conversion of measured results to real temperature required a particular strong-coupling theory, and was therefore model-dependent. This difficulty of model dependence was overcome at the end of year 2006 [14,15], by means of direct measurements of entropy instead of temperature. In this way, both energy E and entropy S were determined without invoking any specific theoretical model. These pioneering and very important model-independent measurements had an accuracy at the level of only a few percent. At the same period, the potential energy of a strongly interacting 40 K gas was also measured at JILA [17]. The temperature was characterized in terms of the noninteracting temperature of an adiabatically equivalent ideal Fermi gas. This is therefore equivalent to an entropy measurement. These sets of experimental data, together with results of another 6 Li experiment at Rice, were analyzed by the present authors [11]. The result was that all the thermodynamic data lay on a single universal curve. This gave the first, very strong evidence for the universal thermodynamics of a strongly interacting Fermi gas [11].
In parallel with these ground-breaking experiments, there have been numerous theoretical studies of the thermodynamics of a strongly interacting Fermi gas. In the absence of exact solutions, the methods used were either strong-coupling perturbation theories [21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40] or ab-initio quantum Monte Carlo (QMC) methods [41,42,43,44,45,46,47,48,49,50]. However, a deeper understanding is made difficult by the absence of a controllable small interaction parameter [51]. The use of standard perturbation theories thus requires infinite order expansions. These typically require truncations of sets of diagrams which cannot be fully justified a priori. Numerically exact QMC calculations are very helpful and can provide unbiased benchmarks, provided that there is an appropriate extrapolation of the lattice results to large lattice size or, in the diagrammatic Monte Carlo case, to zero range potentials.
The first theoretical explanation of the heat capacity at unitarity was given by Chen et al. [13], using a pseudogap theory. In this study an empirical temperature was converted to an approximate real temperature. The present authors subsequently gave a theoretical prediction for the both the homogeneous and trapped equation of state at unitarity [30]. This used a Gaussian pair fluctuation (GPF) theory below threshold, thus extending an approach proposed initially by Nozires and Schmitt-Rink (NSR) [21,29] for the above threshold case. We showed that the conversion of empirical temperature to actual temperature is strongly model-dependent.
Thus, in principle one cannot obtain accurate information about the real temperature from these empirical temperature measurements, without a reliable strongcoupling theory. The model independent measurement of energy as a function of entropy, E(S), by Luo et al. [14,15] was therefore a crucial experimental advance. This provided the first data that could be used to quantitatively compare different strong-coupling theories without any free parameters. One such comparison was performed by the present authors [38], by using different perturbation theories and available QMC results [46]. Even so, it was still impossible to determine the dependence of the energy on temperature E(T ) and of the entropy on temperature S(T ), due to difficulties in determining the absolute temperature T . Moreover, the measurements and comparisons were restricted to the case of a trapped Fermi gas.
Most recently, a general method was developed by Nascimbne et al. at ENS to measure the bulk equation of state of a homogeneous Fermi gas of lithium-6 atoms [19], following a theoretical proposal by Ho and Zhou [52]. The local pressure P (µ(z), T ) or the local thermodynamic potential Ω(µ(z), T ) = −P (µ(z), T )δV (δV is the volume of a cell at the position z) of the trapped gas was directly probed using in situ images of the doubly-integrated density profiles along the long z-axis. The temperature was then determined by using a new thermometry approach employing a 7 Li impurity. The chemical potential could also be determined using the local density approximation, with µ(z) = µ 0 − V trap (z) and the central chemical potential µ 0 being determined appropriately. By introducing a universal function [53] h experimentalists were able to determine h(ζ) with very low noise. Here, ζ ≡ exp (−µ/k B T ), Ω(µ, T ) is the interacting thermodynamic potential and Ω (1) (µ, T ) is the thermodynamic potential of an ideal two-component Fermi gas. This precise measurement allows a direct comparison with many-body theories developed for a uniform Fermi gas [19]. Our main results may be summarized as follows. First, though the theoretical predictions from all the model approximations seem to fluctuate around the experimental data, it turns out that the simplest Gaussian pair fluctuation theory gives the best description of the observed thermodynamic properties, except in the vinicity of the superfluid transition point. Second, using the measured universal functions as a benchmark, we examine the applicability of a quantum cluster (virial) expansion method [54,55]. We find that for a trapped gas, up to the leading interaction effect (second order), the expansion is quantitatively reliable down to T ≃ 0.7T F . This limit can be decreased further to T ≃ 0.4T F , with inclusion of higher order virial coefficients (i.e., up to fourth order). Thus, we demonstrate clearly the usefulness of quantum cluster expansion in the study of a normal, but strongly interacting quantum gas.
Finally, we note that the temperature of a trapped Fermi gas at unitarity is often characterized by a non-interacting temperature T i obtained by an adiabatic sweep to the ideal gas limit. Using the accurate experimental data for S(T ), we calculate the relation T (T i ), which shows an apparent kink at low temperatures. We therefore determine a characteristic temperature (T /T F ) 0 = 0.19±0.02 or (T i /T F ) 0 = 0.16±0.02 for a trapped Fermi gas at unitarity, below which the thermodynamic functions start to deviate from normal Fermi-liquid-like behavior due to pairing effects. This paper is organized as follows: in Sec. 2 we briefly review different strongcoupling perturbation theories and the high temperature quantum virial expansion theory. A comparison for the bulk universal function h(ζ) of a homogeneous Fermi gas at unitarity is presented in Sec. 3. The validity of the quantum virial expansion in the uniform case is discussed. In the Sec. 4, we explain how to reconstruct the trapped universal thermodynamic functions from h(ζ), and present a systematic comparison of different strong-coupling theories with the accurate experimental data, for various thermodynamic functions. We also examine the applicability of the quantum virial expansion to a trapped Fermi gas at unitarity. In Sec. 5, we calculate the actual temperature as a function of the non-interacting temperature at the same entropy. A summary and outlook is given in Sec. 6.

Theoretical review
In this section we review two types of strongly interacting Fermi theories: the strongcoupling perturbation theories and a controllable quantum cluster expansion theory. We consider a two-component Fermi gas with equal spin populations. At ultracold temperatures (< 100nK), the interatomic interactions between atoms with unlike spins can be well-described by an s-wave scattering length a s . For the case of molecule formation with a broad Feshbach resonance, a simple two-species model that neglects the molecular field is very accurate, otherwise a full three-species model is necessary [56,57]. The hamiltonian of the system can then be written as, where ǫ k =h 2 k 2 /(2m) is the fermionic kinetic energy at wave number k, and is the bare contact interaction renormalized in terms of the s-wave scattering length a s .

Strong-coupling perturbation theories
We start with a brief overview of the most commonly used strong coupling theories for a Fermi gas at unitarity. These are approximate many-body T -matrix theories [58,59], since no exact results are known in three dimensions. Such theories typically go beyond BCS theory by including an infinite set of higher order Feynman diagrams. The diagrams included, known as the ladder sum in the particle-particle channel, are still not the complete set of all possible terms in perturbation theory. However, it is generally accepted that a ladder sum is necessary in order to include the strong pair fluctuations in the strongly interacting regime. This is expected to be the leading class of these infinite sets of diagrams [58,59]. However, there are differences in the procedures used to obtain the relevant diagrams that are included. to the thermodynamic potential. The solid line represents the single-particle Green's function, while the dashed line shows the bare contact interaction. Note that, these diagrams are valid at all temperatures, both below and above the normal-superfluid transition. However, below the transition, the single-particle Green's function should be a 2 × 2 matrix.
In more detail, we show in Fig. 1a the diagrammatic structure of the T -matrix [58,59], t(Q), where the sucessive two-particle scattering between fermions with unlike spins is taken into account to infinite order. This forms a ladder structure, with the solid line and dashed line representing, respectively, the single-particle Green function G and the interaction U. Consequently, as an effective interaction the T -matrix can be diagrammatically represented by, by summing all the successive scattering process. In the normal state with contact interactions, the ladder sum can be conveniently calculated as, Here and throughout, Q = (q, iν n ), K = (k, iω m ), and q and k are wave vectors, while ν n = 2nπk B T and ω m = (2n + 1)πk B T (n = 0, ±1, ±2, · · ·) are bosonic and fermionic Matsubara frequencies, respectively.
Different T -matrix theories differ in their choice of the particle-particle propagator χ (Q), and the associated self-energy, where we have introduced an energy-momentum summation, K = k B T ωm k . The subscripts α, β, and γ in the above equations may either be set to "0", indicating a non-interacting Green's function Universal thermodynamics of a strongly interacting Fermi gas: theory versus experiment7 or be absent, indicating a fully dressed interacting Green's function. In the latter case the Dyson equation, is required to self-consistently determine G and Σ. The only free parameter, the chemical potential µ, is fixed by the number equation, N = 2 lim τ →0 + K G(K)e iωmτ [60]. By taking different combinations of α, β and γ, there are six distinct choices of the Tmatrix, for which a notation of (G α G β )G γ will be used. As noted earlier, there is no known a priori theoretical justification for which T -matrix approximation is the most appropriate.
It is important to note that, while having the same diagrammatic structure, the T -matrix above and below the superfluid transition temperature T c are different, due to the use of different Green's functions G 0 or G. In the superfluid phase below T c , the Green's function has to be a 2 × 2 matrix, accounting for U(1) symmetry breaking. Accordingly, an additional parameter, the order parameter, appears.
The simplest choice, (G 0 G 0 )G 0 , was pioneered by Nozires and Schmitt-Rink for a normal interacting Fermi gas [21], with a truncated Dyson equation for the Green's function, i.e., This was shown to be equivalent to including the Gaussian pair fluctuations in the grand thermodynamic potential [22,24], which is shown diagrammatically in Fig. 1b. The NSR theory was extended recently to the broken-symmetry superfluid phase by several authors [25,26,29,61,62]. However, some of these approaches involved additional assumptions to reduce computational difficulties. A full extension of the original idea of Nozires and Schmitt-Rink to the below threshold regime was reported by the present authors [29], with the use of a mean-field (2×2 matrix) BCS Green's function as "G 0 " in the thermodynamic potential in Fig. 1b. In the following we shall refer to this extension as a Gaussian pair fluctuation or GPF approach. Numerical calculations were then performed at the BEC-BCS crossover for the equation of state of a homogeneous Fermi gas. Compared to the zero-temperature QMC simulation for ground state energy [42], we found that the extended GPF approach works extremely well in the superfluid phase. It provides a quantitatively reliable description of the low-temperature thermodynamics of a strongly interacting fermionic superfluid.
In greater detail, in our theory below the critical temperature, the contribution of T -matrix pair fluctuations to the thermodynamical potential takes the form (see, for example, Fig. 1b), where are respectively the diagonal and off-diagonal parts of the pair propagator. Here, G 11 and G 12 , used in Fig. 1b for the single-particle line, are BCS Green's functions with a variational order parameter ∆. Together with the mean-field contribution where the excitation energy , we obtain the full thermodynamic potential Ω = Ω 0 + δΩ. All the thermodynamic functions, including the total energy E and total entropy S, can then be calculated straightforwardly following thermodynamic relations. For consistency, in our formalism we determine the order parameter using the gap equation ∂Ω 0 /∂∆ = 0. Together with the number equation, n = −∂Ω/∂µ, we solve iteratively the two parameters µ and ∆. This theory with bare BCS Green functions in the pair propagators constitutes the simplest universal description of strongly interacting fermions, including the essential contribution from the low-lying collective Bogoliubov-Anderson modes. As such, this type of theory may have useful applications to other types of strongly interacting fermionic superfluids.
The GPF and NSR approximation does not attempt to be self-consistent. More sophisticated strong-coupling theories can be obtained by using dressed Green functions as pair propagators. For example, one may consider a (GG)G approximation, with a fully self-consistent propagator. This was investigated in detail by Haussmann et al. [23,32], both above and below the superfluid transition temperature. One advantage of the self-consistent (GG)G approximation is that the theory satisfies the so-called Φ−derivable approach to the many-body problem due to Luttinger and Ward, in which the exact one-particle Green functions play the role of an infinite set of variational parameters. The (GG)G theory is thus conserving.
An intermediate scheme with an asymmetric form for the particle-particle propagator, i.e., (GG 0 )G 0 , has been discussed in a series of papers by Levin and coworkers [27], based on the assumption that the treatment of fluctuations should be consistent with the simpler BCS theory at low temperatures. Although the (GG 0 )G 0 theory has been explored numerically to some extent [63], a complete numerical solution is difficult. A simplified version was introduced based on a decomposition of the Tmatrix t(Q) in terms of a condensate part and a pseudogap part, leading to the so-called "pseudogap" crossover theory [27]. In the present comparative study, we will include both the pseudogap modification as well as the full (GG 0 )G 0 theory. Due to numerical difficulties, we shall consider the full (GG 0 )G 0 theory in the normal phase only.
It is clear that in the GPF approximation one omits infinite diagrams that are responsible for the multiparticle interactions. The fully self-consistent (GG)G theory and partially self-consistent (GG 0 )G 0 theory attempt to correct for this, by modifying one or more single-particle Green's function in the diagrams. However, the more crucial interaction vertices remain unchanged. For brevity, hereafter we shall refer to the fully self-consistent (GG)G theory and partially self-consistent (GG 0 )G 0 theory as GG and GG 0 theory, respectively.
To close the subsection, we emphasize again that there are currently no general grounds to decide which strong-coupling theory is the most appropriate, due to the absence of a small controllable interaction parameter. However, as we shall see, these approaches do give distinct predictions which can be tested experimentally.

Scale-invariance and universal relation at unitarity
In the unitarity limit, due to the infinitely large scattering length, the interatomic distance becomes the only relevant length scale in the problem. The internal energy and entropy of the system therefore scale like, where T F = ǫ F /k B is the Fermi temperature, and f E and f S are two dimensionless universal functions. The scaling form leads to a well-known scaling identity in free space: which holds as well for any ideal, non-interacting quantum gases. To show this, we note that at unitarity the pressure of the gas can be readily determined from P = − [∂E/∂V ] N,S . From the expression of entropy (15), it is easy to see that holding the entropy invariant is equivalent to fixing the reduced temperature T /T F . Hence, the only dependence of the energy on volume is through the Fermi energy, i.e., E ∝ V −2/3 . Taking the derivative with respect to the volume, one finds that P = 2E/(3V ) or Ω = −2E/3. This simple equation relates the pressure or thermodynamic potential and energy for a strongly interacting Fermi gas at unitarity in the same way as for its ideal, noninteracting counterpart, although the energy would be quite different. The scaling identity Eq. (16) follows naturally from the thermodynamic relations. At this point, only the GPF (NSR) approach and the fully self-consistent GG theory satisfy it, since in both theories we can write down a well-defined thermodynamic potential and then derive from it other thermodynamic quantities in a consistent way. Other strong-coupling theory, more or less, runs into the thermodynamic inconsistencies, and therefore violates Eq. (16). We note that in the superfluid phase, an ad-hoc renormalization of the interaction strength is required in the fully self-consistent GG theory, in order to obtain a gapless phonon spectrum [32]. Hence, the GG theory in the supefluid phase does not appear to be universal without modification. In the calculations with the pseudogap theory and GG 0 theory, we shall obtain the energy and entropy from the chemical potential by integrating the thermodynamic relation n = −∂Ω/∂µ, in order to satisfy Eq. (16). The detailed procedure was outlined in Ref. [38].
In the experimental situation with a harmonic trapping potential V (x), potential energy enters into the total energy. Thus, the scaling identity should be modified accordingly: This is because the potential energy N V equals to the internal energy (3/2) dxP (x) in harmonic traps and then the internal energy is a half of the total energy. The prefactor of 2/3 in Eq. (16) is therefore reduced by a factor of 2.
To prove the equality, we may treat the gas as a collection of many locally equilibrium uniform cells (i.e., using the local density approximation). The local force balance arising from the pressure P (x) and the trapping potential V (x) gives rise to, Taking a inner product of the above equation with x and integrating over the whole space, we readily obtain for a harmonic trap and the integration by part for x · ∇P (x). Hence, the strongly interacting Fermi gas at unitary obeys the same virial theorem as for an ideal quantum gas.

High temperature quantum cluster expansion
At high temperatures there is a controllable, small parameter, given by the fugacity z = exp(µ/k B T ) ≪ 1. This is small because the chemical potential µ diverges to −∞ at large temperatures T . In principle, all thermodynamic properties of a interacting Fermi gas can be cluster expanded in powers of fugacity [54,55], even in the strongly interacting limit. The thermodynamic potential Ω (1) of an ideal, non-interacting uniform Fermi gas takes the form, where λ ≡ [2πh 2 /(mk B T )] 1/2 is the thermal wavelength and b (1) n = (−1) n+1 /n 5/2 is the n-th virial coefficient. We use the superscript "1" to indicate the non-interacting systems. While at the first glance Eq. (19) may be valid at z < 1 only, its applicability is actually much wider. The expansion is meaningful for arbitrary positive values of fugacity through an analytic continuation across the point z = 1. It is then possible to show that, as expected for a non-interacting Fermi gas, In the presence of interactions, the virial coefficients are modified. The thermodynamic potential can instead be written as, where ∆b n = b n − b (1) n . Our key assumption in using the quantum cluster expansion is that the expansion of Ω − Ω (1) might be applicable near to the critical temperature, despite the fact that the fugacity may already be much larger than unity close to the superfluid transition. It is possible to test this conjecture in either BCS or BEC limit, by analytically calculating the virial coefficient ∆b n and hence the radius of convergence of the expansion. We leave this possibility in a future study.
At unitarity, the virial coefficients are temperature independent and are known up to the fourth order: ∆b 2 = 1/ √ 2, ∆b 3 = −0.35501298, and ∆b 4 ≃ 0.096 ± 0.015. The second virial coefficient was already known 70 years ago [54]. The third coefficient was calculated recently by the present authors [55]. This theoretical prediction was confirmed experimentally in the accurate thermodynamic measurements of Nascimbne et al. [19]. These recent experiments were also able to determine the fourth coefficient empirically [19]. With these coefficients, using ζ = z −1 the universal function h(ζ) at high temperatures may be written as, The above discussion for a uniform gas can be easily extended to the case with a harmonic trap: V trap (r) = mω 2 r 2 /2 or, in the general case of a trap with cylindrical symmetry, with ω ≡ (ω 2 ⊥ ω z ) 1/3 . Within the local density approximation, the thermodynamic potential becomes position-dependent through a local chemical potential µ(r) = µ 0 − V trap (r) or a local fugacity z(r) = e βµ(r) . Using the fact that the virial coefficients are constant at unitarity, the total (integrated) thermodynamic potential Ω trap (µ 0 , T ) = drΩ(r) takes the from, where z 0 = e βµ 0 is the fugacity at the trap center. It is easy to see that the nth virial coefficient in a trap, b n,trap = b n /n 3/2 , is much reduced with respect to its uniform counterpart.
The thermodynamic potential of an ideal trapped Fermi gas is now given by, In analogy with the uniform case, we may define a universal function h trap (ζ 0 ) = Ω trap (µ 0 , T )/Ω trap (µ 0 , T ), where ζ 0 = e −βµ 0 is the inverse fugacity at the trap center.
Using the expansion at high temperatures, we find explicitly that, This expression will be used later on in a comparison for a trapped Fermi gas at unitarity.

Comparisons for a uniform Fermi gas at unitarity
We consider now the comparison between theory and experiment for a uniform Fermi gas at unitarity. For this purpose, we calculate the universal function h(ζ) using different perturbation theories and compare the results to the experimental measurement ( Fig.  3a in Ref. [19]) This can be easily done with the known equation of state of the uniform unitarity gas, that is, where f µ and f E are two dimensionless functions that depends on the reduced temperature τ = T /T F only. Numerically, for a fixed reduced temperature τ , we calculate the inverse fugacity ζ and Ω (1) (ζ) /(Nǫ F ). We then obtain,

Perturbation theories
The universal functions h(ζ) obtained in this manner are shown in Fig. 2 for different perturbation theories, and are compared to the experimental data (red squares) [19] without any adjustable parameters. We indicate the experimentally determined superfluid regime by thin cross lines, where ζ c ≃ 0.042. The GG theory and the NSR theory above T c were already compared with experiment by Nascimbne et al. in their experimental paper [19]. The comparison presented here is much more complete. We find that the agreement between experiment and the three T -matrix perturbation theories is very good for a large range of temperatures.
In particular, the simplest Gaussian pair fluctuation theory gives the best quantitative description, though it has a possibly unphysical discontinuity at the superfluid phase transition. In the Gaussian pair fluctuation theory, the universal function does not have values between [ζ c− ] GP F ≃ 0.15 and [ζ c+ ] N SR ≃ 0.22. This non-overlap region is mostly caused by the breakdown of the original NSR theory in the vicinity of transition [30], which predicts a smaller chemical potential and hence a larger inverse fugacity ζ. In addtion, the critical inverse fugacity predicted by GPF, ζ c ≃ 0.2, is significantly larger than the experimental observation [19], ζ c ≃ 0.042. Instead, the selfconsistent GG theory and partially self-consistent GG 0 theory predict more reasonable values for ζ c , although their agreement with the experimental data of h(ζ) is worse than the GPF (NSR) theory. The prediction of available quantum Monte-Carlo simulations is also in close agreement with the experimental data [48,49]. Clearly, it would be useful to have more accurate experimental data at this point, to better understand the nature of the phase transition. On the other hand, the pseudogap theory, as a simplification of the partially selfconsistent GG 0 theory, deviates significantly from the experimental data. It is therefore not able to capture the strong fluctuations at unitarity. However, it is certainly better than the BCS mean-field theory, which completely ignores pairing fluctuations.

Virial expansion comparisons
Let us now focus on the high temperature regime of ζ > 1 or z < 1. A comparison of experimental data to the virial expansion in Eq. (22) has already been carried out by Nascimbne et al. [19]. This led to the confirmation of our theoretical prediction of the third virial coefficient ∆b 3 ≃ −0.35 as well as an experimental determination of the fourth virial coefficient ∆b 4 ≃ 0.096 ± 0.015. The third virial calculation requires an exact solution of a quantum three-body problem, which is known. However, the exact solution of the quantum four-body problem needed for the fourth coefficient, is yet to be theoretically obtained.
Here, we show that this accurate experimental data can serve as a benchmark to determine to what extent the virial expansion is quantitatively reliable. To be concrete, we shall define the criterion of "quantitative" applicability as an agreement within 10% relative error for the function h(ζ) − 1, that is, after the non-interacting background is removed from the universal function h(ζ). For a "qualitative" applicability, we relax the criterion on the relative error to 50%. Fig. 3 compares the virial expansion predictions (up to the fourth virial coefficient) for h(ζ) − 1 with the experimental data, using artificial 50% (a) and 10% (b) relative errors for comparison purposes. We are then able to estimate a critical fugacity, below which the n th -order virial expansion is either qualitatively or quantitatively valid. The result is tabulated in Table I, where the critical fugacities have also been converted to critical temperatures by using the NSR equation of state, which provides an excellent description of the experimental data. For the 2 nd virial expansion that accounts for the leading interaction effect, we determine that for a homogeneous Fermi gas at unitarity, it is quantitatively and qualitatively reliable above T /T F = 1.3 and T /T F = 0.9, respectively.
Order z 50 z 10 (T /T F ) 50  To close this section, we emphasize that, though the experimental data for the universal function h(ζ) are very accurate, without any prior knowledge of the bulk equation of state we still can not determine the temperature from the discrete data. However, for a trapped Fermi gas at unitarity, using the bulk h(ζ) we can indeed determine all the universal thermodynamic functions such as E(T ) and S(T ). This is discussed in detail in the next section. indicate respectively the 50% and 10% relative errors of h(ζ) − 1, in accord with the "qualitative" and "quantitative" criterions as described in the text.

Comparisons for a trapped Fermi gas at unitarity
Let us now turn to the experimental determination of the equation of state of a Fermi gas in a harmonic trap at unitarity, and a comparison of this data with theory. The basic idea has already been outlined by Nascimbne et al. in the Supplementary Discussion part of their paper [19].

Local density approximation
Consider, for example, the total number of atoms, N = drn(r). Within the local density approximation, we may consider the trap to be isotropic trap with a trapping frequency ω ≡ (ω 2 ⊥ ω z ) 1/3 , without loss of generality. Using n(r) = ∂P (r) /∂µ (r) and ∂µ (r) /∂r = −mω 2 r, we find that, Noting that P = P (1) h(ζ) , and integrating by parts, we obtain: Using ζ (r) = e −βµ(r) = ζ 0 exp[mω 2 r 2 /(2k B T )], the integration over radius r can be converted to an integration over the inverse fugacity. One sees that, Recalling that the Fermi energy of a zero-temperature trapped ideal Fermi gas is E F = k B T F = (3N) 1/3h ω, we may rewrite the above equation in dimensionless form, The total energy of the system may be conveniently calculated by using the scaling relation, Eq. (17). Thus, the total energy in a trap is given by Converting to the variable ζ, we find that: The entropy follows directly from the thermodynamic relation S = (E − Ω − µ 0 N)/T . Using the fact that µ 0 /E F = −(T /T F ) ln ζ 0 , we obtain straightforwardly The coupled equations (31), (33) and (34) determine, respectively, the temperature, energy, and entropy of a trapped Fermi gas at unitarity as a function of the inverse fugacity ζ 0 . In the calculations, we discretize the integral over ζ and take the values of h (ζ) solely from the experimental measured data. In this way [19], we avoid the use of any interpolating or fitting function to the experimental data h (ζ). In addition, the statistical error of the experimental data of h (ζ) is reduced. The detailed procedure for these numerical calculations is given in Appendix A. For convenience, we shall refer to this trapped equation of state, re-constructed from h (ζ), as the "experimental" measurement or data for the equations of state of a trapped Fermi gas at unitarity.
It is readily seen that we can calculate the theoretical prediction for the trapped equation of state by using exactly the same local density approximation procedure, combined with a theoretical universal function h(ζ) generated from different strongcoupling theories for a uniform Fermi gas. We note that our numerical procedure of calculating the trapped equation of state is an average procedure integrating over the trap and is quite insensitive to the smoothness of h(ζ). Therefore, even though there is a discontinuity in the theoretical universal function, as in our Gaussian pair fluctuation theory, we obtain a much smoother trapped equation of state. As can be seen from Appendix A, in that case, we simply join linearly between [ζ c− ] GP F and [ζ c+ ] N SR to remove the discontinuity of the universal function in the non-overlap region

Trapped universal thermodynamics: E(S)
In previous work [11], we gave experimental evidence that any strongly interacting Fermi gases at unitarity has universal thermodynamics. The energy and entropy relation E(S) measured on 6 Li and 40 K atomic clouds in three different trapping potentials all fall precisely on a single curve. The trapped equation of state E(S) deduced from the experimental data of h(ζ) by Nascimbne et al. [19] provides an independent check of universality, with a much improved accuracy, in a fourth different set of experimental conditions. This is illustrated in Fig. 4, where we plot the new measurements using green circles. All the four sets of experimental data follow the theoretical prediction given by the simplest GPF approximation. In particular, the difference between our theory and the new measurement is nearly indistinguishable, as shown clearly in Fig.  4b for the low temperature regime. This gives so far the strongest evidence for fermionic universality.
To better visualize the difference between theory and experiment, we follow the strategy used in our previous comparative study [38] and calculate the interaction energy E int = E − E IG , which is the difference of energies between an interacting Fermi gas (E) and an ideal Fermi gas (E IG ) at the same entropy. Fig. 5 shows the interaction energy versus entropy in a harmonic trap as predicted by different strong-coupling theories in comparison with the experimental data reported at Duke [14] and ENS [19]. It is impressive that at this much reduced scale, the new experimental data obtained by Nascimbne et al. at ENS (green empty circles) still appear to be very smooth, suggesting an absolute statistical error of about 0.01NE F or a relative error of one percent in energy, although there may be systematic errors at this level.
This error bar is already much smaller than the difference among different Tmatrix approximations of GPF, GG 0 and GG, which is roughly of the order 0.05NE F . It is clear that below threshold the GPF approach provides the closest prediction to the new accurate measurement below threshold, with a difference at most 0.01NE F . However, well above threshold the fully self-consistent GG theory gives better agreement. We note here that the above threshold GG theory is also universal, without ad-hoc renormalizations.

Trapped universal function h trap (ζ 0 )
In analogy with the comparison between theory and experiment for the universal function h(ζ) for a homogeneous Fermi gas at unitarity, we consider now the comparison for the trapped universal function h trap (ζ 0 ) = Ω trap /Ω (1) trap . Because of the proportionality between energy and thermodynamic potential at unitarity (i.e., the scaling relation), it is convenient to calculate h trap (ζ 0 ) using the expression, where the denominator is simply E (1) /(NE F ) at given inverse fugacity ζ 0 and temperature T . Fig. 6 compares the theoretical predictions for the trapped universal function from different strong-coupling theories, compared with the experimental measurement at ENS (empty squares) [19]. Both the experimental data and the theoretical GPF (NSR) predictions for the trapped universal function h trap (ζ 0 ) are now much smoother after integrating over the trap. Thus, the discontinuity of the normal-superfluid transition in the GPF prediction of h(ζ) disappears completely in h trap (ζ 0 ) , due to trap-averaging. Comparing the different T -matrix fluctuation theories in a trap, one sees that the fully self-consistent GG theory gives a slightly better agreement with the experimental results than the GPF and GG 0 theories at high temperatures, with an inverse fugacity ζ 0 > 2. At low temperatures where ζ 0 ≪ 1, however, the simplest GPF approach provides the best description.

Quantum virial comparisons: h(ζ)
We now examine the applicability of the quantum virial expansion for a trapped Fermi gas at unitarity, by using Eq. (26) for the trapped universal function. We use the same idea as in Fig. 3 and the same criterion for "qualitative" and "quantitative" reliability. In Fig. 7, we report the successive virial expansion as a function of fugacity up to the 4 th order, compared with the experimental data. The estimates of the critical fugacity and of the critical temperature for different orders are tabulated in Table II. To the leading second order, we find that the expansion is quantitatively reliable for temperatures down to T ≃ 0.7T F , much smaller than we found for a homogeneous Fermi gas at unitarity. This much wider applicability is due to the significantly reduced higher order virial coefficients in a harmonic trap, i.e., b 2,trap = b 2 /(2 √ 2). It is readily seen that, with inclusion of higher order virial coefficients, the accuracy of virial expansion can be improved. Up to the known fourth virial coefficient, we find that the bound for quantitative applicability decreases further to T ≃ 0.4T F , which is a typical experimental temperature for a Fermi gas in its normal state. We thus show that the quantum virial  Table 2. Qualitative (50%) or quantitatively (10%) ranges of reliability for different order virial expansions in trapping potentials, as indicated by the subscript.
expansion method is a very useful tool for understanding the properties of a normal, strongly interacting Fermi gas in a harmonic trap.

Thermodynamic functions E(T ) and S(T )
As we mentioned earlier, in addition to the energy-entropy relation E(S), the measurement by Nascimbne et al. [19] was able to re-construct a complete set of thermodynamic functions in harmonic traps, such as E(T ) and S(T ). This provides us with a unique opportunity for a systematic comparison between theory and experiment for a trapped Fermi gas at unitarity, without the use of any fitting functions or adjustable parameters.  Fig. 8 presents the comparison for the total energy as a function of temperature. As anticipated, the simplest GPF approach provides the best quantitative agreement with the experimental data. However, the GPF theory predicts a larger normal-superfluid transtion temperature, (T /T F ) c ≃ 0.27, as indicated by a small bump. In contrast, the less accurate self-consistent GG theory and partially self-consistent GG 0 theory predict that (T /T F ) c ≃ 0.21, which is much closer to the experimental observation of (T /T F ) c ≃ 0.19.
It is interesting to note that the GPF curve in the figure was first calculated by the present authors [30,64] and was compared to the heat capacity measurement reported by Kinast et al. [13]. However, at that time, the temperature was not independently calibrated, due to the absence of a reliable thermometry in the strongly interacting regime. An empirical temperature was used, obtained by fitting the integrated onedimensional density profile to an ideal Thomas-Fermi distribution.
In this earlier comparison, empirical temperatures were converted to actual temperatures using the pseudogap theory [13]. As a consequence, the resulting experimental data of E(T ) appeared to agree well with the pseudogap theory [13]. This is in sharp contrast to what is shown in Fig. 8, where the pseudogap theory clearly fails to account for the strong pairing fluctuations at either low temperatures (T < 0.1T F ) or high temperatures (T > 0.3T F ). We therefore conclude that while the empirical temperature approach provides a rough thermometry, its model-dependence and insensitivity to the actual temperature makes it less useful as a tool for accurate comparison of theory with experiment.   9 shows the comparison between theory and experiment for the entropy as a function of temperature. We find again that the GPF approach gives an overall best fit to the experimental data. Compared to the case of total energy, however, the effect of pairing fluctuations on the entropy is less significant. As a result, all the perturbation theories predict a similar entropy curve. Their difference to the ideal Fermi gas prediction (thin dashed line) is also small. This provides a justification for a recent calibration strategy used for determining the entropy of a weakly interacting Fermi gas [14], in which the entropy of a 6 Li cloud at a magnetic field B = 1200G is assumed to be close to that of an ideal Fermi gas.

Quantum virial expansion for E(T ) and S(T )
Finally, we may check the applicability of quantum virial expansion by using the thermodynamic functions E(T ) and S(T ). This is illustrated in Fig. 10, where the theoretical predictions of a virial expansion up to 4 th order for E(T ) and S(T ) are compared with the experimental data. We may determine directly from the figures the critical temperature related to the reliability of virial expansion. These are in agreement with the values listed in Table II.

Thermometry of a trapped Fermi gas at unitarity
We have noted that the temperature of a strongly interacting Fermi gas is difficult to measure experimentally in ultra-cold atom experiments. Unlike the situation with cryogenic experiments in the past, ultra-cold atoms are completely insulated by a high vacuum from any external reservoir at a known temperature. A useful way to quantify the temperature is to measure a non-interacting temperature T i of an ideal, non-interacting Fermi gas. This can be easily measured from the density profile, before a slow, adiabatic sweep to the Feshbach resonance. Since the entropy of a non-interacting Fermi gas is known, and is unchanged in an adiabatic sweep, this is essentially an entropy measurement. This procedure was first adopted by Regal et al. at JILA [17]. The accurate determination of the universal thermodynamic function S(T ) in the last section presents a model-independent way to re-calibrate the T i thermometry for a trapped Fermi gas at unitarity. By equating S(T ) and S IG (T i ), where S IG is the entropy of an ideal Fermi gas, we can express the temperature T of strongly interacting Fermi gases as a function of the non-interacting temperature T i . The result is reported in Fig. 11. The inset emphasizes the normal-superfluid transition regime. Note that the isentropic conversion to resonance tends to decrease the temperature so that T i is always somewhat below the temperature T at unitarity. This can be seen from the dashed line for which T = T i .
We may identify three different temperature regimes from the figure. At temperatures T i > 0.3T F , the calibration curve is nearly parallel with the equal temperature line. To a good approximation, we find that T ≃ T i + 0.04T F . Below T i = 0.3T F , the unitarity temperature seems to decrease slightly faster with decreasing non-interacting temperature. However, at a characteristic temperature (T i /T F ) c ≃ 0.16 ± 0.02, this trend changes suddenly and we observe that the unitarity temperature now tends to saturate with any further decrease of the non-interacting temperature. This interesting feature is clearly seen in the inset, where we linearly fit the data above and below (T i /T F ) c .
We can identify the sudden change as the deviations of thermodynamic properties away from normal Landau-Fermi-liquid behavior [19,48] and therefore determine a characteristic temperature (T /T F ) 0 ≃ 0.19 ± 0.02. This value agrees very well with the one that deduced from the homogeneous critical temperature by Nascimbne et al. [19] and the condensate fraction measurement by Horikoshi et al. [18]. The characteristic non-interacting temperature (T i /T F ) 0 ≃ 0.16 is also in very good agreement with the measurement at JILA [17] and the recent moment of inertia measurement at Innsbruck [65]. Using the experimental thermodynamic functions E(T ) and S(T ) in Figs. (8) and (9), we obtain (E/NE F ) 0 ≃ 0.68 and (S/Nk B ) 0 ≃ 1.56.

Conclusions and outlooks
In conclusion, from the experimental data measuring a uniform universal function h(ζ) we have deduced a complete set of universal thermodynamic functions for a trapped Fermi gas at unitarity. These accurate experimental results provide a unique opportunity to test quantum many-body theories of strongly interacting Fermi gases.
We have presented such a study by systematically comparing the theoretical predictions from typical strong-coupling theories with the experimental data. The comparison has no fitting functions or adjustable parameters. All the model approximations seem to be fluctuating around and not converging towards the accurate experimental data. We have found that the simple Gaussian pair fluctuation theory pioneered by Nozires and Schmitt-Rink [21,29] provides the best quantitative description for the universal thermodynamic properties of energy and entropy. Our comparison also includes a quantum virial expansion theory (or quantum cluster expansion theory) [54,55]. We have investigated in detail the applicability of the expansion in the quantum degenerate regime.
The experimental universal thermodynamic functions calculated in this work are extremely useful. For instance, the temperature dependence of entropy S(T ) can be used to calibrate accurately the endpoint temperatures obtained from an adiabatic sweep of the magnetic field between the ideal and strongly interacting regimes. Therefore, by measuring an ideal Fermi gas temperature before the sweep and using the curve T (T i ) shown in Fig. 10, one can solve the troublesome thermometry problem for a strongly interacting Fermi gas. From these universal thermodynamic functions, we are also able to determine a characterstic temperature (T /T F ) 0 ≃ 0.19 or (T i /T F ) 0 ≃ 0.16 for a trapped Fermi gas at unitarity, which is responsible for the deviations of thermodynamic properties away from normal Landau-Fermi-liquid behavior due to pairing effects. At these points, our analysis of the experimental universal function h(ζ) provides a new insight on the superbly precise experimental work of Nascimbne et al. [19].
These thorough comparisons between theory and experiment provide a motivation for further developing the challenging many-body theory of a strongly interacting Fermi gases. It is impressive that the simplest Gaussian pair fluctuation approach gives such excellent agreement with the experimental data, especially in the below threshold regime characterized by long-range superfluid order. Yet, it fails to predict the correct normalsuperfluid transition temperature. More work is need to understand the reason for this, but at this stage, we feel that the GPF approximation serves as a good starting point for further theoretical work. Recalling that the GPF approximation includes only the two-body correlations (see, for example, Fig. 1a), a natural way to extend this may be to consider three-body or four-body correlations, in which three or four fermions interact with each other in the scattering process. The thermodynamic potential may be worked out with inclusion of all three-body or four-body scattering matrices. In this manner, we would recover correctly the equation of state predicted by the higher-order (i.e., 3 nd and 4 th order) quantum virial expansion theory at high temperatures. We believe that it will also lead to a more reasonable critical temperature and remove the spurious bend-back structure close to T c as shown in the Gaussian pair fluctuation theory.
On the other hand, the quantum virial expansion gives us another means for theoretical development, from a very different point of view. We have already shown the wide applicability of this expansion for the equation of state of a strongly interacting Fermi gas in harmonic traps, down to temperatures as low as ∼ 0.4T F . We have conjectured that it may be applicable down to the superfluid transition temperature, with inclusion of higher-order virial coefficients. Here, we can also develop a virial expansion for more crucial dynamical properties, such as the dynamic structure factor and single-particle spectral function, as measured recently at Swinburne [66] and at JILA [67]. The quantum virial expansion may therefore solve the troublesome problem of understanding a normal yet strongly interacting Fermi gas, although theoretical predictions beyond third order are not yet available.
Finally, we note that there is huge interest in determining the detailed behavior of a strongly interacting Fermi gas near the normal-superfluid transition. Our comparative study, based on the most recent theoretical and experimental results, may provide useful insights for this future research.