{'fraction_non_alphanumeric': 0.03713071515648677, 'fraction_numerical': 0.014572011256565381, 'mean_word_length': 4.552729803251724, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thermodynamics of the BCS-BEC crossover

We present a self-consistent theory for the thermodynamics of the BCS-BEC crossover in the normal and superfluid phase which is both conserving and gapless. It is based on the variational many-body formalism developed by Luttinger and Ward and by DeDominicis and Martin. Truncating the exact functional for the entropy to that obtained within a ladder approximation, the resulting self-consistent integral equations for the normal and anomalous Green functions are solved numerically for arbitrary coupling. The critical temperature, the equation of state and the entropy are determined as a function of the dimensionless parameter $1/k_Fa$, which controls the crossover from the BCS-regime of extended pairs to the BEC-regime of tightly bound molecules. The tightly bound pairs turn out to be described by a Popov-type approximation for a dilute, repulsive Bose gas. Even though our approximation does not capture the critical behaviour near the continuous superfluid transition, our results provide a consistent picture for the complete crossover thermodynamics which compare well with recent numerical and field-theoretic approaches at the unitarity point.


I. INTRODUCTION
The problem of a two component attractive Fermi gas near a resonance of the s-wave scattering length describing the effective interaction is one of the basic many-body problems which has been brought into focus by the recent realization of molecular condensates in ultra-cold Fermi gases [1,2,3] and the subsequent exploration of the crossover from a Bose-Einstein-Condensate (BEC) to a BCS-like state of weakly bound fermion pairs [4]. Clear signatures for the existence of paired fermion superfluidity with cold atoms have been provided by spectroscopic measurements of the gap [5] and the observation of a vortex lattice on the BCS-side of the transition [6]. The ability of tuning the interaction in cold Fermi gases through Feshbach resonances relies on the resonant coupling of the scattering state near zero energy of two colliding atoms with a bound state in a closed channel [7]. A particularly challenging problem arises right at the Feshbach resonance, where the two-particle scattering length is infinite [8,9]. Precisely at this point and for broad Feshbach resonances, where the range r ⋆ of the effective interaction is much smaller than the mean interparticle spacing [10,11,12,13], the full many-body problem has the Fermi energy ε F as the only energy scale. As pointed out by Ho [14], the thermodynamics of the unitary Fermi gas is then a function only of the dimensionless temperature θ = T /T F . More generally, as emphasized recently by Nikolic and Sachdev [15], the universality also extends to the behaviour away from the Feshbach resonance, as long as the broad resonance condition k F r ⋆ ≪ 1 is obeyed. Thus, for instance, the critical temperature T c /T F for the transition to superfluidity is a universal function of the inverse coupling constant 1/k F a.
A quantitative theoretical understanding of the manybody problem near a Feshbach resonance has been developed recently through numerical calculations. In particular, at zero temperature and for a homogeneous system, fixed-node Green function Monte Carlo calculations provide quantitative results for the gap parameter [16], the equation of state [17], and also the momentum distribution, the condensate fraction and the pair size [18] of the ground-state for arbitrary values of 1/k F a. As expected in the case of an s-wave resonance [19], these quantities all evolve continously as the coupling is varied from the BCS to the BEC-limit. An important ingredient in these results is their account for the repulsive interaction between strongly bound dimers in the BEC-limit with scattering length a dd ≈ 0.60 a > 0 [20]. This interaction is missing in the early qualitative descriptions of the T = 0 BCS-BEC crossover problem by Eagles [21] and Leggett [22], which are based on using the standard BCSgroundstate as a variational Ansatz for arbitrary coupling [23]. Beyond a purely numerical approach, the BCS-BEC crossover problem has recently become amenable also to analytical methods via an ǫ = 4 − d expansion [24]. It is based on the observation [25] that at the unitarity point in d = 4 (i.e. the point where a two-particle bound state appears) the two-component Fermi gas is in fact an ideal Bose gas, because a zero range interaction in d = 4 can bind a state only at infinitely strong attraction. In two dimensions, in turn, binding appears at arbitrary small couplings and the unitary Fermi gas in d ≤ 2 coincides with a non-interacting one [25]. Within a field theoretic description, the physically interesting 3D problem can thus be approached by extrapolating expansions from the upper and lower critical dimensions d = 4 and d = 2 respectively [26]. At finite temperature, numerical calculations are available for the ther-modynamics at the unitarity point. They are based on an auxiliary field quantum Monte-Carlo method for the continuum problem [27] and on a diagrammatic determinant Monte-Carlo method for the negative U -Hubbard model [28]. Field-theoretic results at finite temperature, which open the possibility for controlled and systematic expansions for the crossover thermodynamics have been obtained very recently by Nishida [26] within an expansion around both the upper or lower critical dimension and by Nikolic and Sachdev [15] within a 1/N expansion for a 2N component Fermi gas.
Our aim in the following is to present a self-consistent many-body theory for the thermodynamics of resonantly interacting fermions at arbitrary temperatures and detuning, which directly addresses the physically relevant case of a three dimensional, two-component Fermi gas. The theory is based on a conserving, so-called Φ-derivable approach to the many-body problem, in which the exact one-or two-particle Green functions serve as an infinite set of variational parameters. It is an extension of earlier work by one of us [29,30,31] and employs a combination of the Luttinger-Ward and DeDominicis-Martin approach for obtaining the grand canonical potential and the entropy, respectively. The condition of gaplessness is enforced by a modified coupling constant, thus accounting for the proper low energy behaviour in terms of a Bogoliubov-Anderson mode. We provide quantitative results for the critical temperature, the equation of state and the entropy near the Feshbach-resonance as a function of both T /T F and 1/k F a. In spite of the fact that the critical behaviour at the continuous superfluid transition is not captured correctly in our approach, which gives rise to a weak first order transition, the results provide a quantitative and consistent picture of the crossover which obeys thermodynamic relations at the percent level. Our variational method is complementary both to purely numerical and to field theoretic approaches to the problem. The results can be used e.g. to predict the final temperature reached after an adiabatic ramp across the Feshbachresonance starting deeply in the BEC-regime [32] or to determine the size of the atom cloud in a harmonic trap near unitarity as a function of temperature.
The paper is organized as follows: in Sec. II we introduce our model and the basic many-body formalism necessary for deriving a set of self-consistent equations for the Green and vertex functions which are the variational parameters of the theory. The complete thermodynamics is then determined by integrals of the momentum and frequency dependent Green functions. It is shown that with a modified coupling constant, the theory can be formulated in a way consistent with Ward identities, which guarantees a gapless Bogoliubov-Anderson mode for arbitrary strength of the coupling. In Sec. III we discuss the numerical solution, providing quantitative results for the critical temperature, the pressure, internal energy and the entropy of the BCS-BEC crossover both in the normal and superfluid phase. They are compared both with experimental and theoretical results based on numerical and field-theoretic approaches. Finally in Sec. IV we give a brief summary, and indicate open problems.

II. A MANY-BODY THEORY OF RESONANTLY INTERACTING FERMIONS
In order to describe interacting fermions near a Feshbach resonance, it is in general necessary to include the resonant, closed channel bound state explicitely, e.g. within a Bose-Fermi-resonance model [33,34]. As has been shown for instance by Diener and Ho [11], however, the situation can be simplified in the case of broad Feshbach resonances, where the effective range r ⋆ of the resonant interaction is much smaller than both the background scattering length a bg and the Fermi wavelength λ F . In this limit, which is in fact appropriate for the existing experimental studies of the BCS-BEC crossover problem in 6 Li [5] and in 40 K [1], the problem can be reduced to a single channel Hamiltonian with an instantaneous interaction [10,11,12,13]. The associated effective two-body interaction is thus described by a pseudopotential V (r) ∼ δ(r) (appropriately renormalized, see below) with a strength proportional to the scattering length Here a bg is the off-resonant background scattering length in the absence of the coupling to the closed channel while ∆B and B 0 describe the width and position of the resonance which may be tuned by an external magnetic field B. The interacting Fermi system is thus described by the standard Hamiltonian [∇ψ + σ (r)][∇ψ σ (r)] where ψ σ (r) and ψ + σ (r) are the usual fermion field operators. The formal spin index σ labels two internal degrees of freedom, which in practice are two different hyperfine states. In the approximation, where the effective range of the resonant interaction is taken to zero, the interaction potential can formally be replaced by a delta potential between fermions of opposite spin Its strength g 0 needs to be renormalized for dimensions d ≥ 2 by introducing the scattering amplitude g via For dimensions d ≥ 2 the integral diverges at high momenta. Since the scattering amplitude g is kept constant, the bare interaction parameter g 0 must be taken to zero in the limit where the cutoff diverges. The associated limiting process g 0 → −0 accounts for the replacement of the bare potential (2.3) by a pseudopotential with the proper scattering length. While the formulas are derived for arbitrary space dimensions d, eventually we consider fermions for d = 3. In this case the scattering amplitude g is simply connected to the s-wave scattering length a given in (2.1) by g = 4π 2 a/m. In the following, we consider a homogeneous situation described by a grand canonical distribution at fixed temperature and chemical potential. The thermodynamic properties thus follow from the grand partition function and the associated grand potential which is directly related to the pressure p via Ω = −pV . Within our simplified model, where the range of the interaction is set to zero, the Fermi system is described by three parameters: the temperature T , the chemical potential µ and the s-wave scattering length a. Apart from an overall scale, the thermodynamics thus depends only on two dimensionless ratios. It is convenient to replace the chemical potential µ by the fermion density n = k 3 F /3π 2 , which defines the Fermi wave number k F and the Fermi energy ε F = 2 k 2 F /2m as characteristic length and energy scales. The equilibrium state is then uniquely determined by only two parameters: the dimensionless temperature θ = T /ε F (we choose units for the temperature in which k B = 1), and the dimensionless interaction strength v = 1/k F a. In the special case B = B 0 of an infinite scattering length (the so-called unitarity limit), the parameter v drops out and the resulting thermodynamic quantities are universal functions of θ [14].

A. Luttinger-Ward formalism
The BCS-BEC crossover is controlled by two physical phenomena. The first one is connected with the formation of pairs due to the attractive interaction. The second one is the transition to superfluidity below a certain critical temperature T c . In the BCS-limit, the formation of pairs and the superfluid transition are simultaneous. The transition is driven by the thermal breakup of pairs, i.e. by excitations which may be described by a purely fermionic theory. With increasing strength of the interaction, however, there is an increasingly wide range of temperatures where bound pairs coexist with unpaired fermions. In the BEC-limit, pair formation, as a chemical equilibrium between bound and dissociated atoms, occurs at a temperature scale much higher than the superfluid transition. The latter is driven by collective excitations of a then purely bosonic system. A proper description of the crossover thus requires to account for both bosonic and fermionic excitations simultaneously. Following the formalism developed by Luttinger and Ward [35] for non-superfluid interacting Fermi systems, the grand thermodynamic potential (2.6) can be expressed as a unique functional of the Green function in the form The trace Tr is defined with respect to the formal index X = (r, τ, σ, α) which combines the space variable r, the imaginary time τ , the spin index σ, and the Nambu index α. The interaction between the fermions is described by the functional Φ[G], which can be expressed in terms of a perturbation series of irreducible Feynman-Diagrams where the propagator lines are dressed and identified by the matrix Green function G of (2.7). While the formalism of Luttinger and Ward was originally derived for normal quantum liquids, it is well suited also to describe superfluid systems. Indeed the nondiagonal elements of the matrix Green function G represent the order parameter of the superfluid transition. The minimization of the grand potential Ω[G] as a functional of the Green function G thus incorporates the standard thermodynamic criterion that the order parameter is found by minimizing the thermodynamic potential. The stationarity condition uniquely determines the full matrix Green function G of the interacting system and hence the order parameter. It is important to note, that the thermodynamic potential Ω[G] depends on the exact Green function G. The formalism of Luttinger and Ward thus leads via (2.9) to a self-consistent theory for the matrix Green function G.
Since the Green functions contain information about the full dynamical behaviour via the imaginary time dependence of the Matsubara formalism, the Luttinger-Ward approach not only provides results for the equilibrium thermodynamic quantities but also determines spectral functions and transport properties. In our present work, however, dynamical properties will not be discussed. The functional Φ[G] is defined by an infinite perturbation series of irreducible Feynman diagrams and an exact expression for Φ [G] is clearly beyond what can be done analytically. An approximation which properly describes the formation of pairs, is a ladder approximation [36]. In Fig. 1, the related diagrams of Φ[G] are shown. The ladder approximation is self consistent because the propagator lines are dressed lines which are identified by the matrix Green function G. In the weak coupling BCS regime the ladder approximation becomes exact. For very strong attractive interactions, well above the pairing threshold, the fermion system is a Bose liquid of dilute atom pairs. In this limit the ladder approximation describes the formation of pairs (two-particle problem) exactly, however the interaction between the pairs (four-particle problem) only approximately [20,37]. In particular the resulting dimer-dimer scattering length is given by the Born approximation a (B) dd = 2a.

B. DeDominicis-Martin formalism
An extension of the Luttinger-Ward formalism was given by DeDominicis and Martin [38]. They introduce up to four external fields, which couple to products of one, two, three, and four field operators, and perform the Legendre transformations to the corresponding conjugate variables -the Green functions. For fermion systems only two external fields are relevant which couple to even products of fermion field operators. The related two conjugate variables of the Legendre transformation are the one-particle Green function G and the two-particle Green function G 2 . Within our approach below, the second Legendre transformation is performed explicitely. A more convenient conjugate variable is then the vertex function Γ which is related to G 2 by (2.15) below. Thus, DeDominicis and Martin obtain a thermodynamic potential which is a functional of both G and Γ. More precisely, it turns out that the relevant functional is the entropy S = F (2) where (2.10) (see (61) in the second paper of Ref. 38 and identify G 1 = G, C 2 = −Γ, andC 2 = −Γ therein).Γ is defined in (2.14) below.
The formalism of DeDominicis and Martin is ideally adapted to describe the BCS-BEC crossover because it explicitly deals with the one-particle Green function G, which represent the properties of the single fermions, and the vertex function Γ, which describes the eventually purely bosonic properties of the fermion pairs (both condensed or noncondensed). In particular, a full implementation of their formalism is needed to correctly account for four particle correlation, which is necessary to obtain the exact result a dd = 0.60 a for the dimer-dimer scattering length in the BEC limit.
As in standard thermodynamics, the entropy (2.10) is maximized under the constraints that all conserved quantities are kept constant. For the interacting fermion system defined by the Hamiltonian (2.2) the conserved quantities are the internal energy U = Ĥ and the particle number N = − 1 2 Tr{G}. Evaluating the thermal average of the Hamiltonian (2.2) we find that U can be expressed in terms of G and Γ (see (2.15) and (2.24) below).
Consequently, the entropy F (2) [G, Γ] = S[G, Γ], the internal energy U [G, Γ], and the particle number N [G] are functionals depending on G and Γ. In order to find the maximum of the entropy under the constraint of given average values of the particle number and the internal energy, DeDominicis and Martin [38] consider the functional (2.11) where λ U and λ N are two Lagrange parameters for the two constraints. Alternatively and equivalently, we consider the functional which is the grand thermodynamic potential where the temperature T and the chemical potential µ are the Lagrange parameters. Both functionals (2.11) and (2.12) must be stationary under small variations of G and Γ. In this way, we obtain the stationarity criteria which uniquely determine the one-particle Green function G and the vertex function Γ. In order to simplify the second trace in the entropy functional (2.10) it is convenient to define a modified vertex functionΓ bȳ where the four external propagator lines are amputated only half way (see (46) in the second paper of 38). For a proper definition of the second trace and the related matrix products the four indices of the modified vertex function must be grouped into pairs according tō Γ =Γ (X1X2)(X3X4) . The last term in (2.10), the functional K (2) [G, Γ] (depicted in Fig. 2), is defined by an infinite perturbation series of 2-line irreducible Feynman diagrams, where the propagator lines and the vertices are dressed and identified by the one-particle Green function G and by the vertex function Γ, respectively.
In order to understand the physical meaning of the various contributions to the thermodynamic potential, we note that the Luttinger-Ward formalism and the DeDominicis-Martin formalism are related to each other by a Legendre transformation, in which the bare twoparticle interaction as an external field is transformed into the two-particle Green function G 2 . This Legendre transformation may be interpreted as a renormalization procedure. Since the two-particle Green function G 2 is expressed in terms of the vertex function Γ by the bare two-particle interaction is replaced by the vertex function Γ as the renormalized interaction, which is the many-particle generalization of the scattering amplitude g.
We may compare the functionals (2.8) and (2.10) directly with each other. The first trace in (2.10) is identified by the trace in (2.8), which describe the contribution of single particles to the grand canonical potential and entropy, respectively. The second trace and the functional K (2) [G, Γ] in (2.10) which represent the interaction terms are related to the functional Φ[G] in (2.8). By close inspection (see (62) in Ref. 38) we find that the second trace in (2.10) represents the inverted perturbation series of ladder diagrams. It includes both particle-particle and particle-hole ladders, which describe the scattering and formation of pairs and it also includes bubble diagrams, which describe the screening of the interaction.
These three types of diagrams and mixtures of them arise because the vertices are symmetrized so that each of them can be expressed as a sum of three unsymmetrized vertices. As a result, the self-consistent ladder approximation of the functional Φ[G] shown in Fig. 1 can be reformulated within the formalism of DeDominicis and Martin in the following way: the second trace is approximated by keeping only the particle-particle ladders and the complicated functional K (2) [G, Γ] is set equal to zero. This neglects the screening of the interaction due to particle-hole excitations (see Subsec. III A below) and also the coupling between collective excitations and bound pairs.
In the following subsections we employ the formal-

C. Thermodynamic potentials
The formalism of Luttinger and Ward [35] allows to calculate directly the grand thermodynamic potential Ω. The functional Φ[G] has been evaluated explicitly in Ref. 31. Inserting this result into (2.8) we obtain (2. 16) In this formula the matrix Green function is defined by Knowledge of the matrix Green functions determines the matrix pair propagator via In order to distinguish between fermionic and bosonic functions the fermionic wave vectors and Matsubara fre-quencies are denoted by small letters, while the bosonic wave vectors and Matsubara frequencies are denoted by capital letters. In the second term of (2.16) the anomalous Green function is identified by F (0, 0) = F (r = 0, τ = 0). The formulas are derived for an arbitrary dimension of space d. The volume is assumed to be a cube with edge length L and periodic boundary conditions, where the limit L → ∞ is taken. The strength of the attractive interaction is included by the bare interaction parameter g 0 . The kinetic energy of the atoms ε k = 2 k 2 /2m and the chemical potential µ are implicitly included via the free matrix Green function G 0 (k, ω n ), which is related to the free normal Green function (2. 19) and the free anomalous Green function by a formula which is analogous to (2.17). The tem-perature T is included explicitly by the factors 1/β and implicitly by the Matsubara frequencies ω n and Ω n .
As evident from (2.16), the formalism of Luttinger and Ward, though including the exact single-particle Green function, still contains the bare coupling constant g 0 . In the DeDominicis-Martin formalism the bare coupling is renormalized and replaced by the exact vertex function Γ via a second Legendre transformation. The corresponding functional S[G, Γ] = F 2 [G, Γ], is just the dimensionless entropy as given in ((2.10)) As discussed above we restrict the second trace in (2.10) to the particle-particle ladders and by the nature of our interaction potential (2.3) to s-wave scattering. Furthermore, we omit the 2-line irreducible Feynman diagrams by setting K (2) [G, Γ] = 0. This approximation covers the essential features of the crossover problem namely the formation of pairs and their condensation. Within our ladder approximation, the DeDominicis-Martin formalism thus leads to an expression for the entropy of the form (2.21) The first term is clear. It is directly obtained from the first trace in (2.10). However, the second term resulting from the second trace in (2.10) needs further explanation. From (2.14) and the definition of the pair propagator (2.18) we infer (2.22) The reduction to particle-particle ladders implies that the Nambu indices are pairwise identical. In this way, the four Nambu indices of the vertex function Γ reduce to two Nambu indices. As a result, the vertex function Γ(K, Ω n ) = (Γ αα ′ (K, Ω n )) is a 2×2 matrix in the Nambu space similar to the matrix Green function (2.17). For the formalism of Luttinger and Ward the reduction of the vertex is described in detail in Ref. 31 and also in Ref. 29. Since the second trace of (2.10) is reduced to the particle-particle ladders and the structure of the vertex function is simplified considerably due to s-wave scattering, the prefactors of the terms in the second trace of (2.21) are changed. The factor 1 2 in front ofΓ disappears. Furthermore, the quadratic terms in the second trace cancel.
Another important quantity to consider is the internal energy U . With the help of the delta potential (2.3) we find for the expectation value of the Hamiltonian (2.2), (2.23) The second term contains an average of four fermion field operators which can be expressed in terms of the twoparticle Green function G 2 . Following the formalism of DeDominicis and Martin [38] and using (2.15) the twoparticle Green function can be expressed by four terms. The first three terms represent the three possibilities to factorize the two-particle Green function into products of two one-particle Green functions according to the Wick theorem. These terms provide the Hartree energy, the Fock energy, and the Bogoliubov energy. The fourth term is the connected part of the two-particle Green function and provides the correlation energy. Taking all terms together we obtain the internal energy (2.24) Finally, the particle number N is defined by the average which can be expressed in terms of the normal Green function in the standard form The entropy (2.21), the internal energy (2.24), and the particle number (2.26) are the basic functionals of the formalism of DeDominicis and Martin. The entropy S[G, Γ] is maximized under the constraints that the internal energy U [G, Γ] and the particle number N [G] are constant. In order to do this, the grand thermodynamic potential Ω[G, Γ] is defined by (2.12) where the temperature T and the chemical potential µ are Lagrange parameters. The self-consistent equations for the Green function G and the vertex function Γ are obtained from the stationarity conditions (2.13). Formally, the formalism of DeDominicis and Martin yields a different expression for the grand thermodynamic potential Ω than the formalism of Luttinger and Ward does by (2.16). However, it can be shown that the results are identical if and only if G and Γ satisfy the self-consistent equations (see end of Subsec. II E).
The functionals (2.21), (2.24), and (2.26) do not depend explicitly on the thermodynamic parameters T and µ. While the temperature appears explicitely via the factor β = 1/T and the Matsubara frequency ω n ∼ Ω n ∼ T , a proper rescaling of the functions G → β G, χ → β χ, and Γ → Γ implies that all factors β and T cancel in all three functionals. The temperature T and the chemical potential µ enter only as Lagrange parameters via the constraints. This fact is a general property of the formalism of DeDominicis and Martin. The fermion mass m, the kinetic energy ε k = 2 k 2 /2m, and the interaction parameter g 0 , which determine the microscopic properties of the interacting fermion system, are present only in the internal energy functional (2.24).
An alternative expression for the entropy is obtained from the grand thermodynamic potential of Luttinger and Ward (2.16) according to the standard thermodynamic relation Taking the partial derivative we obtain an expressions which formally differs from (2.21). However, provided that G and Γ satisfy the self-consistent equations, the results for the entropy will be identical. Therefore, both the Luttinger-Ward and the DeDominicis-Martin formalism exactly obey all the standard thermodynamic relations provided the Green functions obey the stationarity conditions (2.9) and (2.13). The equivalence of the different formal expressions in thermal equilibrium is very important for the consistency of our theory and the compatibility of the self-consistent ladder approximation for all thermodynamic quantities. Apart from the entropy, we can also determine the pressure p = −Ω/L d as a functional of the Green function G using (2.16) or (2.12). The dimensionless thermodynamic quantities Ω/N ε F , U/N ε F and S/N will be calculated numerically in Sec. III and discussed in the following sections.

D. Self-consistent equations for the Green and vertex functions
The self-consistent equations for the Green functions follow directly from the stationarity condition (2.9). Inserting the general functional of the Luttinger-Ward formalism (2.8) into this condition we obtain the Dyson equation The self energy Σ is identified by the functional derivative The functional Φ[G] is defined by a perturbation series. The related Feynman diagrams are shown in Fig. 1 for the self-consistent ladder approximation. Inserting the grand thermodynamic potential (2.16) into the constraint (2.9), we obtain an explicit expression for the self energy which is In the first term δ F (τ / ) is the fermionic delta function which is antiperiodic. The order parameter of the superfluid transition is represented by the nondiagonal elements of the matrix In the second term of (2.30) Γ is the matrix vertex function, which is related to the matrix pair propagator χ by Eventually, χ is represented in terms of the matrix Green function by (2.18). Eq. (2.33) is just the Bethe-Salpeter equation in ladder approximation. It is responsible for the fact that the binding of fermion pairs is described appropriately. Taken together, we have now a set of selfconsistent equations for the matrix Green function G and the matrix vertex function Γ which have to be solved numerically.
Alternatively we can derive the self-consistent equations by inserting the functional (2.12) of the formalism of DeDominicis and Martin into the related stationarity conditions (2.13). We obtain the Dyson equation (2.28) from the first condition and the Bethe-Salpeter equation (2.33) from the second condition. In this way we prove that both the Luttinger-Ward formalism and the DeDominicis-Martin formalism are equivalent within our approximation.
Unfortunately, in the present form, the matrix pair propagator χ defined in (2.18) is divergent. While the sum over the Matsubara frequencies is finite, the integral over the wave vector is ultraviolet divergent for dimensions d ≥ 2. For this reason a renormalization is necessary. We define the regularized pair propagator by Inserting this formula into (2.33) we obtain the renormalized Bethe-Salpeter equation The bare interaction strength g 0 is renormalized according to (2.4) and replaced by the scattering amplitude g. For d = 3 dimensions g is expressed in terms of the s-wave scattering length a by g = 4π 2 a/m. The zero range of the interaction between the fermions implies that is infinite. For this reason the order-parameter formula (2.31) must be renormalized, too. Replacing the bare interaction strength g 0 by the scattering amplitude g according to (2.4), we obtain the renormalized formula Here, the integral over the wave vector is finite.

E. Reformulation in terms of mean-field Green functions
In mean-field approximation the self energy Σ(k, ω n ) is replaced by Σ 1 defined in (2.32). Since Σ 1 depends neither on wavevector nor on frequency, the approximation Σ ≈ Σ 1 just describes the formation of a pair condensate within a BCS-type mean-field theory where the destruction of superfluidity is driven by the breakup of pairs. This is the correct description in the weak coupling limit, however for strong coupling the superfluid transition is driven by finite momentum pairs whose contribution is contained in the second term of the self energy (2.30). Inserting the mean-field self energy into the Dyson equation (2.28) we obtain where G 1 is the matrix Green function in mean-field approximation.
If we consider the self-consistent equations and the formulas for the thermodynamic potentials we realize that the spectrum ε k of the fermionic atoms and the chemical potential µ enter the formulas only implicitly via the free matrix Green function G 0 . We can transform the formulas so that G 0 is replaced in favor of the mean-field matrix Green function G 1 . As a result we obtain the Dyson equation is the second term of the self energy (2.30). The other self-consistent equations remain unchanged. The grand thermodynamic potential (2.16) is transformed into (2.41) For a combination of the internal energy (2.24) and the particle number (2.26) we obtain the formula (2.42) The entropy (2.21) depends neither on G 0 nor on G 1 . We transform the formula into (2.43) In the above three formulas we have simplified the terms involving the vertex function Γ by using the Bethe-Salpeter equation (2.33). The grand thermodynamic potential (2.41) was derived using the formalism of Luttinger and Ward [35] while the other two quantities (2.42) and (2.43) were derived using the formalism of DeDominicis and Martin [38]. It is now not hard to see that the above expressions obey the thermodynamic relation which explicitely shows that both formalisms are indeed equivalent yielding the same results for all thermodynamic potentials in self-consistent ladder approximation provided G and Γ satisfy the appropriate stationarity equations.

F. Mean-field approximation
If we insert G 1 for the matrix Green function G and neglect all terms containing the vertex function Γ we obtain the thermodynamic potentials in mean-field approximation. In particular the mean-field grand thermodynamic potential is given by while the mean-field formula for the combination of the internal energy and the particle number are 46) and the mean-field entropy is In order to obtain finite results, we must define the sums over the Matsubara frequencies as described in Appendix A. The sums can be evaluated explicitly. This yields which are the well known results of a BCS variational Ansatz for arbitrary coupling where is an energy constant which after renormalization g 0 → g (see (2.57)) reduces to the BCS condensation energy.
Here E k is the spectrum of the quasiparticles, defined by and n k denotes the Fermi distribution function of the quasiparticles We find that the regularization of the Matsubara-frequency sums described in Appendix A affects only the energy constant E 0 . The regularization has been chosen such that for zero interaction the results for the ideal Fermi gas are obtained which implies E 0 = 0. The other terms in (2.48)-(2.50) are not affected by the regularization.

