The equation of state of ultracold Bose and Fermi gases: a few examples

We describe a powerful method for determining the equation of state of an ultracold gas from in situ images. The method provides a measurement of the local pressure of an harmonically trapped gas and we give several applications to Bose and Fermi gases. We obtain the grand-canonical equation of state of a spin-balanced Fermi gas with resonant interactions as a function of temperature. We compare our equation of state with an equation of state measured by the Tokyo group, that reveals a significant difference in the high-temperature regime. The normal phase, at low temperature, is well described by a Landau Fermi liquid model, and we observe a clear thermodynamic signature of the superfluid transition. In a second part we apply the same procedure to Bose gases. From a single image of a quasi ideal Bose gas we determine the equation of state from the classical to the condensed regime. Finally the method is applied to a Bose gas in a 3D optical lattice in the Mott insulator regime. Our equation of state directly reveals the Mott insulator behavior and is suited to investigate finite-temperature effects.


Introduction
Ultracold gases are a privileged tool for the simulation of model Hamiltonians relevant in the fields of condensed matter, astrophysics or nuclear physics in the laboratory [3]. As an example, thanks to the short-range character of interactions, ultracold Fermi mixtures prepared around a Feshbach resonance mimic the behavior of neutron matter in the outer crust of neutron stars [4,5]. For cold atoms, the density inhomogeneity induced by the trapping potential has long made the connection between the Hamiltonian of a homogeneous system and an ultracold gas indirect. Early experimental thermodynamic studies have provided global quantities averaged over the whole trapped gas, such as total energy and entropy [6,7], collective mode frequencies [8] or radii of the different phases that may be observed in an imbalanced Fermi gas [9]- [11]. Reconstructing the equation of state of the homogeneous gas then requires deconvolving the effect of the trapping potential, a delicate procedure that has not been done so far. However, the gas can often be considered as locally homogeneous (local density approximation (LDA)), and careful analysis of in situ density profiles can directly provide the equation of state of a homogeneous gas [1], [12]- [14]. In the case of two-dimensional (2D) gases, in situ images taken along the direction of tight confinement obviously give access to the surface density [15]- [18] and thus to the equation of state [19]. For three-dimensional (3D) gases, imaging leads to an unavoidable integration along the line of sight. As a consequence, inferring local quantities is not straightforward. Local density profiles can be computed from a cloud image using an inverse Abel transform for radially symmetric traps [20]. A more powerful method was suggested in [13] and implemented in [1,14]: as explained below, for a harmonically trapped gas, the local pressure is simply proportional to the integrated in situ 3 absorption profile. Using this method, the low-temperature superfluid equation of state for balanced and imbalanced Fermi gases was studied as a function of interaction strength [1,14]. In this paper, we describe in more detail the procedure used to determine the equation of state of a spin-unpolarized Fermi gas in the unitary limit [1]. We compare our data with recent results from the Tokyo group [2], and show a significant discrepancy in the high-temperature regime. In the second part, we apply the method to ultracold Bose gases. From an in situ image of 7 Li, we obtain the equation of state of a weakly interacting Bose gas. Finally, analyzing the experimental profiles of a Bose gas in a deep optical lattice [21], we observe clear thermodynamic signatures of the Mott insulator phases.

