Ground-state thermodynamic quantities of homogeneous spin-$1/2$ fermions from the BCS region to the unitarity limit

We experimentally determined various thermodynamic quantities of interacting two-component fermions at the zero-temperature limit from the Bardeen-Cooper-Schrieffer (BCS) region to the unitarity limit. The obtained results are very accurate in the sense that the systematic error is within 4% around the unitarity limit. Using this advantage, we can compare our data with various many-body theories. We found that an extended ${\it T}$-matrix approximation, which is a strong-coupling theory involving fluctuations in the Cooper channel, well reproduces our experimental results. We also found that the superfluid order parameter ${\it \Delta}$ calculated by solving the ordinary BCS gap equation with the chemical potential of interacting fermions is close to the binding energy of the paired fermions directly observed in a spectroscopic experiment and that obtained using a quantum Monte Carlo method. Since understanding the strong-coupling properties of a superfluid Fermi gas in the BCS-BEC (Bose-Einstein condensation) crossover region is a crucial issue in condensed matter physics and nuclear physics, the results of the present study are expected to be useful in the further development of these fields.


I. INTRODUCTION
A many-body system of fermions interacting with swave scattering length a is a fundamental model that extends the ideal Fermi gas model for various interacting Fermi systems. Understanding the ground state properties of fermions in the unitary regime, where the absolute value of the scattering length is larger than the interparticle distance d ∼ k −1 F , is crucial in condensed matter physics and nuclear physics, where k F is the Fermi wave number. For example, the relation between number density n and the internal energy density E with neutronneutron scattering length a NN = −18.63 fm [1] gives the equation of state (EOS) for dilute neutron matter. The EOS characterizes the symmetry energy of nucleons, the neutron skin for neutron-rich nuclei, and the inner structure of neutron stars [2]. Since neutrons have a negative scattering length, dilute neutron matter exists in the interaction range from the Bardeen-Cooper-Schrieffer (BCS) limit to the unitarity limit at different neutron densities. Therefore, many-body physics from the BCS limit to the unitarity limit is common to both condensed matter physics [3] and nuclear physics [4].
Ultracold atomic gases provide an ideal research environment in which we can investigate such manybody Fermi systems universally and systematically [5][6][7]. Many-body systems composed of ultracold atoms have an ideal hierarchy of length and energy, where the interparticle distance and the wavelength of matter are suffi-ciently long compared to the size of the short-range interaction potential. Furthermore, it is possible to tune the scattering length between two fermions using Feshbach resonances [8]. These features realize universal many-body systems, where various physical phenomena are independent of details of particles [9].
The ground state of many-body fermions interacting with an s-wave scattering length is an s-wave superfluid state, which has a nonzero superfluid order parameter ∆ of the paired fermions and fluctuations of the order parameter. The origin of such superfluid fluctuations is the repetition of pair formations and their dissociations as well as non-condensed pairs that are kicked out from the condensate. Far below the superfluid phase transition temperature T c , where thermal excitations are almost absent, many-body corrections to physical quantities are dominated by quantum fluctuations associated with superfluid fluctuations.
At the BCS limit, the BCS mean-field (MF) wavefunction has been considered to be an adequate approximation to describe the ground state properties because the influence of quantum fluctuations appears to be small due to weak interaction between fermions. In this BCS-MF approximation, the magnitude of the order parameter ∆ is equivalent to the binding energy of a paired fermion and to the energy gap in single-particle excitations [10]. The relation between the thermodynamic quantities and ∆ is given by the gap equation and the number equation. The influence of quantum fluctuations on the ground state energy has also been theoretically analyzed at the BCS limit. The thermodynamic quantities of such fermions have been expressed analytically up to the order of (k F a) 2 beyond the MF approximation [11].
In the unitary regime, on the other hand, it is not obvious whether the BCS-MF wave function is an adequate approximation for describing the ground state properties because considerable quantum fluctuations are induced by strong interaction between fermions. There are no analytical formulas to give the order parameter, the binding energy of a paired fermion, the thermodynamic quantities, or the relations among them. Thus far, the binding energy of a paired fermion [12], a single-particle excitation spectrum [13][14][15], and internal energy density E [16], which is the ground state energy per unit volume, have been experimentally determined in the unitary regime. Other thermodynamic quantities, such as pressure [16,17], isothermal compressibility [17,18], speed of sound [19], and contact density [15] were also measured. However, the magnitude of the order parameter and chemical potential µ have not yet been determined experimentally in the unitary regime. While several theories have demonstrated qualitative agreement with these experimental results, quantitative evaluation of these theories has not been achieved because of experimental accuracies, lack of data points, and inhomogeneity effects.
In the present study, we determine the relations among thermodynamic quantities, such as pressure P , number density n, internal energy density E, chemical potential µ, isothermal compressibility κ, and contact density C, for homogeneous spin-1/2 superfluid fermions at the zerotemperature limit from the BCS region to the unitarity limit. All of these quantities were determined within 4% systematic error around the unitarity limit without using any model functions. We provide their thermodynamic functions in universal form as a function of an interaction parameter. Therefore, these functions are applicable to other similar many-body Fermi systems, such as dilute neutron matter. We compared our experimental data with various many-body theories and found that one theoretical model can reproduce our data. We also calculated the superfluid order parameter by solving the ordinary BCS gap equation with the determined chemical potential. We compare the obtained results with the binding energies of paired fermions, which were determined in previous studies through a spectroscopic experiment and quantum Monte Carlo calculation.
The remainder of the present paper is organized as follows. In Section II, we introduce the theoretical framework used in our data analysis. In Section III, we explain our experimental procedure, data analysis, and experimental results. In Section IV, we compare our experimental results with the results of various previous experimental and theoretical studies and discuss calculation of the superfluid order parameter. Finally, in Section V, we conclude this study.