G. Beyond mean-field
In the mean-field approximation, the formation and condensation of fermion pairs occur at the same temperature. This is the well known BCS scenario, which is perfectly captured by the exact solution of the reduced BCS-Hamiltonian. Formally, this solution can easily be extended to arbitrary coupling strengths [23]. At zero temperature, it provides a smooth crossover from the BCS groundstate of highly overlapping pairs to a perfect Bose-Einstein condensate at infinite coupling, similar to the variational Ansatz of Eagles and Leggett [21,22]. At finite temperature, however, superfluidity in this model is destroyed by fermionic excitations, namely the breakup of pairs. The critical temperature is therefore of the same order as the pairing gap at zero temperature, consistent with the well known BCS relation 2∆ 0 /T c = 3.52 in weak coupling. Clearly, such a picture is appropriate for weak coupling, where the transition to superfluidity is driven by the gain in potential energy associated with pair formation. By contrast, for sufficiently strong interactions, the superfluid to normal transition is instead driven by a gain in kinetic energy, associated with the condensation of pre-formed pairs rather than their thermal breakup. The critical temperature is then of the order of the degeneracy temperature of the gas and thus is completely unrelated to the pair binding energy. For a proper description of the BCS-BEC crossover at finite temperature and arbitrary coupling, we therefore need to go beyond mean-field, including excitations, which drive the superfluid order parameter to zero without destroying the bound pairs altogether. On a formal level, this is accomplished by the nontrivial wave-vector and frequency dependent term GΓ in the exact fermion self energy, as given in (2.30).
For the numerical calculation we decompose the thermodynamic potentials into a mean-field part and a correction term according to Ω = Ω 1 + ∆Ω, S = S 1 + ∆S, etc.. The mean-field contributions have been derived in the previous subsection. While in these contributions the Matsubara-frequency sums have been performed explicitly, the integrals over the wave vector remain and must be evaluated numerically. By subtracting the mean-field formulas (2.45)-(2.47) from the general formulas (2.41)-(2.43) we obtain the correction for the grand thermodynamic potential the correction for the combination of the internal energy and the particle number 55) and the correction for the entropy (2.56) In formulas (2.54)-(2.56) the sums over the fermionic Matsubara frequencies ω n converge so that the regularization of Appendix A is not needed. However, the sums over the bosonic Matsubara frequencies Ω n are not well defined and must be regularized. Thus, for a numerical evaluation the formulas (2.54)-(2.56) must be transformed further, which will be done in the next subsection.