Measurement of the local pressure inside a trapped gas
In the grand-canonical ensemble, all thermodynamic quantities of a macroscopic system can be derived from the equation of state P = f (µ, T ) relating the pressure P to the chemical potential µ and the temperature T . P can be straightforwardly deduced from integrated in situ images.
Consider first a single-species ultracold gas, held in a cylindrically symmetric harmonic trap whose frequencies are labeled ω x = ω y ≡ ω r in the transverse direction and ω z in the axial direction. Provided that the LDA is satisfied, the gas pressure along the z-axis is given by [13] where n(z) = dx dy n(x, y, z) is the doubly integrated density profile, µ z = µ 0 − 1 2 mω 2 z z 2 is the local chemical potential on the z-axis and µ 0 is the global chemical potential. n(z) is obtained from an in situ image taken along the y-axis, by integrating the optical density along the x-axis (see figure 1). As described below, if one independently determines temperature T and chemical potential µ 0 , then each pixel row of the absorption image at a given position z provides an experimental data point for the grand-canonical equation of state P(µ z , T ) of the homogeneous gas. The large number of data obtained from several images allows one to perform an efficient averaging, leading to a low-noise equation of state.
This formula is also valid in the case of a two-component Fermi gas with equal spin populations if n(z) is the total integrated density. The method can be generalized to multicomponent Bose and Fermi gases, as first demonstrated on spin-imbalanced Fermi gases in [1,14].

Thermodynamics of a Fermi gas with resonant interactions
In this section, we describe the procedure used in [1] to determine the grand-canonical equation of state of a homogeneous and unpolarized Fermi gas with resonant interactions (a = ∞). We also compare our data with recent measurements from the Tokyo group [1,2]. We then study the physical content of the equation of state at low temperature.

Grand-canonical equation of state
In the grand-canonical ensemble, the equation of state of a spin-unpolarized Fermi gas in the unitary limit can be written as where P (0) (µ, T ) is the pressure of a non-interacting two-component Fermi gas and ζ = exp(−µ/k B T ) is the inverse fugacity. Since P (0) (µ, T ) is known, the function h T (ζ ) completely determines the equation of state P(µ, T ). Let us now describe the procedure used to measure it. The pressure profile of the trapped gas along the z-axis is directly derived from its in situ image using equation (1). The effect of the trap anharmonicity of the optical dipole trap on the pressure measurement is expected to be less than 5%. One still has to know the value of the temperature T and the global chemical potential µ 0 in order to infer h T (ζ ).
We use a small number of 7 Li atoms, at thermal equilibrium in the 6 Li component, as a thermometer. We then extract µ 0 from the pressure profile, by comparison in the cloud's wings with a reference equation of state. For high-temperature clouds (k B T > µ 0 ), we choose µ 0 so that the wings of the pressure profile match the second-order virial expansion [22] (see figure 2(a)): For colder clouds, the signal-to-noise ratio is not good enough, in the region where (3) is valid, to extract µ 0 using the same procedure. We thus rather use the equation of state determined from all images previously treated as a reference, since it is accurate over a wider parameter range than (3) (see figure 2(b)). We then iterate this procedure at lower and lower temperatures, eventually below the superfluid transition. By gathering the data from all images and statistical averaging, we obtain a low-noise equation of state in the range 0.02 < ζ < 5 (see figure 3(a)).
. A wrong choice of µ 0 in this representation corresponds to a translation of the data in abscissa. We adjust µ 0 so that the wings of the pressure profile match a reference equation of state (in red). (a) For high-temperature clouds, we use the secondorder virial expansion (3). (b) For a lower temperature pressure profile, we minimize its distance with the averaged equation of state deduced from higher temperature images (in red) in the overlap region.

Canonical equation of state
In [2], a canonical equation of state E(n, T ) expressing energy E as a function of density and temperature was measured using fits of absorption images taken after a short time-of-flight. In situ density profiles were deduced by assuming a hydrodynamic expansion. The temperature was extracted from the cloud's total potential energy at unitarity, using the experimental calibration made in [7]. In figure 3(b), data from [2] are plotted as E(n, T )/E (0) (n, T ) as a function of θ = T /T F , where n is the total atom density, T F is the Fermi temperature and E (0) (n, T ) is the energy of a non-interacting Fermi mixture.
The comparison between the two equations of state requires expressing our data in the canonical ensemble. The density n = ∂ P/∂µ| T is calculated by taking a discrete derivative, and we obtain the black points in figure 3(b). While the two sets of data are in satisfactory agreement in the low-temperature regime T /T F < 0.4, they clearly differ in the high-temperature regime. The disagreement of the data from [2] with the second-and third-order virial expansions calculated in [22,23] indicates a systematic error in this regime. This is possibly due to a breakdown of hydrodynamics during the time-of-flight as expected at high temperature.

