Pair correlations in the normal phase of an attractive Fermi gas

In a recent paper [Phys. Rev. A 99, 053617 (2019)], the total number of fermion pairs in a spin-balanced two-component Fermi gas of $^6$Li atoms was experimentally probed in the normal phase above the superfluid critical temperature, in order to investigate the sectors of pseudogap and preformed-pair in the temperature-coupling phase diagram. Here, we present a theoretical account of these experimental results in terms of an ab-initio self-consistent $t$-matrix calculation, which emphasizes the role of the pair-correlation function between opposite-spin fermions at equilibrium. Good agreement is found between the available experimental data and the theoretical results obtained with no adjustable parameter.


I. INTRODUCTION
Preformed pairs are meant to be bound states which form above the critical temperature of a fermionic superfluid [1,2]. They are usually associated with the occurrence of a pseudo-gap which can be viewed as a carry-over of the pairing gap in the superfluid phase to the normal phase [3]. Although in the limit of low density and strong fermionic attraction, a preformed pair can be approximately described by a bound state of two fermions of opposite spin, in general it has intrinsically a many-body nature. In order to take into account the many-body character of a pair, it is convenient to describe the pair problem in terms of correlations between the fermions. These correlations are a non-trivial function of temperature, particle density, and the inter-particle coupling.
Preformed pairs were recently studied in an experiment with a spin-balanced two-component Fermi gas of 6 Li in the normal phase [4], where the number of fermion pairs N p was determined by converting all atom pairs to tightly-bound diatomic molecules which afterwards were detected. The pairing fraction N p /N σ (where N σ is the number of all atoms per spin-state) was reported for various temperatures and couplings on the BEC side of the BCS-BEC crossover.
A preliminary theoretical account of the pairing fractions was already presented in Ref. [4], which was obtained by a statistical model of non-interacting atoms and molecules at equilibrium [5,6] as well as by an abinitio diagrammatic t-matrix approach [7]. However, the comparison between experiment and theory presented in Ref. [4] called for further improvements, because the statistical model could not be confidently extended to the crossover region and the t-matrix calculation was lacking refinements which turned out to be important for the crossover region.
Here, we present an improved account of the theoretical approach. We investigate correlations between spinup and spin-down fermions at thermal equilibrium. On the basis of this, we derive a meaningful definition and measure for preformed pairs. We calculate thermodynamic quantities such as the pairing fraction, rather than dynamical quantities such as the pseudo-gap. Nevertheless, the pseudo-gap physics is well contained in our approach. As a consequence, the results of our quantum many-body approach in the crossover region differ significantly from the ones of the statistical atom-molecule model where the fermionic character of the pairs is neglected. In general, we find good agreement between theory and experiment, giving us confidence on the validity of our approach.
Our detailed theoretical interpretation of the experimental data of Ref. [4] and new insights on the separation between the molecular and pseudo-gap regimes are the main results of this paper. In addition, we calculate for a homogeneous Fermi gas (i) the pair correlation function, (ii) Tan's contact (a quantity that sets the overall scale of the pair correlation function), and (iii) the pairing fraction. These three quantities are calculated for different temperatures and couplings across the BCS-BEC crossover. For a trapped system, we also report density profiles and compare them to experimental measurements, and we provide the superfluid critical temperature across the BCS-BEC crossover.
It should be mentioned that the temperature dependence of the contact in the homogeneous case and of the density profiles in the trapped case were already reported in Refs. [8] and [9] within the same self-consistent t-matrix approach of our work, albeit only for the unitary case. We have verified that for this case our results fully agree with the published ones.
The paper is organized as follows. Section II describes the theoretical approach. Section III presents calculated pair fractions N p /N σ for the homogeneous system. Section IV compares these results to the experimental data of Ref. [4] after suitable averaging for the trap. Sec-tion V presents our conclusions. Appendix A discusses the use of conserving approximations for the many-body structure of the pair fraction. Appendix B highlights the circumstances under which the many-body approach to the pair fraction reduces to that of the statistical model. Finally, Appendix C obtains the critical temperature of a trapped low-density Bose gas. Throughout the paper, we set = 1.