H. Renormalization of the thermodynamic potentials
Since the interaction has zero range, the interaction strength g 0 must be renormalized and replaced by the scattering amplitude g according to (2.4). In a first step we renormalize the mean-field formulas of the thermodynamic potentials. In (2.48)-(2.50) the interaction strength g 0 does not occur explicitly. The integrals are thus finite and no renormalization is needed for these formulas. However, the condensation energy (2.51) contains two infinite terms, a divergent integral and the last term with the infinite factor 1/g 0 , which compensate each other. By renormalizing the interaction strength we obtain where both the integral and the last term are now separately finite. Note that the wave vector integrals in (2.48)-(2.50) and in (2.57) are finite in any spatial dimension d with 2 < d < 4.
In a second step we renormalize the correction formulas. In the correction of the grand thermodynamic potential (2.54) we decompose ln[Γ(K, Ω n )/g 0 ] = ln[Γ(K, Ω n )/g] + ln[g/g 0 ]. The separated term ln[g/g 0 ] can be neglected because it does not depend on the Matsubara frequencies Ω n . Following the arguments of Appendix A the Matsubarafrequency sum of this term is zero. Thus, for the correction of the grand thermodynamic potential we obtain the renormalized formula (2.58) Both terms of this formula are now finite. However, Eq. (A4) is needed to evaluate the second term.
In correction (2.55) the second term must be renormalized. This can be achieved by the following sequence of equations (2.59) First, by using the Bethe-Salpeter equation (2.33) we write the integrand as a product of the matrix vertex function Γ and the matrix pair propagator χ. Secondly, we express χ in terms of the matrix Green function G by (2.18), interchange the orders of the integrals and sums, and combine Γ with one of the G into the self energyΣ by (2.40). Finally, we replaceΣ in favor of G and G 1 by using the Dyson equation (2.39). The bosonic integral and sum are transformed into a fermionic integral and sum. Hence, the second term of (2.55) can be combined with the first term. As a result we finally obtain (2.60) By considering (2.38) we explicitly prove Consequently, for the correction of the combination of the internal energy and the particle number we obtain the compact formula which is essential for a stable numerical evaluation of the correction term. The Matsubara-frequency sum is evaluated by using (A3). The wave-vector integral is finite. The correction of the entropy (2.56) is renormalized in an analogous way. Alternatively, we use the thermodynamic relation (2.44). As a result we obtain (2.63) The final results are the mean-field formulas (2.48)-(2.50) together with (2.57) and the correction formulas (2.58), (2.62), and (2.63). In these formulas each term by itself is finite. Eventually, the thermodynamic potentials are obtained by adding the terms together according to Ω = Ω 1 + ∆Ω, S = S 1 + ∆S, etc..