Fermi liquid behavior in the normal phase
Above the superfluid transition and in the low-temperature regime 0.05 < ζ < 0.5, our data are well modeled by a Fermi liquid equation of state where ξ n = 0.51(1) and m * = 1.12(3)m respectively characterize the compressibility of the normal phase extrapolated to zero temperature and the effective mass of the low-lying excitations. The agreement with (4) is better than 5% in a large parameter range 0.33 µ < k B T < 2 µ. Our value of ξ n is in agreement with the variational fixed-node Monte Carlo calculations ξ n = 0.54 in [24], ξ n = 0.56 in [25], and with the quantum Monte Carlo calculation ξ n = 0.52 in [26]. It is surprising that the quasi-particle mass m * is quite close to the free fermion mass, despite the strongly interacting regime. Note also that this mass is close to the effective mass m * = 1.20 m of a single spin-down atom immersed in a Fermi sea of spin-up particles (the Fermi polaron) [1,11,12,25], [27]- [30].

Superfluid transition
The deviation of the experimental data from (4) for ζ < 0.05 signals the superfluid phase transition. This transition belongs to the U (1) universality class, and the critical region is expected to be wide [31] in the unitary limit. Assuming that our low-temperature data belong to the critical region, we fit our data with a function where H is the Heaviside function and α −0.013 is the specific heat critical exponent, measured with a very good accuracy on liquid 4 He [32]. We obtain the position of the superfluid transition ζ c = 0.05, or k B T c /µ = 0.33, in agreement with the value k B T c /µ = 0.32(3) extracted in [1] using a simpler fit function. We thus confirm more rigorously our previous determination of the superfluid transition. In the appendix, we discuss the validity of LDA around the superfluid transition. Under our current experimental conditions, the deviation from LDA is very small.

Thermodynamics of a weakly interacting Bose gas
In this section, we apply equation (1) to the case of trapped Bose gases. Firstly, we test the method by determining the equation of state of a weakly interacting Bose gas [33,34]. We use an in situ absorption image of a 7 Li gas taken from [35] (see figure 4(a)). 7 Li atoms are polarized in the internal state |F = 1, m F = −1 , and held in an Ioffe-Pritchard magnetic trap with ω r /2π = 4970 Hz and ω z /2π = 83 Hz, in a bias field B 0 2 G. The anharmonicity of this magnetic trap is negligible. Thermometry is provided by a gas of 6 Li atoms, prepared in |F = 1 2 , m F = − 1 2 , and in thermal equilibrium with the 7 Li cloud.

Determination of the equation of state
The equation of state of a weakly interacting Bose gas can be expressed, in the grand-canonical ensemble, as where ζ = e −µ/k B T is the inverse fugacity and λ dB (T ) = 2πh 2 /mk B T is the thermal de Broglie wavelength. The pressure profile is calculated using (1). We aim here at measuring g(ζ ). We obtain the global chemical potential value µ 0 = 0.10 k B T by fitting the 7 Li profile in the noncondensed region |z| > 50 µm with a Bose function: Combining the measurement of the pressure profile, the cloud's temperature T and the global chemical potential µ 0 , we obtain the thermodynamic function g(ζ ) plotted in figure 4(b).