II. THEORETICAL FRAMEWORK
A. Thermodynamics of homogeneous fermions from the BCS limit to the unitarity limit at zero temperature We consider homogeneous spin-1/2 fermions at T = 0, which have a number density of n ↑ = n ↓ = n/2 and a chemical potential of µ ↑ = µ ↓ = µ for each spin state. We assume that the interaction between fermions in different spin states is modeled by the s-wave scattering length a. In the following, for convenience, we use a −1 instead of a. We consider the parameter region a −1 ≤ 0, where the Fermi system ranges from the BCS limit (a −1 = −∞) to the unitarity limit (a −1 = 0). Let m be the mass of the fermions, and let be the reduced Plank constant.
The internal energy density E(n, a −1 ) is a function of n and a −1 . Let us write its differential form as The thermodynamic quantity C defined in the above relation is called the contact density [20][21][22][23][24][25]. The pressure P is derived from E(n, a −1 ) as The functional dependence P (µ, a −1 ) has a simple differential form: Isothermal compressibility κ and speed of sound v have the same thermodynamic relations as an ideal Fermi gas under a fixed scattering length:

B. Dimensionless thermodynamic functions
Here, we use the fact that the theory involves only two constants, m and . If we change the unit of length by a factor of λ and that of time by λ 2 , the values of the two constants remain unchanged. Hence, the functional dependence of any dimensionless qualities A on n and a −1 should satisfy A(n, a −1 ) = A(λ 3 n, λa −1 ). This implies that A can be written as a function A(x) of a parameter x proportional to a −1 n −1/3 [9]. Similarly, when a dimensionless quantity B is a function of µ and a −1 , it should be written as a function B(X) of a parameter X proportional to a −1 µ −1/2 .
In the present paper, we choose the two parameters to be grand-canonical interaction parameter: canonical interaction parameter: Here, we define two wave numbers as k µ (µ) = √ 2mµ/ and k F (n) = 2mε F (n)/ , where ε F (n) = 2 2m (3π 2 n) 2/3 is the Fermi energy for ideal fermions. For each thermodynamic quantity, we define an associated dimensionless thermodynamic function by normalizing it with the value for an ideal Fermi gas as follows: number density: f n (X) = n n 0 (µ) , internal energy density: chemical potential: isothermal compressibility: speed of sound: where P 0 (µ) = nεF (n) , and v 0 (n) = (mnκ 0 (n)) −1/2 . An exception is the contact density C, which vanishes for the ideal gas. Here, we choose its dimensionless thermodynamic function as In an experiment with interacting Fermi gas in a trap potential in Section III, we obtain data including various combinations of (µ, P ) with a fixed value of the scattering length a. It is convenient for data analysis to introduce dimensionless versions of µ and P by using normalizing factors solely determined from a as where we define energy as ε a = 2 2ma 2 . According to the dimensionless analysis discussed above, we can see that G is given as a function of χ, i.e., G(χ). Since X and f P are simply related to G(χ) as G(χ) can be directly converted to f P (X). The function f n (X) is derived from f P (X) as The function f E (x) is constructed by The When f µ (x) is given, f κ (x) can be determined by The function f v (x) has a simple relation to f κ (x), The function f C (x) is obtained from the derivative of f E (x) as, Equations (17) through (20) indicate that there is a universal value ξ at the unitarity limit (X = x = 0), Note that ξ is sometimes called the Bertsch parameter [5].
C. Behavior of the dimensionless functions at the BCS limit The asymptotic behavior of f E (x) at the BCS limit is known up to the second order of 1/x [11,26] as Here, we omit the condensation energy term [11]. Moreover, f Asym µ (x) is given by Eq. (19) with f Asym E (x), and f Asym P (X) is obtained using the following relations: The asymptotic behavior of the other thermodynamic quantities at the BCS limit can be obtained from f Asym P (X) using Eqs. (17)  The z axis is defined as the axial direction of the optical trap.