I. Symmetry breaking and Thouless criterion
The interacting fermion system is invariant under the symmetry transformation which is related to a global change of phase of the fermion fields by λ. The superfluid phase breaks this symmetry since the order parameter ∆ is transformed as ∆ → e 2iλ ∆. Clearly, however, the thermodynamic potentials must remain invariant under a global change of the phase both in the normal and in the superfluid state. In the superfluid, the free energy increase associated with a slowly varying phase λ(r) vanishes like (∇λ) 2 . By Goldstone' s theorem, this implies the existence of modes whose energy vanishes in the long wave-length limit. For a neutral superfluid, this is the well known Bogoliubov-Anderson mode. It has a sound like dispersion ω(k) = ck and is physically related to fluctuations of the phase of the order parameter. In technical terms, the existence of zero-energy collective modes can be derived from Ward identities related to the symmetry transformation. By considering the grand thermodynamic potential Ω[G], in Ref. 31 the Ward identity has been derived (see (2.57) in Ref. 31). Here δ λ Σ XX ′ is the variation of the self energy under the transformation (2.64) with an infinitesimal phase change δλ. This quantity may be interpreted as the generalized order parameter of the system. On the other hand, Γ −1 XX ′ ,Y Y ′ is the inverse vertex function. For a short-hand notation the indices X, X ′ and Y , Y ′ are used, which represent a combination X = (r, τ, σ, α) of the space variable r, the imaginary time τ , the spin index σ, and the Nambu index α. According to (2.65) the inverse vertex function Γ −1 may be interpreted as a linear operator which acts on the order parameter δ λ Σ. In the superfluid state the order parameter is nonzero so that the inverse vertex function must have a zero eigenvalue, which is related to a zero energy collective mode. For superfluid fermion systems this fact is known as the Thouless criterion [39].
The Ward identity (2.65) has been derived for the exact theory. However, our present crossover theory is an approximation, based on a certain truncation of the exact functional which enters either into the Luttinger-Ward or the DeDominicis-Martin formalism. In general, such a truncated functional will not obey the Ward identity. Indeed, we find that our inverse vertex function Γ −1 αα ′ (K, Ω n ) obeys instead the equation where (∆ α ) = (∆, ∆ * ) (see (3.56) in Ref. 31). Taking the longitudinal part, this equation correctly describes the smooth evolution from a Ginzburg-Landau type description of weak coupling BCS superfluids to a Gross-Pitaevskii like theory of a dilute, repulsive Bose gas [31]. The transverse part, however, also gives a finite value on the right hand side of (2.66) in the limit K = 0 and Ω n = 0, thus violating the Ward identity by terms of order |∆| 3 . As a result, the Thouless criterion is violated and there is no proper Bogoliubov-Anderson mode in our approach without a further modification (see below). Unfortunately, the violation of the Goldstone theorem for continuous symmetries is a general property of conserving approximations based on the Luttinger-Ward formalism. This problem has been known for a long time for superfluid Bose systems [40] and is sometimes referred to as the 'conserving-gapless dichotomy' [41,42] in the literature. For the exact theory, a Ward identity holds for the inverse matrix boson Green function G B , which reads In the superfluid state the inverse matrix boson Green function has a vanishing eigenvalue. For superfluid boson systems this is known as the Hugenholtz-Pines theorem [43]. Conserving approximations, however, violate the Hugenholtz-Pines theorem. For example, this is true already for the lowest approximation, the well known Hartree-Fock-Bogoliubov theory. In our fermion system for strong attractive interactions v = 1/k F a ≫ 1, the fermions are bound into pairs. These pairs form a Bose system with an effective repulsive interaction which, for a dilute system, is described by the exact scattering length a dd ≈ 0.60 a of the four particle problem associated with dimer-dimer scattering. In the strong coupling limit, therefore, our crossover theory for interacting fermions must converge to an effective theory of repulsively interacting bosons, where both theories are based on the Luttinger-Ward formalism. From the analytical arguments in Refs. 29, 31 and also from our numerical calculations, we find that the crossover theory converges to the Hartree-Fock-Bogoliubov theory quickly for interactions v = 1/k F a > 2. The boson order parameter Ψ B and the matrix boson Green function G B (K, Ω n ) can be identified with the order parameter ∆ and the vertex function Γ(K, Ω n ) according to [29,31] The validity of these relations requires both strong coupling, but also low frequencies and momenta. Indeed, it is only at low energies where the composite particles behave like bosons. At higher frequencies or momenta, the composite nature of the pairs becomes visible. This becomes evident, for instance, in the different behaviour G B ∼ Ω −1 n of a Bose Green function at large frequencies compared to that of the vertex function, which behaves like Γ ∼ Ω −1/2 n as a result of the two particle continuum associated with broken fermion pairs. Clearly, at large coupling constants v ≫ 1, this continuum moves up to very large frequencies of the order of the binding energy We conclude that the violation of the Thouless criterion in our crossover theory is related to the violation of the Hugenholtz-Pines theorem in the Hartree-Fock-Bogoliubov theory for bosons to which our Luttinger-Ward formulation of the fermionic many-body problem converges at large coupling. In the following section, it will be shown that this problem may be solved by an appropriate modification of the coupling constant. In this manner, a self-consistent formulation of the manybody problem is possible which obeys Goldstone's theorem and thus provides a correct description of both fermionic and collective, bosonic excitations along the BCS-BEC crossover.