Analysis of the equation of state
In the region ζ > 1, the data agree with the Bose function g(ζ ) = g 5/2 (ζ ) expected for a weakly interacting Bose gas. The departure from the thermodynamic function of a classical gas g(ζ ) = ζ −1 , and especially the fact that g(ζ ) > 1 above the condensation threshold, is the thermodynamic signature of a bosonic bunching effect, as observed in [36]- [38]. The sudden and fast increase of our data for ζ 1 indicates the Bose-Einstein condensation threshold. In the LDA framework, the chemical potential of a weakly interacting Bose-Einstein condensate reads as follows: where m 7 is the 7 Li atom mass and a 77 is the scattering length describing s-wave interactions between 7 Li atoms. We neglect thermal excitations in the condensed region. Integrating the Gibbs-Duhem relation at a fixed temperature dP = ndµ between the condensation threshold ζ c and ζ < ζ c , and imposing continuity at ζ = ζ c , we obtain the equation of state in the condensed phase: Fitting our data with the function g(ζ ) given by (6) for ζ < ζ c and with g 5/2 (ζ ) for ζ > ζ c , we obtain ζ c = 1.0(1) and a 77 = 8(4) a 0 = 0.4(2) nm. The uncertainties take into account the fit uncertainty and the uncertainty related to the temperature determination. The condensation threshold is in agreement with the value ζ c = 1 expected for an ideal Bose gas, the mean-field correction being of the order of 1% [39,40]. Our measurement of the scattering length is in agreement with the most recent calculations a 77 = 7(1) a 0 [41]. Extending this type of measurement to larger interaction strength Bose gases prepared close to a Feshbach resonance would reveal more complex beyond-mean-field phenomena, provided thermal equilibrium is reached for strong enough interactions.

Mott insulator behavior of a Bose gas in a deep optical lattice
Here we extend our grand-canonical analysis to the case of a 87 Rb gas in an optical lattice in the Mott insulator regime. By comparing experimental data with advanced Monte Carlo techniques, it has been shown that in many circumstances the LDA is satisfied in such a system [42]. We analyze the integrated density profiles of the Munich group (see figure 2 of [21]).

Realization of the Bose-Hubbard model with ultracold gases
Atoms are held in a trap consisting of the sum of a harmonic potential V h (x, y, z) and a periodic potential, V 0 (sin 2 (kx) + sin 2 (ky) + sin 2 (kz)), created by three orthogonal standing waves of red-detuned laser light at the wavelength λ = 2π/k = 843 nm. The atoms occupy the lowest Bloch band and realize the Bose-Hubbard model [43]:Ĥ with a local chemical potential µ(r) = µ 0 − V h (r). The index i refers to a potential well at position r i , J is the tunneling amplitude between nearest neighbors, and U is the on-site interaction, U and J being a function of the lattice depth [3]. The slow variation of V h (r) compared with the lattice period λ/2 justifies the use of LDA.
We consider here the case of a large lattice depth V 0 = 22E r , for which J 0.003 U ∼ 0, and assume that the temperature is much smaller than U . In this regime, the gas is expected to form a Mott insulator: in the interval µ ∈ [( p − 1)U, pU ], where p is an integer, the atom number per site remains equal to p, and the density is equal to n = p(2/λ) 3 . Integrating the Gibbs-Duhem relation between 0 and µ, we obtain that the pressure P is a piecewise linear function of µ:

Determination of the equation of state
We use a series of three images from [21], labeled a, b and c, with different atom numbers N a = 1.0 × 10 5 , N b = 2.0 × 10 5 and N c = 3.5 × 10 5 (see figure 5(a)). The integrated profiles n(z) are not obtained using in situ absorption imaging but rather using a tomographic technique, providing ∼1 µm resolution. The pressure profile is then obtained using equation (1). Each image i = a, b and c plotted as P as a function of − 1 2 mω 2 z z 2 provides the equation of state P(µ) translated by the unknown global chemical potential µ 0 i . By imposing that all images correspond to the same equation of state (in the overlapping µ/U region), we deduce the chemical potential differences between the different images µ 0 b − µ 0 a = 0.56 U and µ 0 c − µ 0 b = 0.61 U (see figure 5(b)). Gathering the data from all images, we thus obtain a single equation of state, translated by µ 0 a , which is still unknown. We fit these data with a function translated by µ 0 a from the following function, capturing the Mott insulator physics: with µ 0 a , δµ 1 , δµ 2 , n 1 , n 2 and n 3 as free parameters. The value µ 0 a = 1.51 U yielded by the fit thus corresponds to the condition P → 0 when µ → 0. Once it is determined, we obtain the equation of state of the Bose-Hubbard model in the Mott regime, plotted in figure 6.