A. Data acquisition
We have an experimental apparatus that can produce dual Bose-Einstein condensates (BECs) of paired 6 Li (fermion) and one spin state of 7 Li (boson) in the hybrid trap of an optical dipole trap (ODT) and a magnetic trap (MT) [27]. A schematic drawing of the experimental apparatus is shown in Figure 1. Detailed information about our experimental apparatus and experimental procedure is provided in a previous study [27]. In this experiment, only 6 Li atoms were used. The shape of the trapping potential U trap is analytically given without using a harmonic approximation in terms of the parameters of the laser power, the wavelength, the beam waists of the ODT, and the magnetic curvature of the MT. The trap has an elliptic symmetry and can be written The ODT is produced by a focused 1070-nm laser beam.
The 1/e 2 beam radii are (w 0x , w 0y ) = (43.5, 46.9) µm, and the final laser power after evaporative cooling is P ODT =45 mW. The depth of the ODT is 38 uK, and the trapping frequency realized by the ODT is (ω x , ω y , ω z ) = 2π × (250, 230, 1.3) Hz. The MT is produced by the magnetic curvature of a bias magnetic field for the Feshbach resonance and produces magnetic curvature ω mag = 2π × 0.24 √ B Hz in the z direction, where B is the produced bias magnetic field in Gauss. The effective trapping frequency along the z direction is given by ω z,eff = ω 2 z + ω 2 mag . Typically, 2 × 10 5 degenerate fermions are prepared in this experiment. 6 Li atoms have a broad s-wave Feshbach resonance for collisions between the two lowest spin states: |1 ≡ |m L = 0, m S = −1/2, m I = +1 and |2 ≡ |m L = 0, m S = −1/2, m I = 0 , where m L , m S , and m I are, respectively, the electronic orbital angular momentum projection, the electronic spin projection, and the nuclear spin projection. The scattering length is given by a(B) = a bg 1 + Wres B−Bres with parameters of a bg =−1582a 0 , B res =832.18 Gauss, and W res =262.3 Gauss, where a 0 is the Bohr radius [28]. Then, the Fermi system is in the BEC region (a > 0) at B < B res . On the other hand, the Fermi system is in the BCS region (a < 0) at B > B res . At a(B = B res ) = ±∞, the Fermi system reaches the unitarity limit.
In the present study, we are interested in the ground state properties of the interacting fermions from the BCS limit to the unitarity limit. In order to prepare highly degenerate Fermi gas (t = k B T /ε F (n) ≪ 1) in the interaction region, we produced an almost pure molecular BEC consisting of two fermions in the state of |1 and |2 at 777 Gauss in the BEC regime. We then changed the scattering length adiabatically by sweeping the magnetic field to 832.18∼1050 Gauss in 1.5 s.
We can accurately obtain their density distribution in situ using the imaging techniques developed in a previous study [29]. We measured the optical depth (OD) of the trapped fermions at each magnetic field, as shown in Figure 2a. From the OD and the absorption crosssection σ abs , we evaluated the column densityn(x, z) = +∞ −∞ n(x, y, z)dy byn(x, z) = OD(x, z)/σ abs . We determined the effective value of σ abs within 4% uncertainty [29]. Note that this uncertainty results in one of the systematic errors in our final results.
The local pressure of the trapped fermions can be calculated from the column densityn(x, z) and the trapping potential U trap (ρ, z) using the following formula [29][30][31]: which is derived from the Gibbs-Duhem equation [Eq.
(3)], the local density approximation (LDA), i.e., µ = µ 0 − U trap , and the inverse Abel transformation of the column density, i.e., n(ρ, z) ∂x dx. This calculation can be carried out without the knowledge of the global chemical potential µ 0 of the trapped system. Finally, we can obtain pressure P as a function of the trapping potential U trap by relating P (ρ, z) and U trap (ρ, z) at each position. Figure 2b shows a typical example of P (U trap ) at 882.86 Gauss. The data are averaged at 241 potential  We extracted P (U trap ) and n(U trap ) in the highly degenerate region, where the number density is larger than half of the peak density for the following data analysis, because t = k B T /ε F (n) ∝ n −2/3 . The red curves plotted in Figures 2b and 2c show the regions. We acquired data sets of {P (U trap ), n(U trap ), a(B) −1 } at 44 magnetic fields of from 832.54 to 1,050 Gauss for the BCS region and of 832.18 Gauss for the unitarity limit. The 44 magnetic fields were chosen in order that the canonical interaction parameter x n(U trap ), a(B) −1 obtained at B overlaps one another. In this experiment, we repeated data acquisition eight times at each magnetic field to decrease statistical errors.
We hereinafter express physical quantities at B by adding superscript B, for example, P B , n B , and U B trap .