J. Modified coupling and gapless Bogoliubov-Anderson mode
In the following, our aim is to modify the theory in a way which is consistent with the Thouless criterion, giving rise to a gapless Goldstone mode in the whole regime of coupling strengths. If we require the Thouless criterion a further equation will be added to the self-consistent equations for the Green and vertex functions in Subsec. II D. However, then another equation must be discarded or a further parameter must be introduced. We find that the Bethe-Salpeter equation (2.35) and the orderparameter equation (2.37) can not be satisfied together, if (2.70) is required. For this reason we modify the theory by introducing a modified scattering amplitude g mod , which is determined by the modified order-parameter We have solved the self-consistent equations together with the Thouless criterion (2.70) and the modified orderparameter equation (2.71). The numerical effort is much less for the modified theory than for the original theory. Since the scattering amplitude g is related to the scattering length a, we obtain a modified dimensionless interaction strength v mod = 1/k F a mod . We find a difference δv mod = v mod − v in the range between 0.0 and −0.1.
In order to obtain a consistent theory, we must check that the modification is compatible with the Luttinger-Ward formalism. We must find a modified grand thermodynamic potential Ω mod [G], so that the condition for stationarity (2.9) yields the self-consistent equations with the modified order-parameter equation (2.71). For this purpose we consider the second term of (2.16) which reads We replace this term by the modified term is the modified order-parameter equation. The modified interaction strengthsg 0,mod =g 0,mod (|∆|) and g 0,mod = g 0,mod (|∆|) depend on the order parameter |∆|, are not equal, and differ from the bare interaction strength g 0 . In order to apply the stationarity condition (2.9) we must consider the variation of (2.73) with respect to G. Since the modified parameterg 0,mod (|∆|) depends implicitly on G via (2.74), the chain rule of differential calculus must be applied. Eventually, the variation of (2.73) must have the form By comparing the resulting terms with (2.75), we obtain the differential equation . (2.76) On the other hand, Eq. (2.66) implies that the Thouless criterion holds without modification if |∆| = 0. Thus, we findg 0,mod = g 0,mod = g 0 for |∆| = 0 (2.77) which is an initial condition for (2.76). Eq. (2.76) can be integrated together with (2.77). We obtain .
(2.78) The thermodynamic state of the interacting fermion system in the superfluid state is therefore determined by three parameters. We may choose the order parameter |∆|, the chemical potential µ, and the interaction strength g 0 for these parameters. Hence, the modified interaction strengths g 0,mod = g 0,mod (|∆|, µ, g 0 ) and g 0,mod =g 0,mod (|∆|, µ, g 0 ) are functions of these parameters. While g 0,mod (|∆|, µ, g 0 ) is uniquely determined by (2.74) and the other self-consistent equations, g 0,mod (|∆|, µ, g 0 ) depends on the path in the parameter space when the integral (2.78) is calculated. Since g 0 and µ are external parameters of the theory, for a correct formulation of the modification these parameters must be kept constant.
The modification is compatible also with the DeDominicis-Martin formalism. In this case the internal energy U [G, Γ] includes the term (2.72) which must be modified according to (2.73). The modification of the coupling constant g 0 described by (2.76)-(2.78) is derived in an analogous way.
Eqs. (2.72)-(2.78) describe the modification of the crossover theory in terms of the bare interaction parameters g 0 , g 0,mod , andg 0,mod . A renormalized version of the modification is obtained, if we replace the bare parameters by the renormalized scattering amplitudes g, g mod , andg mod according to (2.4). Eqs. (2.76)-(2.78) are valid also for the renormalized scattering amplitudes without changes. From (2.78) we obtain . (2.79) The renormalized modified order-parameter equation is defined by (2.71). In order to obtain the modified formulas of the renormalized thermodynamic potentials in Subsec. II H only a single change is needed. We must replace the energy constant (2.57) by .  While δv mod is obtained directly from (2.71) by solving the self-consistent equations, δṽ mod is obtained by evaluating the integral in (2.81) numerically. As a result we obtain modifications which are restricted to the interval − 0.1 δv mod < δṽ mod < 0 for |∆| > 0 . (2.82) In Fig. 3 the modifications δv mod and δṽ mod are shown as red solid curve and blue dashed curve, respectively, for T = 0 and |∆| = |∆ 0 |. Clearly, the modifications are largest in the crossover region close to the unitarity point. At finite temperature for increasing T the order parameter |∆| and the modifications δv mod and δṽ mod decrease together. Eventually, for |∆| = 0 the modifications are δv mod = δṽ mod = 0.
In the previous subsection we have argued that for strong attractive interactions the fermions are bound into pairs. Our crossover theory for the interacting fermion system then converges to a Luttinger-Ward type theory for interacting bosons. Since the modified version of our theory obeys the Ward identity, its strong coupling limit necessarily leads to a description of dilute, repulsive bosons which has the correct linear spectrum of excitations at low energies. It turns out, that the limiting theory here is the Luttinger-Ward version of the socalled Shohno theory [44,45] which is equivalent to the more well known Popov approximation. While we have not been able to derive the Shohno-Popov theory analytically from the Luttinger-Ward functional of the original fermionic model, we find a quick convergence numerically in all thermodynamic quantities for dimensionless couplings v > 2. Considering the entropy, in particular, the Shohno-Popov theory gives rise to the standard expression for the entropy of a noninteracting gas of bosonic quasiparticles with the standard distribution function The corresponding spectrum of excitation energies (2.85) has the well known form of a Bogoliubov spectrum with a temperature dependent condensate density n B,0 = |Ψ B | 2 and a positive Bose-Bose scattering amplitude g B . Within our approximation, we have g B = 2g, i.e. the exact dimer-dimer scattering length a dd ≈ 0.60 a is replaced by its Born approximation result a (B) dd = 2a [46]. The effective mass and chemical potential take their obvious values m B = 2m and µ B = 2µ + ε b where ε b = 2 /ma 2 is the two-particle binding energy on the BEC-side of the crossover, where a > 0. The order parameter Ψ B is given by (2.68). Other thermodynamic quantities are obtained using more complicated expressions, which are not presented here in detail [47]. As will be shown explicitely in the following sections, the numerical results for the critical temperature or the entropy converge quickly to that of the Shohno-Popov theory for coupling strengths v > 2.