Observation of Mott insulator behavior
After fitting the value of µ 0 a , the other parameters resulting from the fit exhibit the characteristic features of incompressible Mott phases. The occupation number in the first Mott region is n 1 = 0.9(1) atom per site and the size is δµ 1 = 0.9(1)U . The second Mott region occupation number is n 2 = 2.0(1) and its size is δµ 2 = 1.1(1)U . Finally, the third Mott region occupation number is n 3 = 3.1 (1). These values agree with the theoretical values n i = i and δµ i = U , in the T = 0 and J = 0 limits.

Estimation of finite-temperature effects
The equation of state deduced from the experimental data is also suited for investigating finitetemperature effects. Since sites are decoupled in the regime J U , k B T considered in this study, the finite-temperature equation of state is easily calculated from the thermodynamics of a single site [44,45]: Fitting now the experimental data with (8) and T and µ 0 a as free parameters, we deduce k B T = 0.09 +0.04 −0.09 U. This value is in agreement with a direct fit of the density profiles and number statistics measurements [46]. Firstly, this temperature is significantly smaller than the temperature k B T * 0.2 U at which the Mott insulator is expected to melt [44]. Secondly, this temperature should be considered as an upper limit because of its uncertainty on the low-temperature side. Indeed, the finite resolution of the images tends to smear out the sharp structure associated with Mott insulator boundaries, leading to an overestimation of the actual temperature. To overcome this limit, the spin-gradient thermometry proposed in [47] could be employed.

Summary and concluding remarks
To summarize, we have shown on various examples of Fermi and Bose gas systems how in situ absorption images can provide the grand-canonical equation of state of the homogeneous gas. This equation of state is obtained up to a global shift in chemical potential and we have given several examples for its determination. The method relies on the LDA, which is satisfied in many situations, but notable exceptions exist such as the case of the ideal Bose gas. The equation of state given by this procedure allows a direct comparison with many-body theories. Although we have here illustrated this method on a single-component Bose gas and a spin-balanced Fermi gas, it can easily be generalized to multi-component gases. For instance, the phase diagram and the superfluid equation of state of spin-imbalanced Fermi gases were obtained in [1,14]. We expect this method to be very useful in the investigation of Bose-Bose, Bose-Fermi and Fermi-Fermi mixtures. Finally, the equation of state of a Bose gas close to a Feshbach resonance may reveal thermodynamic signatures of beyond-mean-field behavior in Bose-Einstein condensates [48].

Appendix. Validity of local density approximation (LDA)
Let us now discuss the validity of LDA around the superfluid transition in our experiment. Along the z-axis, the correlation length ξ diverges around the transition point z = z c according to ξ ∼ k −1 F |(z − z c )/z c | −ν , where ν = 0.67 is the correlation length critical exponent, directly measured in [49], and in agreement with ν = (2 − α)/3. LDA is expected to become inaccurate in the region z c − δz < z < z c + δz, where δz is given by [31,50] δz ∼ ξ(z c + δz), i.e. δz ∼ z c (k F z c ) −1/(1+ν) .
z c is of the order of the cloud size along z, and is much larger than k −1 F , which is of the order of the inter-particle distance. Given the parameters of our experiments, (k F z c ) −1/(1+ν) ∼ 1% and the size δz where LDA is invalid is very small. Given the noise of our data (a few per cent), the deviation from LDA is thus negligible. Investigating the critical behavior at the superfluid transition, such as measuring the critical exponent α, would be an interesting development for this method, as proposed in [50].