B. Evaluation of fκ(x)
The dimensionless isothermal compressibility f κ (x) can be constructed model-independently from data sets trap , B) from these sets using the thermodynamic relation of Eq. (4) and the definitions of x(n, a −1 ) and κ 0 (n). We collected f B κ (x B ) from the BCS region to the At the unitarity limit, our data indicate that f κ (0) = 2.66 (7), which corresponds to the universal value of ξ = 0.375 (10). This value is very consistent with the accurate experimental value of 0.376(4) determined by measuring thermodynamics for homogeneous unitary Fermi gases [17]. In the BCS region, our data approach the theoretical asymptotic behavior of f κ (x).
Experimentally, it is impossible to prepare fermions at zero temperature, even if we start from a molecular BEC. Therefore, experimental temperature T exp always has a nonzero value, i.e., t exp > 0. Since the superfluid critical temperature t c decreases monotonically to zero from the unitarity limit to the BCS limit as a function of the interaction parameter x, it is inevitable that t c (x) intersects t exp somewhere between the two interaction limits. Here, we define the transition point as t exp = t c (x * ).
Based on another experiment, as described in Appendix A, the transition point was determined to be x * = −1.14, and the temperature parameter of the fermions was estimated to be 0.06 t exp 0.1 for our experimental condition. We indicate the superfluid tran-sition point x * with a vertical dash-dot line in Figure 3. In the region of x > x * , the fermions are in the superfluid state. A cusp in f κ (x) appears around the transition point, as shown in Figure 3, which is similar to that observed in a previous study [17].
According to the previous experiment [17] and a theoretical calculation [32], our temperature range of 0.06 t exp 0.1 is low enough that we can consider our Fermi system as the ground state at the unitarity regime, because a large amount of fermions are in the superfluid state, which has zero entropy.