III. NUMERICAL RESULTS
Following the detailed discussion of the formalism used to describe the thermodynamics of attractively interacting fermions at arbitrary coupling and temperature, we now present numerical results which cover both the normal and superfluid regime. These results require a solution of the self-consistent equations determining the Green function G and the vertex function Γ, which are scalars above, and two-by-two matrices below the critical temperature. An iteration procedure is performed where a numerical Fourier transformation is needed to transform the functions back and forth. Since the Green function G, the vertex function Γ, and the related func-tionsΣ and M are singular at small values of r and τ and also exhibit significant variation over several orders of magnitude, the numerical Fourier transformation is quite challenging. In practice, the variables need to be discretized on logarithmic scales. Standard procedures like the fast Fourier transformation are therefore not applicable. The basic principles of our special numerical Fourier transformation are described in Appendix B.

A. Critical temperature
The crucial quantity which determines the overall structure of the phase diagram is of course the critical temperature T c for the transition to a superfluid. This temperature is known analytically only in the extreme BCS-and BEC-limit. In the BCS-limit k F |a| → 0, where the average distance between the fermions is much larger than the magnitude of the scattering length, the standard solution of the gap-equation for an attractive pseudopotential gives a critical temperature is exponentially small on the characteristic scale of the Fermi energy. Since typical Fermi temperatures in cold gases are of the order of micro-Kelvin, the BCS-regime is in practice hardly attainable in these systems.
The leading order corrections to the BCS-result in an expansion in the small parameter k F |a| ≪ 1 have been determined a long time ago by Gorkov and Melik-Barkhudarov [48]. They arise from induced interactions, where one fermion sees the polarization in the Fermi gas due to a second fermion. The density induced interaction changes the dimensionless coupling constant N (0)g = 2k F a/π of the BCS-theory to [49] to where N (0) = mk F /2π 2 2 is the standard density of states per spin at the Fermi energy. Since the additional contribution to the two-body scattering amplitude g < 0 is positive, the induced interactions weaken the attractive interaction between two fermions in vacuum and lead to a reduction of the transition temperature by a factor (4e) −1/3 ≈ 0.45. The nonanalytic dependence of the BCS-transition temperature on the dimensionless coupling constant k F a thus give rise to a finite change in the prefactor in (3.1) from the BCS value 0.61 to 0.28, even though the contribution of induced interactions is of order k F a compared to the bare interaction. On the BEC-side, the zeroth order result for the critical temperature is obtained from the value obtained for an ideal Bose gas with density n B = n/2 and mass m B = 2m. The leading corrections to this result arise from the residual interactions between the strongly bound bosonic dimers. As shown by Petrov et al. [20,50], these interactions can be described by a positive dimerdimer scattering length a dd ≈ 0.60 a. With the quite plausible assumption, that the total potential energy in a dilute gas of dimers is the sum of its two-body interactions, the scattering length of the four fermion problem determines the corresponding interaction constant in the theory of a weakly interacting Bose gas in the regime of a small gas parameter n 1/3 B a dd ≪ 1, where Bogoliubov theory is applicable. The exact dependence of the critical temperature of the dilute, repulsive Bose gas on the interaction strength has been calculated only in recent years. To lowest order in the interaction, the shift is positive and linear in the scattering length [51], with a numerical constant c ≈ 1.31 [52,53]. As a result, the evolution of the critical temperature in the homogeneous case as a function of the dimensionless coupling constant v = 1/k F a necessarily exhibits a maximum, since the asymptotic ideal Bose gas result is approached from above. Such a maximum has been found in the early calculations of T c along the BCS-BEC crossover by Nozieres and Schmitt-Rink [54] and by Randeria et al. [55]. The precise height and location of this maximum, however, has not been determined so far in a quantitatively reliable manner. Given that our present theory exhibits a first order transition, there is a range of multivaluedness of the thermodynamic potentials as a function of temperature. This regime is bounded in Fig. 4 by the upper and lower T c curves respectively. The lower T c curve (shown as the red dashed line) which is monotonic in v coincides with the T c curve previously calculated [30] by implementing the Thouless criterion coming from the normal fluid side. In a situation, where a true first order transition is expected, we would need to perform a Maxwell construction to obtain the proper transition line. As was discussed above, however, the first order transition is an artefact of the approximations involved.
In particular the spectrum of excitations right at T c is free particle like in our approximation rather than ω K ∼ K 3/2 [56]. In order to determine the proper critical temperature within our approximation , we have used two essentially equivalent criteria: the fact that the exact entropy is continuous at T c suggests that our best approximation for the critical temperature is where the jump in the entropy between the two branches characterising the superfluid and the normal regime has a minimum. Essentially the same value is obtained by defining T c through the criterion that it is the maximum temperature at which the order parameter ∆(T ) is nonzero. Remarkably, these criteria lead to a critical temperature (shown as θ (upper) c in Fig. 4) which exhibits a maximum on the BEC-side of the crossover around v ≈ 1 as expected on general grounds. Moreover, our theory predicts the correct asymptotic functional form (3.4) of the T c -enhancement in the BEC limit v ≫ 1. Even though the dimer-dimer scattering length a (B) dd = 2 a and the prefactor c ≈ 0.58 of our approximate Popov-type theory differ from the exact values a dd = 0.60 a and c ≈ 1.31, respectively, the agreement of our theory with the exact result is very good (see Fig. 4).
A quite sensitive test of the quantitative reliability of our present result for the critical temperature at arbitrary coupling is provided by a comparison with the recent, rather precise numerical results right at the unitarity point by Burovski et al. [28]. In fact, our result for the dimensionless ratio T c /ε F ≈ 0.16 of the critical temperature in units of the bare Fermi energy, which is one of the universal numbers of the BCS-BEC crossover problem (see Subsec. III B below), agrees precisely with the numerical results of Burovski et al. within the given error bars. As will be shown below, a similar rather precise agreement is obtained with other thermodynamic quantities, except for the chemical potential. Thus, even in the absence of a small parameter which would allow to control our theory systematically in the crossover regime, the agreement with the numerical results at unitarity gives us confidence that the approach is quantitatively reliable at arbitrary coupling strengths. In Fig. 5 the temperature evolution of the entropy is shown for various coupling parameters v. Here the multivalued character is clearly seen which reflects the firstorder transition. Furthermore, three-dimensional plots of the entropy and of the pressure are presented in Figs. 6 and 7, respectively. In both figures a rather sharp drop is observed in the crossover region from weak coupling v −1 (fermionic regime) to strong coupling v +1 (bosonic regime). In the weak coupling limit v ≪ −1 the results of the nearly ideal Fermi gas are approached which are defined by the BCS formulas (2.48)-(2.50) and (2.57). On the other hand, in the strong coupling limit v ≫ +1 the results of Shohno's mean-field theory are approached. While the strong-coupling entropy is defined by (2.83), the other thermodynamic quantities are defined by more complicated formulae [47]. For v > 1.0 and very low temperatures the pressure is nearly zero which reflects a special property of Shohno's mean-field theory of weakly interacting bosons. At high temperatures T ≫ ε F the entropy, the pressure, and the related thermodynamic quantities approach the Boltzmann limit. to the well known BCS result for v ≪ −1.0. In the opposite limit of strong coupling the behaviour can be derived from µn → −∆ 2 /2g which reflects the fact that the fermion chemical potential in the strong coupling limit is governed by the potential (i.e. binding) energy. This yields ∆(T = 0)/ε F → (16/3π)v with the square root behaviour clearly visible in Fig. 8. Near T c the gap function displays the multivalued behaviour characteristic of a first-order transition.
At low temperatures the entropy has to vanish, in accordance with the third law of thermodynamics. The way it does, is in fact universal along the full BCS-BEC crossover. Indeed, at low temperatures, the twocomponent Fermi gas is in a superfluid state, independent of the strength of the attractive interaction. On quite general grounds therefore, the low lying excitations above the ground state are sound modes of the Bogoliubov-Anderson type. These modes give rise to an entropy which vanishes like T 3 for arbitrary coupling strength.
The associated sound velocity c is constant at low T and may be determined from the pressure via mc 2 = ∂p/∂n. Fig. 9 displays (c/v F ) 2 at T = 0 as a function of coupling strength with v F the Fermi velocity. The dilute interacting Fermi gas limit (c/v F ) 2 = (1 + 2/(πv))/3 and the BEC limit (c/v F ) 2 = k F a dd /(6π) for a dd = 0.60 a are represented by the blue squares and the green triangles respectively. The red triangles are obtained by extending the expression of the ground state energy of a dilute weakly interacting Fermi gas [57,58] with the help of a Pade approximation to the strong coupling regime [8,9]  and the thermodynamic identity Obviously the present crossover theory provides a very good description of the equation of state and sound velocity except in the regime v > 1, where our results underestimate both the pressure and its density dependence. In principle we should be able to independently obtain c from the low entropy asymptotics (3.5). Our numerical results are consistent with S(T ) ∼ T 3 , however they are not precise enough at such low temperatures, to extract the sound velocity in this manner.