II. THEORETICAL APPROACH
The theoretical approach that we set up to account for the experimental results of Ref. [4] on the pair fraction builds on the following ingredients: (i) The definition of the many-body propagator for composite bosons introduced in Appendix A of Ref. [10]; (ii) The formalism developed in Ref. [11] to calculate the pair correlation function of opposite-spin fermions also in the normal phase; (iii) The experience recently nurtured in Ref. [7] on the fully self-consistent solution of the t-matrix approach to a Fermi gas with an attractive inter-particle interaction.
This Fermi gas is made to span the BCS-BEC crossover by varying the (dimensionless) coupling parameter (k F a F ) −1 , where k F = (3π 2 n) 1/3 is the Fermi wave vector associated with the number density n and a F is the scattering length of the two-fermion problem [12]. In practice, the crossover between the BCS and BEC regimes is exhausted within the range −1 (k F a F ) −1 +1 about unitarity where (k F a F ) −1 = 0. In the following, we shall mostly be interested in the coupling region 0 (k F a F ) −1 +1.5 on the BEC side of unitarity for which the experimental data of Ref. [4] are available.
A. Outline of the theoretical expressions to be related with the experimental data Strictly speaking, a pair of spin-up and spin-down fermions can be regarded as a purely bosonic entity only in the BEC regime and at sufficiently low temperatures. In all other cases, one should search for correlations between fermions and define the occurrence of pairing accordingly. Adopting this point of view, which applies also to the so-called Cooper pairs in the BCS regime, is definitely required on the BEC side of unitarity in the normal phase, where the experimental data reported in Ref. [4] were collected. To this end, a suitable definition is needed of what would loosely speaking be referred to as a "preformed pair" in the normal phase of a fermionic superfluid. This definition should be based on a quantum many-body approach where fermions are the elementary constituents of the theory, with no a priori reference to the preformed pairs themselves.
We begin by introducing the bosonic propagator where x = (r, τ ) groups the spatial position r and imaginary time τ , Ψ B (r) is a bosonic field operator, T τ the time-ordered operator, and · · · a thermal average taken at temperature T [13]. In terms of this propagator, the total number of bosons is given by where q is a wave vector, Ω ν = 2πν/β (ν integer) a bosonic Matsubara frequency with β = (k B T ) −1 and k B the Boltzmann constant, and η = 0 + . In the last line of Eq. (2) a homogeneous system has been assumed, for which one may simply write N p = V n p where V is the volume occupied by the system and n p the boson density.
To the extent that the bosonic entities we are considering are made up of fermion pairs, the bosonic operator Ψ B (r) has to be related to its fermionic counterparts ψ σ (r), where σ = (↑, ↓) is the spin projection. This can be achieved by setting where φ(ρ) is a suitable function that should itself embody the correlations within a fermion pair we are after.
On physical grounds, at sufficiently low temperature in the BEC regime it is reasonable to take φ(ρ) as the (normalized) bound-state wave function of the fermionic two-body problem in vacuum, namely, where ρ = |ρ|, whose Fourier transform reads As already mentioned, the definition (3) together with the expression (4) was originally used in Ref. [10] to describe condensed composite bosons well below the superfluid transition temperature T c with fermions treated within the mean-field approximation [14]. The same combination of the expressions (3) and (4) was then utilized in Ref. [4], aiming to account for the quantity N p of Eq. (2) on the BEC side of unitarity in the normal phase above T c , even up to a few times the Fermi temperature T F . In addition, in this case fermions were treated within the self-consistent t-matrix approach [7], with a further trap averaging to comply with the experimental procedure of Ref. [4].
To account for the experimental data of Ref. [4] in a comprehensive way, however, the function φ(ρ) with which the projection is performed in Eq. (3) should acquire a more general form than the expression (4), which is expected to be valid only in the BEC regime at low temperature. Accordingly, in what follows (cf. Section III-A) we will replace the expression (4) by a more general form obtained from the pair correlation function studied in Ref. [11], a form which can thus be utilized even past unitarity towards the BCS regime and up to a temperature of even several times T F .
In addition, we shall see below (cf. Section II-B) that in the diagrammatic expansion of the expressions (1) and (3) one should also retain an "unbound" term that was disregarded in the analysis of Ref. [4] since it is negligible in the BEC limit.
It turns out (cf. Section IV-B) that both these refinements (namely, the inclusion of the above unbound term and the improvement of the expression (4) in terms of the pair correlation function) improve the comparison with the experimental data of Ref. [4], especially just on the BEC side of unitarity. This comparison will also make it possible to distinguish between the pseudo-gap and the molecular regimes mentioned in the Introduction. Specifically, we argue that the molecular regime should be reached when the unbound term contributes in a negligible way to the quantity N p of Eq. (2).