C. Construction of G(χ)
In order to determine the other thermodynamic quantities, it is necessary to determine the chemical potential µ or contact density C along with P and n. In this experiment, we apply an almost homogeneous magnetic field to the gas, which implies that the Fermi gas has a homogeneous scattering length. In this case, it is impossible to determine C from the gradient of P relative to a −1 under the same µ for each data.
The local chemical potential µ B is given by trap under the LDA, where µ B 0 is the global chemical potential at magnetic field B. If we have experimental data corresponding to the fermions in the ground state at the BCS limit, we can determine µ B 0 in such a way that f P (X) matches the asymptotic theory at the BCS limit. In this case, µ B 0 toward the unitarity limit can be determined by iterative fitting from the BCS limit with the same principle as demonstrated in a previous study [33]. However, as shown in Figure 3, it is difficult to judge whether our data reaches the BCS limit, and the data contains some finite temperature effects at x < x * because the fermions are in the normal state there. Consequently, this method cannot be used to determine µ B for our data.
Although the value of µ B 0 for each value of B is unknown, we can still determine the relative difference among the various values of B using the general properties of dimensionless parameters described in Section II B. Equation (15) and the LDA suggest that all of the experimental data acquired at various Feshbach magnetic fields, B, should be related through a common function G(χ) as where . We plot (χ, G) with different choices of µ B 0 for all data in Figure 4. In the case of µ B 0 = 0 for all B, the plots do not overlap each other as the curves (i). However, we can overlap the plots model-independently by simply adjusting each value of µ B 0 from the BCS limit (B = 1, 050 G) to the unitarity limit, as in the curves (ii). This leaves only a single ambiguous parameter, namely, the arbitrary choice of µ B=1050G 0 . As a result, the experimentally constructed value G(χ) deviates from the true value G true (χ) by a constant offset as G(χ) = G true (χ) + ∆G, where ∆G is independent of χ.
We can derive a lower and an upper bound on the offset value in the following way. The contact density should be positive in the BCS-BEC crossover (C > 0), and C is related to f P (X) by C 15π f ′ P (X), which can be derived from Eq. (3). Then, f P (X) must have positive gradient f ′ P (X) > 0 at an arbitrary interaction parameter X. This restriction gives the lower bound on G(χ) through Eq. (16) as The border of this condition is plotted by the dash-dot red curve in Figure 4a. This condition should be satisfied for the entire experimental curve G(χ). The condition is most restrictive at χ max = −1.0×10 −4 , which is the maximum value achieved in this experiment. Thus, the lower bound is given by G(χ max ) > 0.13, which gives a lower bound on the experimental range of X from Eq. (16) as X min > −2.77.
In addition, the experimental data for f P must satisfy the condition of f P > f P (X min ) from the requirement of f ′ P (X) > 0. This condition gives an upper bound on G(χ) as for χ < χ max . In order to use the above inequality, we still need a lower bound on f P (X min ). Here, we invoke the theoretical asymptotic behavior to derive a bound that is better than the obvious bound of f P (X min ) > 1. We confirmed that calculations of f Asym P (X) up to the first and second orders of x −1 in Eq. (23) give almost the same values at X = −5, which is f P (X = −5) = 1.1. By assuming that this value is reliable, we have the bound f P (X min ) > 1.1, because X min > −5. Using this bound, the condition of Eq. (28) is represented by the dashed red curve in Figure 4a. We confirmed that the entire experimental curve G(χ) lies below the red curve if and only if G(χ max ) < 0.14.
Having established that 0.13 < G(χ max ) < 0.14, we chose the offset such that G(χ max ) = 0.135, which is shown as curve (iii) in Figure 4a. This leads to a systematic error of |∆G| < 0.05. Note that this constant uncertainty in the offset value becomes negligible toward the unitarity limit, because the value of G(χ) around the unitarity limit is approximately four orders of magnitude larger than that around the BCS limit, as shown in Figure  4. In terms of parameter X, the relative systematic error |∆G|/G(χ) is bounded by 0.05|X| 2 because G(χ) = |X| −2 from Eq. (16).  We converted G(χ) to f P (X) according to Eqs. (16), and derived the other dimensionless functions from f P (X) using the thermodynamic relations of Eqs. (17 through 20). Figure 5 shows the results. The red regions indicate systematic errors caused by uncertainty ∆G for determining the offset value of G(χ), which depend on X. The blue regions indicate systematic errors caused by uncertainty of absorption cross-section σ abs , which do not depend on X. In the interaction region of X > −0.18 or x > −0.12, the former errors become one order smaller than the latter. Therefore, all of the thermodynamic quantities have been determined by the uncertainty of σ abs as δfP fP = δfn fn = δσ abs σ abs = 4% and δfE fE = δfµ fµ = 2 3 δσ abs σ abs = 3% around the unitarity limit. Additional systematic errors caused by finite temperature effects in the normal state can be estimated to be around 1% from the experimental temperature t exp 0.1 and t 2 dependence of thermodynamic quantities.
We found that the superfluid transition point corresponds to X * = −1.33 at x * = −1.14 as shown in Figure  5(f). The dotted bars at X * and x * indicate the superfluid transition points. Changes in f P , f n , f E , and f µ and a cusp in f C around the points can be seen. Such critical behavior can be explained as appearing the condensation energy at the normal to the superfluid transition point [34]. ( In the present paper, we analyzed data under the assumption of zero temperature. Therefore, reliable data ranges are the superfluid region at X > X * and x > x * . The quantitative evaluation of critical behaviors requires finite-temperature analysis. This is a subject for future research.

IV. DISCUSSIONS
A. Comparison with previous experiments and many-body theories

Isothermal compressibility: fκ(x)
The obtained dimensionless compressibility is shown by the red circles in Figure 6, along with previously obtained experimental data and a theoretical curve. The vertical axis follows a logarithmic scale.
At the unitarity limit, f κ (0) = 2.66(7) obtained in the preset study is consistent with a previously reported experimental value (2.66(3)), which was determined accurately by measuring the thermodynamics for homogeneous unitary Fermi gases [17]. The results for two different experiments are shown around the unitarity limit. One experiment measured density fluctuations [18]. The compressibility was determined based on density fluctuations according to the fluctuation-dissipation theorem. The other experiment measured the speed of sound propagating in the Fermi gas [19]. The speed of sound can be easily converted to isothermal compressibility, as shown in Eq. (21). While the results of these experiments are not the exact values for a homogeneous system because they were measured for a trapped system, the results of the present study show good qualitative agreement. The theoretical curve was calculated using an extended T -matrix approximation (ETMA) developed in a previous study [35,36] without using any fitting parameters. The theory excellently reproduces our experimental results in the superfluid region.

Pressure: fP (X)
The determined dimensionless pressure is shown by the red curve in Figure 7, along with previous experimental data and a theoretical curve. The width of the red curve indicates the overall systematic error. The vertical axis follows a logarithmic scale.
Our value of f P (0) = 4.35 (17) at the unitarity limit is consistent with the experimental value of 4.34 (7), which is obtained by f P (0) = ξ −3/2 , with ξ = 0.376(4) determined using spin-balanced unitary Fermi gases in a previous study [17]. Our values from the BCS region to the unitarity limit agree qualitatively with the values determined using spin-imbalanced Fermi gases in a previous study [16], whereas our results are slightly larger. The ETMA [35,36] again reproduces our experimental results.

Internal energy density: fE (x)
The determined dimensionless internal energy density is shown by the red curve in Figure 8, along with previous experimental data and various theoretical curves. The width of the red curve indicates the overall systematic errors.  [19], density fluctuations [18], and thermodynamic measurements at the unitarity limit [17], respectively. The light-blue curve indicates the theoretical calculation of ETMA [35,36]. The vertical dash-dot line indicates the superfluid transition point x * obtained in the present study.
At the unitarity limit, our value of f E (0) = 0.375(10) agrees with the experimental value of 0.376(4) [17]. In a previous experiment [16], an approximated curve of f Padé P (X) was prepared using the Padé approximation.
An approximated curve f Padé E (x) was obtained from f Padé P (X) using Eqs. (17) and (18), as shown by the dotted green curve. We found our data to be inconsistent with this approximated curve f Padé E (x). The results of the ETMA [35,36] and the quantum Monte Carlo (QMC) method [32] agree with our result within the error bars. A fixed-node diffusion Monte Carlo (FNDMC) [37] method and the Nozières and Schmitt-Rink (NSR) theory [38] provide larger values around the unitarity limit. A Luttinger-Ward (LW) approach [39] shows smaller values than our data around the unitarity limit.

Contact density: fC(x)
The obtained contact density is shown by the red curve in Figure 9, along with previous experimental data and various theoretical curves. The width of the red curve indicates the overall systematic error. The vertical axis follows a logarithmic scale. The vertical dash-dot line indicates the superfluid transition point at x = x * .
In the previous experiment [15], contact density was The vertical axis follows a logarithmic scale. The red curve indicates the experimental data of the present study. The width indicates the total systematic error. The green circles and the black square indicate the experimental values obtained using spin-imbalanced Fermi gases [16] and spin-balanced unitary Fermi gas [17], respectively. The light-blue curve indicates the theoretical calculation of the ETMA [35,36]. The vertical dash-dot line indicates the superfluid transition point X * of the present study. measured using fermions in the normal state at t = 0.18 (2) t c by radio frequency spectroscopy (RF line shape) and photon emission spectroscopy (PES). The temperature parameter is approximately three times higher than in our experimental condition. Nevertheless, the results of the present study agreed quantitatively within the error bars, indicating that the short-range correlation is insensitive to both temperature and whether the many-body state is the normal state or the superfluid state.
We compared our results with the results of the FNDMC method [37], the NSR theory [40], the LW approach [41], T -matrix approximation (TMA) [42], and ETMA [35,36]. All of these results agreed quantitatively with our data in the BCS region.

B. Superfluid order parameter
In several strong-coupling theories, the order parameter ∆ and chemical potential µ are calculated for given values of n and a by an equation obtained by combining a gap equation and a number equation. The gap equation is often approximated by the ordinary BCS gap equation [11,35,36,38,43], whereas the number equation is derived from the single-particle Green's function or the thermodynamic potential, which includes quantum fluc- The width indicates the total systematic error. The black square indicates the experimental data obtained using spinbalanced unitary Fermi gas [17], and the dotted green curve is an approximated curve based on pressure measurements of [16]. The light-blue curve indicates the ETMA [35,36]. The open brown squares indicate data obtained by the QMC method [32]. The long dashed green curve was obtained by the FNDMC method [37]. The dash-dot purple curve was obtained by the NSR theory [38]. The dashed black curve was obtained by the LW approach [39].
tuations. In this case, ∆ can be uniquely determined by the chemical potential of interacting fermions at T = 0. If this theoretical approach is adequate, we can calculate ∆ by solving the gap equation with the chemical potential determined in the present study. Note that it is not obvious that the superfluid order parameter should satisfy the gap equation in the unitary regime. Nonetheless, it is interesting to evaluate ∆ from thermodynamic quantities using the ordinary BCS gap equation and to compare it with values obtained by other methods.
The BCS gap equation at T = 0 is given by where the function I 1 (µ/∆) is defined in a previous study [44]. When we define the dimensionless superfluid gap as the gap equation can be expressed in dimensionless form, relating x and f ∆ as The red curve indicates the experimental data of the present study. The width indicates the total systematic errors. The blue triangles and the green inverse triangles indicate experimental values determined by radio frequency spectroscopy (rf line shape) and photon emission spectroscopy (PES) [15]. The light-blue curve is the ETMA [35,36]. The long dashed green curve indicates the results obtained by the FNDMC method [37]. The dash-dot purple curve indicates the results obtained by the NSR theory [40]. The short dashed black curve indicates the LW approach [41]. The brown diamonds indicate the results obtained by TMA [42].
Note that ∆ and µ determined by the gap equation of Eq. (29) have the universal ratio of ∆/µ = 1.16 at the unitarity limit.
We substituted f µ (x) determined in the superfluid range (x > x * ) into Eq. (31) and calculated f ∆ as a function of x. The result is shown by the red curve in Figure 10. This curve essentially represents the relation among ∆, n, and a under the assumption that the gap equation holds true. The width of the red curve indicates the overall systematic error. The obtained values are close to the binding energy of the paired fermions directly observed by a spectroscopic experiment [12] as well as that obtained by the quantum Monte Carlo method [45,46]. Consequently, we found that the binding energy of the paired fermions can be simply estimated by substituting f µ (x) into the gap equation. Note that the above observation does not indicate whether the magnitude of the order parameter obeys the gap equation, because the order parameter can deviate from the binding energy in the unitary regime. Direct measurement of the order parameter, such as the measurement of the Higgs mode, is desired in this regard [47].  [44]. The green triangle and the green inverse triangles indicate the theoretical values obtained by the QMC method [45,46]. The blue circles are experimental values obtained by quasiparticle spectroscopy [12].

V. CONCLUSION
We determined various thermodynamic quantities in their dimensionless forms from the BCS region to the unitarity limit for homogeneous spin 1/2 fermions in the superfluid state at the zero-temperature limit. All of the quantities were determined within systematic errors of 4% around the unitarity limit based on standard thermodynamic relations, the scale invariant property, and the local density approximation. In particular, we determined the f E (x) for internal energy density directly from experimental data without using an approximated model function. The dimensionless function represents a universal property shared by various physical systems, and helps to construct the EOS for pure neutron matter, for example.
The thermodynamic quantities obtained here are valuable for studying the strong-coupling properties of a superfluid Fermi gas in the BCS-BEC crossover region. We evaluated various many-body theories using the obtained thermodynamic quantities. These theories agree qualitatively with our results, but the measurement results for f E (x) allowed us to perform more quantitative comparison among those theories. The ETMA, which is a strong-coupling theory involving fluctuations in the Cooper channel, provided the closest results to those of the present study. The details of the ETMA can be found in previous studies [35,36].
We also found that ∆ obtained by substituting the chemical potential of interacting fermions into the ordi-nary BCS gap equation has a value close to the binding energy of the paired fermions. This is new information regarding the relation between thermodynamic quantities and the binding energy in the BCS region close to the unitarity limit. Future experiments, such as experiments involving the measurement of the Higgs mode of the order parameter, will reveal the magnitude of the order parameter as well as the relation between the order parameter and the binding energy. know neither the experimental temperature t exp nor the function for the critical temperature t c (x). Since the local density has a peak value at the bottom of the trap, the ratio t exp /t c takes the minimum value there. Therefore, CF takes finite values when fermions satisfy t exp /t c < 1 at the bottom, and CF becomes zero at t exp = t c (x * ).
We carried out the measurement of CF as follows. We prepared ultracold fermions at various Feshbach magnetic fields using the same experimental procedure as in the present study. Instead of measuring the in-situ col-umn density of the fermions, we measured the centerof-mass (CM) momentum distribution of the paired fermions [27]. We turned off the ODT and the magnetic field simultaneously within 10 µs. At that time, the paired fermions were converted to molecules because the Fermi system was changed to the BEC regime. These molecules expanded with the original CM momentum distribution of the paired fermions during time-of-flight (TOF). We measured the momentum distribution after a TOF of 11 ms by absorption imaging. Figure 11 shows the CF evaluated by fitting a bimodal function to the measured momentum distribution. In order to estimate the magnetic field B * , at which CF becomes zero, we fit an empirical function of CF = a(B * − B) b · Θ(B * − B) to the data with fitting parameters B * , a, and b. The fitting result yielded B * =951 Gauss. We then evaluated the value of x * by x * = 1 kF (n)a(B * ) with the peak density and the scattering length at the magnetic field. In this way, the superfluid transition point was determined to be x * =−1.14. In the region of x > x * , fermions are in the superfluid state.
Next, we estimated the experimental temperature t exp . At the unitarity limit (832.18 Gauss), we observed a condensate fraction of 0.6, as shown in Figure 11. This suggests that t exp < 0.1 at x = 0, according to a previous experimental result [17]. At x * = −1.14, the critical temperature was calculated to be t c (x * ) = 0.06 [39]. Since t exp = t c (x * ), we estimated t exp ∼ 0.06 at x = −1.14. At x < −1.14, it is difficult to estimate the temperature parameter from our experimental data. According to a theoretical work of Ref. [48], the temperature parameter at the trap center decreases during an adiabatic change from the BEC region to the BCS region along the isentropic curve. Consequently, we estimated the temperature parameter to be 0.06 t exp 0.1 from around the BCS region to the unitarity limit.