B. Thermodynamics in the unitarity limit
After presenting the results for the critical temperature and the thermodynamics at arbitrary coupling, we now turn to a more detailed discussion of the unitarity limit, where the scattering length is infinite. This particular line in the phase diagram has received a lot of attention recently. In particular, precise numerical results are available at this point [28], which provide a sensitive test of analytical approaches to the crossover problem.
As has been mentioned before, the Fermi gas at infinite scattering length v = 0 is rather special since the only relevant length and energy scales remaining in the problem are the Fermi wave length set by the density and the Fermi energy ε F , provided we remain within the zero range pseudopotential approximation. The free energy has a simple scaling form which will be determined below, are the chemical potential and the internal energy in units of the Fermi energy or the entropy per particle at T c . In addition, also the gap for single particle excitations or the condensate fraction at zero temperature are universal at the unitarity point. Fig. 10 shows the temperature dependence of the internal energy calculated in two different ways. The solid line is our numerical result which is compared with the internal energy (depicted as the dashed line) as obtained from the numerically calculated pressure p = −Ω/V via the scaling relation U = 3pV /2 valid at the unitarity point. Our numerical results display perfect scaling above T c . The scaling violation below T c is a consequence of the modification of the theory. In order to preserve the conserving nature of our theory while obeying the Thouless criterion an extra length scale a mod had to be introduced leading to a modified dimensionless interaction strength v mod = 1/k F a mod with δv mod = v mod − v in the range between 0.0 and −0.1 with v mod = 0 for v = 0 (see Subsecs. II I and II J for details). Fig. 11 displays the behaviour of the chemical potential µ(T ) as a function of temperature. Using µ(T ) in a local density approximation µ = µ h n(r), T (r)/T F + V (r) (3.9) with µ h the chemical potential of the homogenous case we can calculate the density profiles of harmonically trapped ultracold gases at unitarity [59]. We have also checked the convergence of our µ(T ) to the high temperature expansion obtained by Ho and Mueller [60] which however only occurs for T ≫ ε F . Note that below T c the chemical potential µ(T ) is an increasing function of T . This perhaps counterintuitive result can be understood quite easily from the fact, that the low temperature thermodynamics is determined by the Bogoliubov-Anderson mode. As argued in the previous section, this leads to an entropy which vanishes with a power law ∼ T 3 . Fig. 12 displays the entropy at unitarity as a function of temperature. Now, at a given volume, there is a Maxwell relation of the form which connects the temperature dependence of the chemical potential to the density dependence of the entropy. Using the universal result (3.5) for the low temperature entropy, this relation shows that at low temperatures the chemical potential exhibits a T 4 dependence with a prefactor determined by ∂µ ∂T = 3S 2V mc 2 s ∂ 2 p ∂n 2 > 0 . (3.11) Obviously, this argument is not confined to the unitarity   [8,9]  Our value of β ∼ −0.640 is very close to β = −0.67 obtained via simply Pade approximating the weak coupling result for the ground state energy [8,9] [16,17] (see Table I). Evidently, there is still considerable uncertainty in both the experimental and theoretical results. A promising route in the direction of thermometry for trapped gases is provided via the reversible adiabatic (isentropic) sweeps [5,61] from the BEC limit. In Fig.  13 we depict the resulting changes in temperature when moving across the unitarity limit for the homogenous  [27]), ε = 4 − d expansion (Nishida and Son [24,26]), Borel-Padé approximation connecting an expansion in ε = 4 − d and one in ε = d − 2 [26]) and a 1/N expansion (Nikolić and Sachdev) [15] at T = Tc.  case. For the trapped case this cooling mechanism was first advocated by Carr et al. [32] and recently quantitatively refined by Hu et al. [66]. Finally to facilitate quantitative comparison with various Quantum Monte Carlo results we have collected available data from the literature presented in Table I and II. Apart from the value for T c which is explicitely quoted in the paper by Bulgac (with errors) we have estimated the remaining quantities from their presented results and utilized scaling to fill in the missing data below. The T = 0 results are fixed node QM results by Astrakharchik et al. [17] and Carlson et al. [16]. Note that our result for ∆/ε F is close to the value ∆ GMB /ε F = (2/e) 7/3 = 0.49 obtained by a naive extrapolation of the Gorkov Melik-Barkudarov result to k F a = ∞.
At T c our results are in very good agreement with those of Burovski et al. except for the value of the dimensionless chemical potential µ/ε F and that of the entropy per particle at T c . Now Burovski et al. have obtained their values for the pressure p/nε F and the entropy S/N indirectly from the internal energy and the chemical potential by using 3pV = 2U right at unitarity and the Gibbs-Duhem relation. The different results for the chemical potential then entail the considerable discrepancy in the value of S/N at T c . Within our numerical scheme, the chemical potential is the most directly -via (2.26) -obtainable quantity among the thermodynamic data. In light of the excellent agreement of all other quantities with the numerical results of Burovski et al., the discrepancy for the chemical potential is thus quite surprising. Indeed, we believe that our values for both the chemical potential and the entropy, for which the validity of the Gibbs-Duhem relation and of 3pV = 2U at unitarity have been checked independently, are rather close to the exact results. This point of view is supported by considering the evolution of the entropy per particle right at the critical temperature as a function of the dimensionless coupling. In the BCS limit, the entropy associated with single particle excitations can be calculated from the exactly soluble reduced BCS-Hamiltonian and is given by the standard mean-field expression (2.50). At the critical temperature this entropy coincides with that of an ideal Fermi gas (3.13) Since the ratio T c /T F is exponentially small in the weak coupling limit, the entropy (3.13) associated with fermionic excitations is dominant compared to the contribution arising from the collective Bogoliubov-Anderson mode. Indeed, extrapolating the corresponding low temperature entropy (3.5) associated with collective excitations up to the critical temperature gives rise to a contribution of order (T c /T F ) 3 , which is negligible compared to (3.13).
At very large coupling strengths, the strongly bound fermion pairs form an eventually ideal Bose gas, for which the entropy per particle right at T c can again be determined analytically. Recalling, that the number of bosons N B = N/2 in this limit is just half the number of fermions, we obtain a universal number S(T c )/N = 5ζ(5/2) 4ζ(3/2) = 0.6417... . (3.14) As is evident from Fig. 14, where the complete evolution of the ratio S(T c )/N is shown as a function of the dimensionless coupling parameter v, the limiting value of the ideal Bose gas is in fact not far from the entropy which is obtained from the Shohno-Popov theory of noninteracting bosonic quasiparticles in the range v > 1, according to (2.83). It is interesting to note, that the entropy per particle right at T c exhibits a maximum as a function of the coupling constant of order S(T c )/N ≈ 0.78 around the same coupling, where the critical temperature exhibits a maximum. Considering the smooth evolution of S(T c )/N as a function of v, the value S(T c )/N ≈ 0.16 at unitarity, which is deduced from the results of Burovski et al., appears to be far too small. On the other hand, the result S(T c )/N ≈ 0.99 obtained by Bulgac et al. seems to be too high.