B. Diagrammatic approach to the pair fraction
We pass now to describe the diagrammatic approach that we have adopted for the calculation of the expressions (1)-(3). Although we are interested in the normal phase above T c which the experimental data of Ref. [4] are restricted to, we find it convenient to adopt the Nambu representation of the fermionic field operators in terms of which the diagrammatic approach for the superfluid phase below T c is usually formulated [15]. This is mainly because the concept of fermion pairing originates from the superfluid phase [14], from which it can be extrapolated to the normal phase in the context of the BCS-BEC crossover [12] under suitable circumstances, like in the present case. In addition, through the Nambu representation (6) one finds it easier to deal with the issue of conserving approximations for a fermionic superfluid [16]. This proves important when one selects the set of diagrams that would describe at best the physical problem of interest, with the condition that their numerical implementation remains affordable. We shall discuss this issue in Appendix A.
In terms of the Nambu representation (6), one writes for the fermionic single-particle Green's function and for the fermionic two-particle Green's function with the short-hand notation 1 = (r 1 , τ 1 , ℓ 1 ) and so on, where the Nambu index ℓ = (1, 2) refers to the upper or lower component in the expression (6). Here, G 2 is related to the two-particle correlation function which satisfies the Bethe-Salpeter equation [10,16,17] where is an effective two-particle interaction with Σ the fermionic self-energy. Equation (10) can be formally solved in terms of the many-particle T-matrix, defined as the solution to the equation [10,16,17] by writing The above equations hold quite generally, regardless of the specific approximation for the kernel Ξ defined in Eq. (11). In particular, to the BCS approximation Σ BCS for the self-energy there corresponds the kernel: where only the off-diagonal terms of the BCS self-energy have been retained following a common practice. In the expression (14), τ 3 is the third Pauli matrix [15], x 1 = (r 1 , τ 1 ) and so on, and v(x is the attractive fermionic interaction. For the ultra-cold Fermi atoms of interest, one takes v(r 1 − r 1 ′ ) = v 0 δ(r 1 − r 1 ′ ) of the contact form, where the (negative) strength v 0 is further eliminated in favor of the scattering length a F through a standard regularization procedure [12].
In addition, it will be shown in Appendix A that, due to the specific identification of the Nambu indices relevant to Eq. (15), the many-particle T-matrix of Eq. (12) which solves the Bethe-Salpeter equation for L can be built only in terms of the effective two-particle interaction Ξ of the form (14) [18]. This leaves us with the freedom of endowing the fermionic single-particle Green's function G of Eq. (7) with a suitable additional selfenergy Σ to be selected on physical grounds, without being forced to introduce at the same time related additional terms in the kernel Ξ via Eq. (11).
With these considerations in mind, we have selected this additional self-energy of the form of the fully selfconsistent t-matrix approach, whose performance in the normal phase above T c has been recently tested against those of the non-self-consistent as well as of other partially self-consistent t-matrix approaches [7], with the result that the fully self-consistent one performs best at least as far as thermodynamic quantities are concerned. To the extent that the quantity N p given by the expression (2) of interest here is itself a thermodynamic quantity (consistently with the fact that no analytic continuation from Matsubara to real frequencies is required to calculate it), this choice for Σ within the fully selfconsistent t-matrix approach appears to be adequate for our purposes. In addition, the BCS self-energy Σ BCS , which has served to obtain the kernel Ξ BCS of Eq. (14), vanishes identically in the normal phase and no longer needs to be considered in what follows. For a homogeneous system, we can further make use of the Fourier representation and rewrite Eq.(15) as: where ω n = (2n + 1)π/β (n integer) is a fermionic Matsubara frequency (the conventions for the Nambu indices are specified in Fig. 10 of Appendix A). The expression (16) will be utilized in Eq. (2) to obtain the number of pairs N p . Solving then for the many-particle T-matrix of Eq. (12) as described above and entering the result in Eq. (13) for L, Eq. (16) reduces eventually to the form: Here, are "form factors" associated with the particle-particle bubble where j = (1, 2), and is the particle-particle propagator in the normal phase where is the regularized particle-particle bubble [12]. We emphasize again that the fermionic single-particle Green's functions G entering the expressions (18) and (20) are meant to be obtained within the self-consistent t-matrix approach in the normal phase [7]. What is still left to be specified is the form of the wave function φ(p) that enters Eq. (18). We have already mentioned that, in the theoretical diagrammatic approach to N p presented in Ref. [4], φ(p) was taken of the form (5) corresponding to the fermionic two-body problem. With this choice, however, meaningful results could be obtained only towards the BEC edge of the BEC side of the unitary region. To overcome this limitation, here we adopt a more general form for φ(p) which will be obtained from the pair correlation function, as discussed in Section III-A below.
In addition, in Ref. [4] the first term on the right-hand side of Eq. (17) was not retained. As anticipated in Section II-A, this term will be referred to as the "unbound" term as opposed to the "bound" term discussed below. Here, we are going to keep this "unbound" term and show that it gives a non-negligible contribution to N p , through the pairing correlations contained both in the fermionic single-particle Green's function G and in the wave function φ(p) that enter the expression (18) with j = 2. Accordingly, through this term spin-↑ and spin-↓ fermions correlate with each other indirectly via their separate interaction with the environment.
In contrast, the second term on the right-hand side of Eq. (17) is referred to as the "bound" term, because in this case spin-↑ and spin-↓ fermions correlate with each other directly through their inter-particle attractive interaction. The result for N p obtained from this term will be shown to reduce to that of the statistical model of atom-molecule equilibrium introduced in Refs. [5,6], past the BEC side of the unitary region and for not too high temperatures above T c . The reasons for the success of the statistical atom-molecule model in this sector of the phase diagram will be discussed in Appendix B.

C. Single-particle Green's function
As discussed in Section II-B, the single-particle Green's function G(p, ω n ) that enters the expressions (18) and (20) is taken within the fully self-consistent t-matrix approach. It then reads: where G 0 (p, ω n ) = [iω n − ξ(p)] −1 is the non-interacting counterpart with ξ(p) = p 2 /(2m) − µ (m being the fermion mass and µ the chemical potential) and is the self-energy with Γ(q, Ω ν ) given by Eqs. (19) and (20). The chemical potential is eventually obtained from the fermionic density n σ via the relation where n ↑ = n ↓ = n/2 like in Ref. [4]. The numerical calculation of the expressions (21)-(23) will be implemented by taking advantage of the detailed procedures recently reported in Ref. [7]. In addition, the strong-coupling (BEC) limit of the expressions (21)-(23), together with that of the expressions (2) and (17)- (20), will be examined in Appendix B, to determine under what circumstances the results for n p and n σ obtained by our diagrammatic quantum manybody theory reduce to those of the statistical model of atom-molecule equilibrium developed in Refs. [5,6].

III. RESULTS FOR A HOMOGENEOUS GAS
In this Section, we implement the calculation of the bosonic density n p obtained from Eqs. (2) and (17) for a homogeneous gas, as a function of coupling and temperature. The information gathered in this way will be used in Section IV when dealing with a trapped gas, by performing a trap average within a local-density approach. At that point it will be possible to compare the theoretical results with the experimental data of Ref. [4].
The main ingredients of the calculation of n p are the single-particle Green's function G(p, ω n ) and the wave function φ(p) that enter Eqs. (18)- (20). The calculation of G(p, ω n ) was already considered in Section II-C. It thus remains to consider the calculation of the wave function φ(p), as discussed next.

A. Pair correlation function
Our interpretation of the experimental data of Ref. [4] rests on the occurrence of correlations between spin-up and spin-down fermions at equilibrium. The preliminary theoretical account of those experimental data presented in Ref. [4] took the wave function φ(p) entering Eq. (18) of the form (5) associated with the fermionic two-body problem. This form, however, proves able to account for the correlations between spin-up and spin-down fermions only in the BEC regime of coupling and at low enough temperature. As anticipated in Section II-B, we now consider a more general form for φ(p) which is obtained from the pair correlation function This function contains information about correlations between fermions of opposite spins at a distance ρ = |ρ| apart. This quantity was studied in detail in Ref. [11] throughout the BCS-BEC crossover, both in the superfluid phase below T c and in the normal phase above T c .
Here, we consider the formalism of Ref. [11] above T c and rephrase it in terms of the fully self-consistent t-matrix approach that was summarized in Section II-C.
Within the fully self-consistent t-matrix approach, the expression (24) for g ↑↓ (ρ) can be cast in the form [11]: (26) Here, the fully self-consistent G's are considered, while in the original Ref. [11] non-interacting G 0 corresponding to the non-self-consistent approximation were utilized.
It was also shown in Ref. [11] that g ↑↓ (ρ) given by the expression (25) recovers the short-range behavior related to Tan's contact C [20][21][22] such that lim ρ→0 (4π) 2 C ρ 2 g ↑↓ (ρ) = 1 irrespective of coupling and temperature. We have reproduced here these analytic results within our fully self-consistent t-matrix approach, with the numerical values of C obtained in agreement with Ref. [7].
Examples of the spatial profiles of the pair correlation function g ↑↓ (ρ) are shown in Fig. 1, for several couplings and temperatures above T c . Reported in each inset are also the respective values of the contact C, from which the numerical values of g ↑↓ (ρ) can be explicitly reconstructed. Note the oscillatory behavior of g ↑↓ (ρ), which is present on the BCS side at low temperatures but quickly fades away either by moving towards the BEC side or by increasing temperature. Due to this oscillatory behavior, g ↑↓ (ρ) may acquire negative values which correspond to a weaker correlation with respect to the uncorrelated value n ↑ n ↓ = (n/2) 2 [11]. This behavior, however, will not affect our argument below, whereby the oscillations about zero (whenever present) will be averaged out.
It can be further verified from the expression (25) that, in the BEC limit and at sufficiently low temperatures, g ↑↓ (ρ) reduces to the product of the density n σ = n/2 of a single fermionic species times the square of the wave function (4) corresponding to the fermionic two-body problem. This suggests that the function φ(p), to be utilized in the form factors (18), can be extracted from the pair correlation function g ↑↓ (ρ) also away from the BEC limit and at high temperatures. To this end, we adopt the following strategy.
We begin by fitting the spatial profiles of the function where b is a parameter that depends on coupling and temperature (note that the function (28), too, has unit value at ρ = 0). We then take the square root of the expression (28) to extract φ(ρ), and multiply the result by a suitable normalization factor N , thus writing: with where erfc(z) is the complementary error function of (complex) argument z [23]. Note that the two-body wave function (4) is recovered for b → 0. Finally, we take the Fourier transform of the expression (29) and obtain the desired result: where p = |p|. This expression recovers Eq. (5) in the limit b → 0. Figure 2 shows the behavior of the parameter b obtained in this way, over a wide range of coupling and temperature relevant to the experiment of Ref. [4]. In particular, for sufficiently high temperature and irrespective of coupling, b is expected to become proportional to mkB T is the thermal wavelength. To evidence this linear behavior of b vs T at high temperature, the inset of Fig. 2 plots the derivative of b with respect to T for the same temperature range and couplings of the main panel. In all cases, we have found that, at high temperature, this derivative is well reproduced by the expression kB 2m ∂b ∂T = 0.25 − 0.175(k F a F ) −1 T F /T . The fitting function φ(ρ) given by Eq. (29) focuses on the short-range part of the pair-correlation function g ↑↓ (ρ) given by Eq. (24), which is dominated by the intrapair correlations of relevance here. It thus disregards a possible long-range part of g ↑↓ (ρ) which may include correlations between spin-↑ and spin-↓ fermions belonging to different pairs (although this long-range part does not occur within the t-matrix approach adopted here).

B. Pair fraction
We are now in a position to calculate the pair density n p given by together with the fermionic density n σ given by Eq. (23), for a homogeneous system as a function of coupling and temperature.
To begin with, Fig. 3 compares the pair fraction n p /n σ at T c over a wide range of the coupling (k F a F ) −1 , as obtained by the fully self-consistent and non-self-consistent t-matrix approaches. As for other thermodynamic quantities [7], also in this case the fully self-consistent approach proves superior to the non-self-consistent one, to the extent that the ratio n p /n σ should never exceed unity. Accordingly, from now on results obtained by the fully self-consistent approach will only be presented. In addition, the use of the two-body form (5) for φ(p) in the form factors (18) is seen to lead to unstable results upon entering the unitary regime with (k F a F ) −1 +1. Abandoning the two-body form (5) in favor of the expression (31) associated with the pair correlation function is thus expected to yield a definite improvement over the theoretical analysis made in Ref. [4] when accounting for the experimental values of the pair fraction for the trapped system (cf. Section IV-B below). In Fig. 4 the pair fraction n p /n σ is shown over a wide range of temperature and a selected number of couplings across unitarity. In particular, this figure compares the results obtained by including (full lines) or neglecting (dashed lines) the "unbound" term represented by the term -F 2 on the right-hand side of Eq. (17). One sees that inclusion of this unbound term over and above the bound term (represented by the second term on the righthand side of Eq. (17)) leads to substantial differences, especially in the unitary regime at low temperature. The unbound term was not included in the diagrammatic approach to the pair fraction presented in Ref. [4]. It will be shown in Section IV-B that the agreement with experimental data will be definitively improved by its inclusion.
In preparation for this comparison, Fig. 5 shows three contour plots where a given value of the pair fraction n p /n σ is seen to evolve in the T -vs-(k F a F ) −1 phase diagram. Similarly to what was done in Fig. 4, for each of the three values of n p /n σ here reported the numerical results have been obtained by including (full lines) or neglecting (dashed lines) the unbound term in Eq. (17). In all cases, the difference between these two sets of results turns out to be substantial as soon as entering the unitary regime with (k F a F ) −1 +1. This implies that, in this regime of most physical interest, the fermionic character of the constituent particles reveals itself. As a consequence, this counting has to rely on filtering the occurrence of fermionic correlations, and not merely on signaling the presence of bound pairs which would instead apply to the molecular regime with (k F a F ) −1 +1.
To confirm this point of view, Fig. 5 also shows for comparison the contour plots of n p /n σ corresponding to the statistical model (dotted lines), as obtained from the law of mass action where ε 0 = (ma 2 F ) −1 is the two-body binding energy, which results from the integrals in Eq. (B11) of Appendix B by neglecting ±1 in the denominators therein. It turns out that the results of the statistical model coincides with those of the quantum many-body approach that includes only the bound term, but only at most up to (k F a F ) −1 ≈ 0.6 after which the molecular regime with the two-body wave function (4) loses its meaning.

IV. RESULTS FOR A TRAPPED GAS AND COMPARISON WITH EXPERIMENTAL DATA
The results obtained in Section III for n p given by Eq. (32) and for n σ given by Eq. (23) refer to a homogeneous system. In order to compare with the experimental data of Ref. [4], these theoretical results need to be averaged over the trap that contains the Fermi gas.

A. Trap average
When considering a Fermi gas trapped in an anisotropic harmonic potential of the type one can adopt a local-density approach and obtain the total number N p of pairs and the total number N σ of fermions in the trap in the following way. One first replaces the fermionic chemical potential µ entering the single-particle Green's function G(p, ω n ) of Eq. (21) by µ → µ − V (r), thereby obtaining the local function G(p, ω n ; r). One then replaces G(p, ω n ) → G(p, ω n ; r) everywhere this function occurs, namely, in the expressions (17) T F ) is the value of n(r = 0) for the non-interacting gas at T = 0 within a local-density approximation.
for fermions. Finally, one integrates the expressions of the local densities n p (r) and n σ (r) obtained in this way over the spatial variable r, to get the total number of pairs N p and the total number of fermions N σ with spin σ. The value of the fermionic chemical potential µ for the trap is eventually determined for given coupling and temperature by solving for µ as a function of N σ . In practice, in the experiment of Ref. [4] typical values of ω x = ω y range from 2π × 300 Hz to 2π × 1.6 kHz, while ω z = λω x = 2π × 21 Hz (with λ < 1).
Profiles of the total fermionic isotropic density n ′ (r ′ ) obtained in this way are shown in Fig. 6, for several couplings across unitarity and temperatures in the normal phase. The coupling parameter (k F a F ) −1 associated with the trap is expressed in terms of k F = √ 2mE F , where E F = ω 0 (3N ) 1/3 is the Fermi energy of the trap and N = N ↑ + N ↓ is the total number of fermions. (In the experiment of Ref. [4], typical values of N range from 3 × 10 4 to 3 × 10 5 .) The values of the critical temperature T c for the trap case, reported in Fig. 6 only for three specific couplings, can be obtained throughout the whole BCS-BEC crossover. This information is important also to verify whether the experimental values of the pair fraction in the trap of Ref. [4] were measured in the normal phase. To calculate T c for the trap, we adopt again a localdensity approach and define a local Fermi temperature T F (r) such that k B T F (r) = 3π 2 n(r) 2/3 /(2m). This implies that the local Fermi temperature, like the density n(r), has its maximum value at r = 0, to which there corresponds a minimum value of T /T F (r) for given temperature T . Accordingly, the central portion of the cloud density is where superfluidity is first established upon lowering the temperature from the normal phase.
To obtain T c for the trapped system, we then apply the Thouless criterion in terms of the particle-particle propagator (19) in the normal phase, where now µ(r = 0) = µ − V (r = 0) = µ is the fermionic chemical potential for the trap calculated at the critical temperature T c . Details on how the variables (T c , µ c ) have been determined by solving the Thouless criterion in conjunction with the density equation are given in Appendix B of Ref. [7]. Figure 7 shows the results of our calculation for T c in the trap across the BCS-BEC crossover. The results of the fully self-consistent t-matrix approach (full line) are also compared with those of its non-self-consistent counterpart (dashed line). While the two calculations essentially coincide with each other in the BCS regime (k F a F ) −1 −1, they differ considerably on the BEC side of unitarity. We attribute this difference to the occurrence of a residual interaction between composite bosons in the BEC regime (k F a F ) −1 +1, which is present within the fully self-consistent but absent within the nonself-consistent calculation [7].
To make a check on the results of our numerical calculation, also shown in Fig. 7 are the results for T c (dasheddotted line) obtained for a low-density trapped Bose gas with a residual interaction specified by the scattering length a B (cf. Appendix C), where for internal consistency the (approximate) value a B = 1.16a F that results from the fully self-consistent t-matrix approach [7] was considered. In this way, we can confirm quantitatively the effects of a B on T c for the trapped system in the BEC regime, which are contained in the fully self-consistent tmatrix approach. For comparison, the inset reports additional bosonic calculations for: (i) a B = 2.0a F which corresponds to the residual bosonic interaction being treated at the level of the fermionic exchange diagrams [24]; (ii) a B = 0.75a F when the T -matrix for the dimer-dimer scattering built on these exchange diagrams is further considered [24]; (iii) The exact value a B = 0.6a F obtained either by a numerical solution of the four-body Schrödinger equation [25] or by a full diagrammatic treatment in the zero-density limit [26].
Finally, it should be mentioned that the value T c /T F = 0.2074, which we have obtained at unitarity by the fully self-consistent calculation, coincides with that obtained in Ref. [9] by the same approach. However, our calculation for T c is extended to the whole BCS-BEC crossover while that of Ref. [9] was limited to unitarity only.

B. Comparison between theory and experiment
A first quantity to be compared with the experimental data of Ref. [4] is the so-called axial density n a (z) where z runs along the main axis of the trap, which is obtained by integrating the full density n(x, y, z) over the radial directions x and y. Specifically, the experimental profiles n a (z) can be compared with their theoretical counterparts n ′ a (z ′ ), obtained by integrating over x ′ and y ′ the isotropic profiles n ′ (r ′ ) = n ′ (x ′ , y ′ , z ′ ) (like those shown in Fig. 6) and then performing the rescaling n a (z) = λ 2/3 n ′ a (λ 2/3 z) .
(37) Figure 8 shows this comparison for three sets of values of temperature, coupling, and anisotropy λ. In all cases, excellent agreement results between the experiment and the quantum many-body approach with no adjustable parameter. The figure shows also the comparison with the statistical atom-molecule model, for which notable deviations from the experiment occur, as expected, for low temperature and close to unitarity. Finally, Fig. 9 presents the comparison of the pairing fraction N p /N σ obtained by our ab initio quantum manybody calculation with the experimental data of Ref. [4] over the temperature-coupling phase diagram (where k F and T F now refer to the trapped system). The comparison is made for three characteristic values of N p /N σ . In all cases, good agreement is obtained between theory and experiment (we emphasize that the theoretical results have been obtained with no adjustable parameter).
In particular, this comparison shows that the contribution of the unbound term significantly improves the agreement of our calculations with the experimental data, despite the presence of the trap which acts to suppress the contribution of the unbound term (which is evident by comparing Figs. 5 and 9). This suggests that the experimental data probe indeed the pairing correlations between spin-up and spin-down fermions as defined by our formalism.
From this comparison one can argue that the crossover, between the pseudo-gap regime (where the fermionic character of the constituent particles matters) and the molecular regime (where only the presence of bosonic pairs is relevant), sets in about where the theoretical results for N p /N σ , obtained with and without the unbound term, start departing from each other. This argument cannot be made in terms of the statistical atom-molecule model [4], that misses the contribution of the unbound term.

V. CONCLUDING REMARKS
In this paper, we have provided a detailed account of a theoretical approach to interpret the experimental data reported in Ref. [4] in a quantitative way. By this approach, from the data Ref. [4] we have been able to unravel how the occurrence of pairing correlations between spin-up and spin-down fermions at equilibrium develops, as a function of temperature in the normal phase and of coupling on the BEC side of unitarity. What we claim to have learned from this is how the pseudo-gap regime (where fermions matter) and the molecular regime (where only composite bosons matter) separate from each other. This should be considered rather remarkable, since this result was extracted from experiment [4] where an equilibrium quantity was measured (i.e. the number of fermion pairs) and not a dynamical quantity (the excitation gap).
From the theoretical side, to account for the experi-mental data we have taken advantage of several favorable circumstances. On the one hand, since the number of fermion pairs in a Fermi gas undergoing the BCS-BEC crossover is an equilibrium quantity, it can be accounted for quite well in terms of the fully self-consistent t-matrix approach [7]. On the other hand, this physical quantity that was measured experimentally by its own nature does not require one to endow the theory with a series of complicated Aslamazov-Larkin and Maki-Thompson diagrams, which should otherwise be included to fulfill conservation criteria when addressing dynamical response functions [16], to the extent that the single-particle selfenergy is treated within the fully self-consistent t-matrix approach. In addition, our emphasis here on fermionic correlations has drawn on our previous experience on the pair-correlation function in the normal phase, which was addressed in detail in Ref. [11] within the non-selfconsistent t-matrix approach and here extended to the fully self-consistent one.
Along these lines, future perspectives, that could reinforce our argument about the evidence for the separation between the (fermionic) pseudo-gap and the (bosonic) molecular regimes, may hinge on the possibility of extending the measurements of the ratio N p /N σ towards unitarity at temperatures close enough to T c .
In addition, to highlight experimentally the relevance of the correlations induced indirectly by the environment between spin-↑ and spin-↓ fermions, which are embodied in the "unbound" term in the expression (17), it could be worth to consider repeating the experiment of Ref. [4] by replacing the harmonic trap with a box trap along the lines of Ref. [27]. In this way, one should be able to amplify the difference between the values of the pair fraction obtained with and without the inclusion of the unbound term, as one may anticipate by comparing the results of Fig. 5 for the homogeneous case with those of Fig. 9 for the trapped case.
It is, finally, interesting to draw a physical connection between our finding about the indirect correlations established between spin-↑ and spin-↓ fermions through their environment and the recent results of Ref. [28] about the way the quark-gluon structure of a nucleon bound in an atomic nucleus is modified by the surrounding nucleons. In both cases, it is the environment that plays an important role in modifying the properties of what would be a bound system in isolation.  In Section II-B we have argued that only the form (14) of the effective two-particle interaction Ξ is of relevance for the calculation of the bosonic propagator G B (q, Ω ν ) of Eq. (16) (and thus of the quantity N p of experimental interest). We have also anticipated that the reason for this is to be found in the specific sequence of Nambu indices appearing in the expression (15) from which Eq. (16) is derived. Here, we show specifically how the diagrammatic contributions to Ξ, that would derive from the tmatrix approach for the fermionic self-energy Σ, cannot modify this result. Under different circumstances, like for the calculation of the density and spin response functions, on the other hand, the diagrams for Ξ corresponding to the Aslamazov-Larkin (AL) and Maki-Thomson (MT) contributions would instead result from the t-matrix approach for Σ (see, e.g., Fig. 3 of Ref. [29]). In our case, the importance of introducing the t-matrix approach for Σ arises from the need of obtaining an accurate description of the thermodynamic properties of the Fermi gas in the normal phase [7].
Probably the simplest way to convince oneself that the AL-type and MT-type contributions to Ξ, which would result from the t-matrix self-energy taken below T c , do not contribute to the expression (15) of the pair propagator G B once carried over to the normal phase above T c , is to draw these contributions in a diagrammatic way. This is done in Fig. 10. Here, the series of ladder diagrams that approximate the many-particle T-matrix in the brokensymmetry phase is reported in panel (a), while the corresponding t-matrix self-energy is shown in panel (b). For simplicity, in these diagrams only the Nambu indices have been explicitly indicated, while the space and imaginary time variables are not reported since they are not essential to the following argument. The crucial point is that for the T-matrix of panel (b) only combinations with Nambu indices ℓ L = ℓ ′ L and ℓ R = ℓ ′ R occur, owing to the inter-particle interaction of the contact form that we have adopted (cf. also Ref. [10]). In addition, only combinations with ℓ L = ℓ R and ℓ ′ L = ℓ ′ R will survive when these diagrams are extrapolated to the normal phase. As a consequence, a typical example of MT contribution is shown in Fig. 10(c), while a typical example of AL contribution is shown in Fig. 10(d). In all cases, it turns out that at least two single-particle Green's functions with off-diagonal Nambu indices would be required to match these contributions to Ξ with the Nambu indices appearing in the expression (15). Since the off-diagonal (anomalous) single-particle Green's functions vanish in the normal phase above T c , the MT-and AL-type contributions to Ξ vanish, too, and do not affect the expression (15) which is relevant for the calculation of N p above T c . This proves our statement. It is interesting to determine under what physical circumstances the expressions for the total number of bosons N p and for the total number of spin-σ fermions N σ of our fully quantum many-body approach reduce to those of a statistical model of a fermion-boson mixture at equilibrium [5,6].
To this end, we consider a homogeneous system, for which N p = V n p and N σ = V n σ are expressed in terms of the respective densities. By our quantum many-body approach, n p is given by Eq. (32) with G B given by the expression (17), while n σ is given by the expression (23). To recover the physics of a fermion-boson (or atommolecule) mixture, one requires the fermionic coupling to be sufficiently strong in the BEC regime and the temperature sufficiently low, for the internal structure of the composite bosons (dimers) to become irrelevant.
In this limit, the fermionic chemical potential µ becomes the largest energy scale of the problem and is written in the form µ B = 2µ + ε 0 , where ε 0 = (ma 2 F ) −1 is the dimer binding energy and µ B the dimer chemical potential [12]. The expression (26) then reduces tõ which, together with the expression (5) for φ(p) appropriate to this limit, yields the the following approximate form for the form factors (18) [10]: This implies that, in the BEC limit where a F → 0 + , the "unbound" term F 2 vanishes faster than F 1 and can thus be neglected in the expression (17). In addition, in the same limit the particle-particle propagator Γ(q, Ω ν ) of the "bound" term in the expression (17) acquires the polar form [12]: Combining these results together, one gets eventually for the bosonic density: in terms of the Bose-Einstein distribution of argument ξ B (q) = q 2 4m − µ B . To determine n σ in the BEC limit at sufficiently low temperature, we consider the expression (23) where we expand the single-particle Green's function (21) in series of the self-energy Σ G(p, ω n ) ≃ G 0 (p, ω n )+ G 0 (p, ω n ) Σ(p, ω n ) G 0 (p, ω n ) + · · · (B5) where G 0 (p, ω n ) = [iω n − ξ(p)] −1 is the non-interacting single-particle Green's function, by again relying on the fact that the fermionic chemical potential µ entering ξ(p) = p 2 /(2m) − µ is the largest energy scale in the problem. We thus obtain: n σ ≃ dp (2π) 3 1 β n e iωnη G 0 (p, ω n ) + dp (2π) 3 1 β n G 0 (p, ω n ) 2 Σ(p, ω n ) + · · · ≡ n (0) σ + n (1) σ .
Here, n (0) σ = dp (2π) 3 1 β n e iωnη G 0 (p, ω n ) = dp (2π) 3 1 e βξ(p) + 1 (B7) coincides with the density n f of fermions (atoms) expressed in terms of the Fermi-Dirac distribution of argument ξ(p), and owing to the approximate form for the self-energy (22) which is valid in this limit. With the polar approximation (B3) for Γ(q, Ω ν ) and the further approximate result (cf., e.g., Section 3.1 of Ref. [12]) dp (2π) 3 n σ = n f + n p = dp (2π) 3 1 e βξ(p) + 1 At this point, the fermionic chemical potential µ can be eliminated from Eq. (B11) by fixing the value of n σ therein, with the bosonic chemical potential µ B = 2µ+ε 0 following in a consistent way. There remains to find an explicit connection with the expressions of the fermion-boson (atom-molecule) model, which were obtained in Refs. [5,6] in the classical limit and used in Ref. [4] to account for the experimental data in the BEC regime of the phase diagram. To this end, we consider the classical limit of the expressions (B11) by neglecting ±1 in the denominators, and perform the trap average by replacing µ → µ− V f (r) and µ B → µ B − V p (r) and integrating over the space variable r, similarly to what was done in Section IV-A. Here, V f/p (r) = 1 2 M f/p ω 2 x x 2 + ω 2 y y 2 + ω 2 z z 2 (B12) is the (anisotropic) harmonic oscillator potential commonly considered for ultra-cold gases, with M f = m for fermions (atoms) and M p = 2m for bosons (molecules).
The results for the total number of fermions N f and the total number of bosons N p then become: N f ≃ dr dp (2π) 3 e −β p 2 2m +V f (r)−µ = k B T ω 0 3 e µ/kB T (B13) and N p ≃ dr dq (2π) 3 e −β q 2 4m +Vp(r)−µB = k B T ω 0 3 e µB /kB T (B14) where ω 0 = (ω x ω y ω z ) 1/3 is the average trap frequency (cf., e.g., Refs. [30,31]). From these results it follows that from which, by replacing ω 0 = E F /(6N σ ) 1/3 where E F is the Fermi energy for the trap, one recovers the expression reported in Appendix A of Ref. [4]. More generally, N p and N f for the trapped case could be obtained in closed form directly from Eqs. (B4) and (B7), in terms of Li 3 (e βµB ) for bosons and Li 3 (−e βµB ) for fermions (where Li n (z) is the poly-logarithmic function of index n and argument z). The expression (B15) generalizes to a harmonically trapped system the law of mass action valid for a homogeneous system [32].
Finally, it is worth summarizing what is lost when passing from the fully quantum many-body approach to its simplified version obtained above. To get this simplified version, in Eq. (17) we have (i) neglected the "unbound" term −F 2 (q, Ω ν ), (ii) approximated φ(p) in the expression (18) for F 1 (q, Ω ν ) by the two-body form (5) and taken µ = −ε 0 /2 therein with ε 0 ≫ k B T , and (iii) approximated Γ(q, Ω ν ) by the polar form (B3); while in Eq. (23) we have performed the expansion (B5) with the typical approximations that apply to the BEC limit at low temperature when µ is the largest energy scale in the problem. None of these approximations, however, is valid either away from the BEC limit when approaching unitarity at any temperature, or in the BEC limit itself for sufficiently high temperature. In both these cases, the fermionic nature of the "preformed pairs" manifests itself and only fermionic correlations remain physically relevant. On physical grounds, the results of the quantum many-body approach and of the statistical fermionboson model differ from each other to the extent that the latter bears essentially on the chemical reaction (dimer ←→ spin-↑ + spin-↓) for molecules that break up into atom pairs and vice-versa, with no regard on the way the molecules are formed by the laws of quantum mechanics and on the effects that the surrounding environment might exert on them through inter-particle collisions.
In this context, it is interesting to explicitly verify to what extent the results of the quantum many-body approach (Q) and of the classical statistical model (C) differ from each other in the BEC limit of the homogeneous k F a F