IV. DISCUSSION AND CONCLUSION
In conclusion let us summarize what has been achieved, mention shortcomings of the present approach and indicate possible future extensions.
The formal basis of our results is a self-consistent, conserving theory, which is based on an approach due to Luttinger-Ward and DeDominicis-Martin, in which the exact one-or two-particle Green functions serve as an infinite set of variational parameters. In order for this approach to provide consistent thermodynamic results it is essential that the Green functions satisfy self-consistency conditions which reflect the stationarity of the appropriate thermodynamic potentials. Approximate formulations, in which free Green functions are replaced by full ones according to a choice of G 0 G 0 , GG 0 or GG will in general not obey conservation laws or exact thermodynamic identities, in contrast to the Φ derivable formulation presented here. The stationarity conditions were also crucial for the proof of thermodynamic equivalence of the Luttinger-Ward with the DeDominicis-Martin formalism on the level of our approximate functional for the grand canonical potential or the entropy. In fact, to our knowledge, the theory presented here is the first concrete application of the DeDominicis-Martin formulation to the fermionic many-body problem.
An important point, we want to emphasize, is the necessarily self-consistent nature of the formalism. Indeed, within the Luttinger-Ward or the DeDominicis-Martin formulation an approximate functional for the grand canonical potential Ω[G] or the entropy S[G, Γ] is made stationary by determining the space-and time-dependent Green and vertex functions from the variational conditions (2.9) and (2.13) respectively. The solution of these equations necessarily leads to a self-consistent mutual dependence of the various Green functions. Self consistency is thus reached precisely at the stationary point of these functionals. At this point, equations (2.9) and (2.13) are valid, conditions which are necessary for the theory to give consistent thermodynamics, as pointed out e.g. in the context of equation (2.27).
A well known shortcoming of conserving approximations is the dichotomy with the gapless nature of the collective modes, which reflects the broken continuous symmetry of the superfluid state. For the present theory the formal reason for this dichotomy is a violation of the Ward identity resulting from the global gauge symmetry of the exact theory. In order to overcome this problem an extension of the theory was introduced which forces the gapless nature in the symmetry broken phase while remaining Φ-derivable at the same time to maintain the conserving property.
We have provided quantitative results for essentially all thermodynamic properties at temperatures below half the Fermi temperature, thus covering the relevant regime of the degenerate gas. Overall our results agree remarkably well with recent numerical calculations at the unitarity point giving confidence that our approach is quantitatively reliable over the full range of couplings between the BCS-and the BEC-limit. In particular, we provide concrete predictions for a number of universal ratios characterizing the unitary Fermi gas both at T = 0 and at T = T c .
The extensive numerical work entering the solution of the stationarity constraints and thermodynamic potentials is reflected most clearly in the three dimensional plots of the entropy Fig. 6, pressure Fig. 7 and the order parameter Fig. 8. Most noteworthy are the quite abrupt change from fermionic to bosonic character for v in the interval −1 < v < +1 which are most obvious in the entropy and pressure and the quick convergence to a Shohno-Popov theory of noninteracting bosonic quasiparticles for v > 1.
An initially unexpected result, which is clearly visible in the numerical data, is the fact that our superfluid phase transition is weakly first order, instead of being continuous as it should be. The origin of this failure to capture the critical behaviour correctly is found in the Shohno-Popov theory, which is obtained from our approach in the limit v ≫ 1. The Shohno-Popov theory of a dilute, repulsive Bose gas generalizes the Bogoliubov theory to finite temperatures. It takes into account the thermal depletion of the condensate by including the effect of bosons with finite momentum K in the particle number equation. Long ago Reatto and Straley [44] analyzed Shohno's theory in a self-consistent formulation and obtained a first-order superfluid transition. Physically, the origin of the associated entropy jump is the collapse of the single particle spectrum right at the transition. Indeed, within the Shohno-Popov theory, the single particle spectrum changes from initially linear to initially quadratic on raising the temperature through T c . As a result, the density of states is changed from a ε 2 dependence below T c to the free particle √ ε result right at and above T c . The associated drastic increase in the available phasespace leads to a jump in the entropy.
For a purely bosonic system, a proper treatment of the behaviour near the critical point was recently given by Baym and coworkers [51]. Baym and Holzmann [68] showed that a change of the spectrum for long wavelength excitations occurs right at T c . This hardening of the spectrum (the low K spectrum is of the form K α with α < 2) leads to the required reduction in the density of states to render the superfluid transition continuous. The subtle low K correlations necessary for this change in spectrum are clearly missing in our self-consistent approach.
The BCS-BEC crossover being continuous however implies that the first order result also pertains to the v ≪ −1 limit of our theory. We have checked that at the transition the discontinuities of all thermodynamic quantities are ∼ exp(−C|v|) for v ≪ −1 [47]. The associated difficulties of a proper treatment of bosonic excitations do not occur in the reduced BCS hamiltonian which neglects collective modes altogether, resulting in a continuous superfluid transition. To correctly account for the critical regime ∆T /T c → 0 our theory would need to be extended to treat the feedback between different bosonic modes accurately. Bickers and Scalapino [69] have shown that this requires the incorporation of single particle self consistency and two particle self consistency on the same level of approximation. This may be achieved via so called parquet resummations. Currently, however, the inclusion of these contributions appears extremely challenging. A systematic and analytically accessible description of the crossover which is uniformly valid in both the normal and superfluid regime and which gives a proper account of the critical behaviour is provided by a 1/N -expansion as recently shown by Nikolić and Sachdev [15]. This method can in fact be extended in a straigthforward manner to the case of unbalanced spin populations, a subject which has attracted a lot of attention very recently [64,70].
In our formulas of the thermodynamic potentials most sums over Matsubara frequencies are not well defined. The functions which are summed do not decay to zero fast enough so that the Matsubara-frequency sums diverge. However, this problem can be fixed. To do this we first perform a Fourier back transformation to obtain a function in terms of the imaginary time τ . Then we take the limit τ → −0 or τ → +0 which is finite and well defined.
We must distinguish between fermion and boson functions which have different Matsubara frequencies.
Fermion functions are of the type A(k, ω n ) = A(k, ω n ) B(k, ω n ) −B(k, ω n ) * A(k, ω n ) * (A1) where A(k, ω n ) may be either A(k, ω n ) = − ln[G(k, ω n )] or A(k, ω n ) = [G 0 (k, ω n ) −1 G(k, ω n ) − 1]. (Note that the lower row of the matrix (A1) has the opposite sign than the lower row of the matrix (2.17). The reason ist that in the terms of the thermodynamic potentials always an even number of fermion Green functions is multiplied together.) In this case we define where we assume that A(K, τ ) is real.
In some terms of our formulas the fermion function A(k, ω n ) or the boson function A(K, Ω n ) is proportional to the unit matrix 1. In this case the Fourier backtransform A(K, τ ) is δ F (τ / ) or δ B (τ / ), respectively. Hence, the related Matsubara-frequency sums (A2) or (A4) are zero.

APPENDIX B: NUMERICAL FOURIER TRANSFORMATION
The special numerical Fourier transformation has been invented long time ago by one of the authors in a different context in order to solve the mode-coupling equation for the liquid-glass transition [71]. In this case relaxation phenomena are considered on a logarithmic time scale over many decades, starting at microscopically short times and extending up to very long macroscopic times. Thus, a Fourier transformation is needed which can handle functions with features on logarithmic time and frequency scales extending over ten and more decades. Clearly, a standard fast Fourier transformation can not be applied because a constant step width would be needed. Rather the function to be transformed has been discretized on a logarithmic scale and interpolated by cubic spline polynomials. Since for polynomial functions the Fourier integrals can be evaluated exactly, we end up with a transformation formula which depends on the spline coefficients of the function.
Later this special numerical Fourier transformation has been extended to transform Matsubara Green functions in order to solve the self-consistent equations for the BCS-BEC crossover [30]. Here, three-dimensional spatial Fourier transformations of isotropic functions and discrete Fourier sums with Matsubara frequancies were considered. These Fourier transformations are used also in the present paper for the numerical calculations. Only a few modifications and optimizations have been made over the years. The basic principles of the special numerical Fourier transformation are descibed in the appendix of Ref. 30. Here we present the fundamental formulas in order to make the numerical method available for applications.
In order to perform a discrete Fourier transformation the following sum must be evaluated where x is a discrete variable with constant step width ∆x. In this formula and in the following formulas the sum over x is defined as a trapezoid sum. This means that the first term and the last term in the sum are multiplied by a factor 1 2 , respectively. The continuous Fourier transformation is defined by a related integral which is obtained from (B1) in the limit ∆x → 0.
We assume that the function values are known in a finite subset of points x j according to f (x j ) = a j where j = 0, 1, . . . , N . The points x j cover the whole interval between x min and x max on a logarithmic scale so that x min = x 0 < x 1 < . . . < x N −1 < x N = x max . Consequently, the Fourier sum (B1) can be divided into a sum of N trapezoid sums according to Now, we assume that the function is given by the cubic spline polynomial if x is located in the interval x j ≤ x ≤ x j+1 . The spline coefficients a j , b j , c j , and d j are calculated numerically. Inserting the cubic spline polynomial (B3) into the formula (B2) we find that the trapezoid sums within the curved brackets can be evaluated exactly. Thus, as a result we obtain the Fourier transform (B5) By construction a cubic spline function and its first two derivatives are continuous. These facts imply the following continuity conditions which may be used to regroup the terms in (B4). Consequently, as a result we obtain the alternative formula f (k) =J (0) (k) e ikxN a N − e ikx0 a 0 The terms with spline coefficients a j , b j , and c j have cancelled for j = 1, 2, . . . , N −1. In the limit k → 0 the functions (B10) diverge according to J (n) (k) ∼ |k| −(n+1) . For this reason, the alternative formula (B9) can be applied numerically only for large k so that |kx j | 1 for all j = 0, 1, . . . , N . On the other hand, the functions (B5) are finite in the limit k → 0 so that the formula (B4) can be applied numerically for small k where |kx j | 1 for all j = 0, 1, . . . , N . In practice we use a combination of both formulas (B4) and (B9). Which formula is used for a particular j we decide by considering the value of |kx j | and comparing this value with 1. In this way we obtain a special numerical Fourier transformation which is stable and reliable for points x j and k l distributed on a logarithmic scale over many decades.
We have derived our special numerical Fourier transformation for discrete variables x with a finite constant step width ∆x. The continuous Fourier transformation is obtained easily and naturally by taking the limit ∆x → 0 in the functions (B5) and (B10) which is well defined.
In order to transform the Green and vertex functions forward and backward, we need two kinds of Fourier transformations. First we transform between the Matsubara frequencies and the imaginary time variable. In this case we can apply a continuous (forward) and a discrete (backward) Fourier transformation (B1) directly. Secondly we transfrom between the wave vector and the spatial coordinate in d = 3 dimensions. Since the functions are spherically symmetric, an integration over the angles can be performed, so that the resulting transformation integrals are one dimensional depending only on radial variables, a radial wave number and a radial space coordinate, respectively. For d = 3 the transformation integrals can be recast into a one-dimensional continuous Fourier transformation so that our special numerical Fourier transformation (B1) can be used once again.
In practice we use N = 300 points for all variables. The values of the wave numbers and the values of the space coordinates are distributed on logarithmic scales over six decades, respectively. The Matsubara frequencies are distributed on a logarithmic scale over about twelve decades. The imaginary time variables are distributed appropriately over a finite interval with two logarithmic scales, one for each boundary.
The Green and vertex functions are singular and have slowly decaying long tails. For this reason, reference functions must be subtracted which remove the singularities and the long tails. The reference functions are derived from free Green functions and the two-particle scattering amplitude (T matrix). For these reference functions analytical expressions must be available in all Fourier representations. The difference functions f (x) which are eventually transformed by our numerical method (B1) must be smooth in x and decay accoording to f (x) ∼ x −2 or faster for |x| → ∞.