Correlations in the low-density Fermi gas: Fermi-Liquid state, Dimerization, and BCS Pairing

We present ground state calculations for low-density Fermi gases described by two model interactions, an attractive square-well potential and a Lennard-Jones potential, of varying strength. We use the optimized Fermi-Hypernetted Chain integral equation method which has been proved to provide, in the density regimes of interest here, an accuracy better than one percent. We first examine the low-density expansion of the energy and compare with the exact answer by Huang and Yang (H. Huang and C. N. Yang, {\em Phys. Rev.\/} {\bf 105}, 767 (1957)). It is shown that a locally correlated wave function of the Jastrow-Feenberg type does not recover the quadratic term in the expansion of the energy in powers of $\a0\KF$, where $\a0$ is the vacuum $s$-wave scattering length and $\KF$ the Fermi wave number. The problem is cured by adding second-order perturbation corrections in a correlated basis. Going to higher densities and/or more strongly coupled systems, we encounter an instability of the normal state of the system which is characterized by a divergence of the {\em in-medium\/} scattering length. We interpret this divergence as a phonon-exchange driven dimerization of the system, similar to what one has at zero density when the vacuum scattering length $\a0$ diverges. We then study, in the stable regime, the superfluid gap and its dependence on the density and the interaction strength. We identify two different corrections to low-density expansions: One is medium corrections to the pairing interaction, and the other one finite-range corrections. We show that the most important finite-range corrections are a direct manifestation of the many-body nature of the system.


I. INTRODUCTION
The study of ultracold quantum gases is interlinked with controlling magnetic Fano-Feshbach resonances and thereby changing the effective interparticle interaction by many orders of magnitudes [1,2], see also Refs. [3][4][5]. This makes ultracold Fermi gases a convenient tool to study the behavior of a degenerate fermionic many-body system [6] over a wide range of interaction strengths, in particular fermionic superfluidity [7]. Indeed, changing the magnetic field across a resonance makes it possible to continuously tune the gas from the Bardeen-Cooper-Schrieffer (BCS) state of Cooper pairs to a Bose-Einstein condensate (BEC) of weakly bound molecules. After the first observations of a molecular BEC of fermions [8][9][10], this so-called BCS-BEC cross-over has been widely studied experimentally [11][12][13][14][15][16][17] and theoretically [18][19][20][21], see also reviews Refs. 22 and 23. The observation of quantized vortices on both sides of the BCS-BEC crossover provided an unambiguous proof of superfluidity by fermionic pairing [24]. Recent work investigated the effect of partial polarization of a two-component Fermi gas on the Fermi liquid parameters [25], the nature of the transition from a BCS state to a state of a molecular BEC [26], and quantifying the superfluid fraction in a Fermi gas by means of second sound [27].
In this work we concern ourselves with the low-density properties of homogeneous Fermi gases at zero temperature. We will use a square-well and a Lennard-Jones interaction potential for our study. Changing the interaction strength (coupling constant) of the respective potential changes the scattering length for two-body scattering a 0 , which we refer to as vacuum scattering length, as op-posed to the in-medium scattering length a introduced later. When a 0 < 0, i.e. the interaction is effectively attractive, one expects BCS type pairing of particles with opposite spin. As a 0 → −∞, a low energy resonance of the two-body problem generates bound dimers. This unitary limit, where the only relevant length scale is the inverse Fermi wave number k −1 F , marks the border between the BCS and BEC regime. Increasing the interaction strength further, the Cooper pairs become weakly bound bosonic molecules. The singularity of the vacuum scattering length signals the transition.
We are in particular interested in structural quantities such as the energetics, distribution functions, the stability of the system against spinodal decomposition and dimerization, and possibly BCS pairing. We utilize a quantitative method of microscopic many-body theory to determine correlation effects, i.e. effects beyond the weak coupling or mean-field approximations [28] that are routinely applied at low densities.
For example, the ground state energy depends in the low-density limit only on the dimensionless parameter k F a 0 [28,29]. We are interested here in the question in which parameter range this "universal" behavior persists. Another example for correlation effects is the pair distribution function, g(r). In correlation-free mean-field approaches, g(r) is equal to the distribution function of the non-interacting Fermi gas, g F (r). As we will see below, g(r) deviates, in particular for spin-antiparallel particles, substantially from g F (r) when the absolute value of the scattering length becomes large.
In the limit of weak interactions, the system can be described by a BCS type wave function. For systems where the weak-coupling approximation does not ap-ply (for example for Lennard-Jones type interactions), the pairing gap can be obtained by an extension of the Jastrow-Feenberg variational method, the correlated BCS (CBCS) method [30,31], reviewed in section III A. See also Refs. 32 and 33 and 34 for a similar implementation of the same ideas. The CBCS theory takes into account short-ranged correlations analogously to the theory for normal systems, and supplements these by the typical BCS correlations. In its essence, CBCS provides a recipe for calculating an effective interaction that enters the standard BCS formalism. An alternative way to deal with the problem that implements a full Fermi Hypernetted Chain (FHNC) summation for large gap parameters has been suggested by Fantoni [35,36]; unfortunately, the approach uses a normalization of the correlated BCS wave function which leads to divergences for optimized or otherwise long-ranged correlations. It was therefore replaced in Ref. 30 by the method reviewed in Sec. III A.
Our paper is organized as follows: In Sec. II we will review briefly review the basics of the correlated basis functions (CBF) method. We call this approach "generic" many-body theory because the same equations can be derived from Green's functions approaches [37], from Coupled Cluster theory [38] and from an extension of density functional theory which includes pair correlations. We evaluate in section II B the low-density limit and show that the exact formula [29] is not reproduced by the Jastrow-Feenberg and/or the "fixed node" approximation. To correct this problem, we review in section II C perturbation theory in a correlated basis and show that second-order CBF corrections must be added.
In Sec. III A we will review the CBCS theory. It is seen that the theory can be formulated in exactly the same way as ordinary BCS theory, CBCS simply provides a means for deriving weak, effective interactions from a strong, bare interaction. Upon closer inspection, the mapping of the bare interaction to an effective interaction is closely related to the transition from the bare interaction to the T -matrix used in the low-density expansion of BCS theory [39,40].
In Sec. IV, we present our results for the energy, the pair distribution function, the in-medium scattering length, and the gap energy as function of Fermi wave number k F and the vacuum scattering length a 0 . We show that the dynamical correlations can be characterized by three regimes: For short distances, r ≈ σ, these correlations are, of course, determined by the interaction. An intermediate regime is dominated by two-body scattering where correlations decays as 1/r. A third, asymptotic regime for r larger than the average interparticle distance is dominated by many-body effects where the correlations decay as 1/r 2 .
We find an instability of the FHNC-Euler-Lagrenge (FHNC-EL) solutions accompanied by a divergence of the in-medium scattering length a, which we interpret as onset of dimerization. This happens well before the divergence of the vacuum scattering length a 0 and is caused by the induced interaction mediated by phonon exchange.
Thus, dimers can be formed at finite density even if the bare potential does not have a bound state. Wave functions of the type normally used only in the "BEC" regime [41][42][43] would be more appropriate in this regime.
In the regime a < 0 we solve the CBCS gap equation and show that deviations from the simple BCS approximation can be separated into contributions from the induced interaction, i. e. the density dependence of the in-medium scattering length, and from the non-negligible momentum dependence of the pairing interaction in the CBCS gap equation.

A. Variational wave functions
We start our discussion with the Jastrow-Feenberg theory for a strongly interacting, translationally invariant normal system. As usual, we assume a non-relativistic many-body Hamiltonian in terms of a local, phenomenological interaction v(r).
The method starts with a variational ansatz for the wave function [44] Ψ 0 (r 1 , . . . , where Φ 0 (r 1 , . . . , r N ) is a model state, normally a Slater-determinant for fermions and Φ 0 (r 1 , . . . , r N ) = 1 for bosons. There are basically two ways to deal with this type of wave function. In quantum Monte Carlo studies, the wave function (2.2) is referred to as "fixed node approximation" and an optimal correlation function F (r 1 , . . . , r N ) is obtained by stochastic means [41][42][43][45][46][47][48][49][50][51][52], decomposition into n-body correlations u n (r 1 , . . . , r n ) is then, of course, not necessary. Alternatively, one can use diagrammatic methods, specifically the optimized Fermi-hypernetted chain (FHNC-EL) method for the calculation of physically interesting quantities. These diagrammatic methods have been successfully applied to such highly correlated Fermi systems as 4 He and 3 He at T = 0 [53], they are naturally expected to work much better in the low density systems of interest here. In fact, we have shown in recent work [54] that even the simplest version of the FHNC-EL theory is accurate within better than one percent at densities less that 25 percent of the ground state density of liquid 3 He. Diffusion Monte Carlo calculations typically use a parametrized Jastrow-Feenberg ansatz for importance sampling, where the parameters are optimized by variational Monte Carlo calculations. JF theory makes explicit use of the form (2.3). It has been shown [37] that triplet correlations contribute to the ground state energy only in fourth order of the interactions. Even in strongly interacting quantum fluids like the helium liquids, triplet correlations contribute no more than five to ten percent to the ground state energy [53,55] in both isotopes, they are completely negligible below approximately 25 percent of the respective equilibrium densities. It is therefore admissible, at the low densities we are concerned with here, to identify the Jastrow-Feenberg approximation with the fixed-node approximation in quantum Monte Carlo calculations. Monte Carlo calculations based on lattice Hamiltonians have been performed for the unitary Fermi gas, where this nodal restriction was not required [48,51].
The correlations u n (r 1 , . . . , r n ) are obtained by minimizing the energy, i.e. by solving the Euler-Lagrange (EL) equations The evaluation of the energy (2.6) for the variational wave function (2.3) and the analysis of the variational problem is carried out by cluster expansion and resummation methods. The procedure has been described at length in review articles [56] and pedagogical material [31]. In any approximate evaluation of the energy expectation value, it is important to make sure that the resulting equations are consistent with the exact variational determination of the correlations. It has turned out that the (Fermi-)hypernetted chain hierarchy ((F)HNC) of approximations is the only systematic approximation scheme that preserves the properties of the variational problem [44].
Here, we spell out the simplest version of the equations that is consistent with the variational problem ("FHNC-EL//0"). These do not provide the quantitatively best implementation [53], instead, they provide the minimal version of the FHNC-EL theory. In particular, they contain the relevant physics, namely the correct description of both short-and long-ranged correlations. They also are the minimal implementation of the theory that gives the correct expansion of the ground state energy in powers of (k F a 0 ) for the wave function Eq. (2.2).
In the FHNC-EL//0 approximation [57], which contains both the random phase approximation (RPA) and the Bethe-Goldstone equation in a "collective" approximation [31], the Euler equation (2.6) can be written in the form where S(k) is the static structure factor of the interacting system, t(k) = 2 k 2 /2m is the kinetic energy of a free particle, S F (k) is the static structure of the noninteracting Fermi system, and is what we call the "particle-hole interaction". For further reference, we have above defined the "Clark-Westhaus effective interaction" v CW (r) [56]. As usual, we define the Fourier transform with a density factor, Auxiliary quantities are the "induced interaction" and the "direct-direct correlation function" Eqs. (2.7)-(2.11) form a closed set which can be solved by iteration. Note that the Jastrow correlation function (2.3) has been eliminated entirely.
The pair distribution function can generally be written as where, roughly speaking, Γ dd (r) describes dynamic, short-range correlations, g F (r) = 1 − 1 ν ℓ 2 (rk F ) is the pair distribution function of the non-interacting Fermi gas, and ℓ(x) = 3 x j 1 (x) is the Slater exchange function. C(r) describes the combination of statistical and dynamic correlations. In leading order in the dynamic correlations we haveC where (∆X ee ) 1 (k) is represented by the two leading order exchange diagrams shown in Fig. 1. The energy per particle is, in this approximation [53,54], JF , 2m is the Fermi energy of non-interacting particles, ν is the degree of degeneracy of the single particle states; in our case we have generally ν = 2. The term t (3) JF is the three-body term of the Jackson-Feenberg kinetic energy. The function Γ dcc (r 1 ; r 2 , r 3 ) is the sum of all three-point diagrams that have an exchange path connecting points r 2 and r 3 and no exchange lines attached to point r 1 which is dynamically connected in such a way that there exists a path between r 1 and each of the other two external points that does not go trough the third external point. The term is normally numerically very small, we keep it here for the purpose of deriving the low-density expansion. To obtain the correct low-density limit, we retain all contributions to t JF with two factors Γ dd (r): The term t (3b) JF cancels exactly the contribution to e Q originating from the second diagram in (∆X ee ) 1 (r); the terms e The term t CW is recognized as the three-body term of the "Clark-Westhaus" form of the kinetic energy. To summarize, the total energy has the form For further reference, we also spell out the pair distribution functions in the spin-parallel and the spinantiparallel channel: where [. . .] F (r) indicates the Fourier-transform (2.9). To leading order in the density, the term ℓ 2 (rk F ) is the only term that reflects Fermi statistics, whereas the factor [1 + Γ dd (r)] describes dynamical correlations. g ↑↑ (r), g ↑↓ (r), and g(r) are normalized such that both go to unity for large r.

B. Low-density limit
In the limit of low densities, the equation of state and related quantities depend only on the vacuum swave scattering length a 0 and the Fermi wave number k F . For example, the energy per particle has the expansion [28,29] (2.23) Note that the expansion (2.23) is strictly valid only for a 0 > 0, for attractive potentials the superfluid condensation energy must be added.
The locally correlated wave function (2.2) is not exact, and the question arises whether it recovers the expansion (2.23). It is plausible that this is not the case: The calculation of the third term in the expansion (2.23) makes explicit use of the form of the energy denominator in second order perturbation theory [28]. The local correlation operator corresponds to a "collective approximation" in which, among others, the particle-hole propagator is approximated by a collective mode.
Our task is to express the variational energy expression (2.19) to second order in the vacuum scattering length a 0 . One can deal with this task in two ways: One is to permit hard-core interactions, the other, somewhat simpler, approach is to assume a weak interaction that has a Fourier transform. In this case, one can parallel the derivation of Ref. 37 for fermions.
We will show the details of the calculation in Appendix A 1, here we discuss only the essential steps: The vacuum scattering length is determined from the zero-energy scattering equation The scattering equation has the asymptotic solution Multiplying Eq. (2.24) with ψ(r) and using the identity gives a relationship that will be useful later: We notice that the induced interactionw I (k) as defined in Eq. (2.10) is of second order in the interaction. To leading order in the density we can also expand Eq. (2.7) and obtain from Eq. (2.11) the solutioñ In addition to calculating the energy contributions (2.19) for the correlation function (2.29), we must expressṽ CW (0+) in terms of the scattering length becausẽ v CW (0+) is calculated with the optimal correlation function (2.11) of the many-body problem at finite density, and not with the solution of the zero-density scattering equation (2.24). In Appendix A 1 we will prove the relationshipṽ Collecting all results, one finds (2.31) see Appendix A 1 for details of the calculation. The numerical factor 2430284/1576575 = 1.5415 is to be compared with the factor 4(11 − 2 ln 2)/35 = 1.098 of Eq. (2.23). We emphasize again that the result (2.31) also applies to the "fixed-node" approximation in Monte Carlo calculations because the terms where that approximation deviates from our expansion are of at least fourth order in the potential strength. To get the exact result, one must go beyond local correlation operators; this is done by perturbation theory in a correlated basis generated by the correlation operator F (r 1 , . . . , r N ) described in the next section II C. The situation is analogous to the case of the high-density limit of the correlation energy of the electron gas. With local correlations one obtains for the logarithmic term 0.05690 ln r s Ry [59] instead of the exact value 0.06218 ln r s Ry [60,61]. This deficiency is, for the electron gas, removed by second-order CBF theory [62], and we will see that the same is true here.

C. Elements of Correlated Basis Functions
We have seen above that a locally correlated wave function (2.2) does not produce the exact low-density limit (2.23) of the ground state energy. As mentioned above, the problem can be cured by applying second-order perturbation theory with correlated basis functions (CBF theory). We will also need the basic ingredients of CBF theory for examining the superfluid state. We review here this method only very briefly, details may be found in pedagogical material [31] and review articles [53,56]; the diagrammatic construction of the relevant ingredients is given in Ref. 63.
CBF theory uses the correlation operator F to generate a complete set of correlated and normalized N -particle basis states through are not orthogonal, perturbation theory can be formulated in terms of these states [44,64].
For economy of notation, we introduce a "secondquantized" formulation of the correlated states. The Jastrow-Feenberg correlation operator in (2.3) depends on the particle number, i.e. F = F N (1, . . . , N ) (whenever unambiguous, we omit the corresponding subscript). Starting from the conventional a † k , a k that create and annihilate single particle states, new creation and annihilation operators α † k , α k of correlated states are defined by their action on the correlated basis states: According to these definitions, α † k and α k obey the same (anti-) commutation rules as the creation and annihilation operators a † k and a k of uncorrelated states, but they are not Hermitian conjugates. If Ψ m is an Nparticle state, then the state in Eq. (2.33) must carry an (N +1)-particle correlation operator F N +1 , while that in Eq. (2.34) must be formed with an (N − 1)-particle correlation operator F N −1 .
In general, we label "hole" states, which are occupied in |Φ o , by h, h ′ , h i , . . . , and unoccupied "particle" states by p, p ′ , p i , etc.. To display the particle-hole pairs explicitly, we will alternatively to Ψ m use the notation Ψ p1...p d h1...h d . A basis state with d particle-hole pairs is then For the off-diagonal elements O m,n of an operator O we sort the quantum numbers m i and n i such that |Ψ m is mapped onto |Ψ n by From this we recognize that, to leading order in N , any O m,n depends only on the difference between the states |Ψ m and |Ψ n and not on the states as a whole. Consequently, O m,n can be written as matrix element of a d-body operator (2.37) (The index a indicates antisymmetrization.) Key quantities for the execution of the theory are diagonal and off-diagonal matrix elements of unity and [31,63] of the matrix elements of H ′ m,n into the off-diagonal quantities W m,n and N m,n and diagonal quantities H m,m .
To leading order in the particle number, the diagonal matrix elements of H ′ ≡ H − H o,o become additive, so that for the above d-pair state we can define the CBF single particle energies is an average field that can be expressed in terms of the compound diagrammatic quantities of FHNC theory. According to (2.37), W m,n and N m,n define d−particle operators N and W, e.g.
Diagrammatic representations of N (1, 2, . . . , d) and W(1, 2, . . . , d) have the same topology [63]. In homogeneous systems, the continuous parts of the p i , h i are wave numbers p i , h i ; we abbreviate their difference as q i . In principle, the N (1, 2, . . . , d) and W(1, 2, . . . , d) are non-local d-body operators. In the next section, we will show that we need, for examining pairing phenomena, only the two-body operators. Moreover, the low density of the systems we are examining permits the same simplifications of the FHNC theory that we have spelled out in Sec. II A. In the same approximation, the operators N (1, 2) and W(1, 2) are local, and we have [53] N (1, 2) = N (r 12 ) = Γ dd (r 12 ) The most straightforward application of CBF theory is to calculate corrections to the ground state energy. In second order we have, for example The magnitude of the CBF correction is normally comparable to the correction from three-body correlations [53]; it is also important to note that there are significant cancellations between the two terms in the numerator. We will show in appendix A 1 that the CBF correction (2.44) corrects the coefficient of the third term in the expansion (2.31) and leads to the exact low-density limit (2.23).

A. General derivation
We show in this section how the variational theory is adapted to the description of the superfluid state. We restrict ourselves here to the simplest case of S-wave pairing and show how the effective interactions, which enter phenomenological theories as parameters, may be calculated from first principles. This section reviews the derivations of Refs. 30 and 31.
The BCS theory of fermion superfluidity generalizes the Hartree-Fock model {|Φ m } by introducing a superposition of independent particle wave functions corresponding to different particle numbers The coefficient functions u k and v k , known as Bogoljubov amplitudes, describe the distortion of the Fermi surface due to the pairing phenomenon.
To deal with strongly interacting systems, adequate provision must be made for the singular or near-singular nature of the two-body interaction v(r) for small r. To build the required geometrical correlations into the microscopic description of the system, we can define a correlated BCS state, incorporating both short-ranged and BCS correlations. We are faced with a formal mismatch, which prevents us from simply applying the correlation factor F N to BCS , since the former is defined in the N -particle Hilbert space and the latter is a vector in Fock space with indefinite particle number. The most natural way to deal with this is first projecting the bare BCS state on an arbitrary member of a complete set of independent-particle states with fixed particle numbers, then applying the correlation operator to that state, normalizing the result, and finally summing over all particle numbers. We must then distinguish between correlation operators and normalization integrals corresponding to different particle numbers N . Thus, the correlated BCS state is written as The trial state ( To derive the relevant equations we consider the expectation value of an arbitrary operatorÔ with respect to the superfluid state: One may pursue cluster-expansion and resummation methods of expectation values (3.3) for the superfluid trial state (3.2), this has been done successfully for the one-and two-body density matrices corresponding to a slightly different choice of the correlated BCS state [35] which exhibits, unfortunately, divergences for optimized correlation functions. We do not follow this route, but instead consider the interaction of only one Cooper pair at a time. The error introduced by this is of order ξ = (∆ F /ǫ F ) 2 , where ∆ F is the superfluid gap energy. We will demonstrate below that this quantity is indeed small in the regime where the wave function (3.2) is appropriate.
To implement the decoupling approximation, it is sufficient to retain the terms of first order in the deviation v 2 k − v 2 0,k and those of second order in u k v k . The calculation of Ĥ − µN for correlated states [30] is somewhat tedious, we only give the essential steps and the final result. It is convenient to introduce the creators and annihilators of correlated Cooper pairs, In terms of these quantities, the expectation value of an operatorÔ is, to leading order in the amplitudes v 2 In Eqs. oo the expectation value of the operatorÔ in the N -particle correlated ground state. Insertinĝ H − µN forÔ into the expansion (3.5), where µ is a Lagrange multiplier (the chemical potential) introduced to adjust the average particle number N s = N , we recover the effective interactions, overlap integrals (2.42), and single-particle energies (2.41) of section II C, e.g.
(3.8) Accordingly, we may write the energy of the superfluid state in the form with the "pairing interaction" With the results (3.9) and (3.10), we have arrived at a formulation of the theory which is formally identical to the BCS theory for weakly interacting systems. Upon closer inspection (see the next section) we will see that our formulation corresponds to a BCS theory formulated in terms of the scattering matrix [65]. The correlation operator serves here to tame the short-range dynamical correlations, the effective interaction W(1, 2) is just an energy independent approximation of the T -matrix.
We may proceed now in the conventional way to determine the Bogoljubov-amplitudes u k , v k , by variation of the condensation energy (3.9) to compute the superfluid condensation energy, or to investigate the local stability of the normal ground state by second variation. Minimization of the energy expectation value determines the BCS amplitudes u k , v k . The CBCS gap equation becomes The conventional, i.e. uncorrelated, BCS gap equation [66] is retrieved by replacing the effective interaction P kk ′ by the matrix elements of the bare interaction. Note that our "decoupling approximation" simply means that we assume that the pairing interaction P kk ′ does not depend on the Bogoljubov amplitudes. This does not assume that the gap is small compared with the Fermi energy.

B. Analysis of the gap equation
In the local approximations appropriate for low densities, the effective interaction is given by Eqs. (2.43). The pairing matrix element is expressed in terms of the Fourier-transformsW(k) andÑ (k): The remaining arguments are standard, cf. Ref. 65 and 66: If the gap at the Fermi surface is small, we can replace the pairing interactionW(k) by its s-wave matrix element at the Fermi surface, Then we can write the gap equation as which is almost identical to Eq. (16.91) in Ref. 65. In particular, the second term has the only function to regularize the integral for large k ′ . We can, therefore, immediately conclude that the zero temperature gap is, in this approximation, given by The low-density limit is then obtained by identifying a with the vacuum scattering length a 0 : Of course, our equations (3.10), (3.11) are much more general: At low densities, the subtraction termi.e. the second term in the square bracket of Eqs. (3.12) and (3.14) is important to regularize the integral for large momentum transfers. At higher densities, where the finite range of the interaction provides that momentum cutoff, the same term is negligible since the energy numerator term (|e k − µ| + |e k ′ − µ|) is zero at the Fermi momentum.
By comparison with the low-density limit and Eq. (2.28) we interpret the constant a ≡ m 4πρ 2W (0+) (3.18) as an "in-medium" scattering length. Hence, at finite densities, one expects two types of corrections: (i) Medium corrections: The effective pairing interac-tionW(k) is related toṽ CW (k) through Eqs. (2.7)-(2.11) which leads to Because of Eq. (2.30) and the fact that the induced interaction w I (r) is of second order in the interaction, we conclude that that In the same order, non-local contributions to the pairing interaction (2.43) [63] contribute. These can be identified with particle-hole ladder diagrams and vertex corrections. Topologically, one of these diagrams corresponds to the polarization correction identified by Gorkov et al. [67]. Moreover, similar to our analysis of the low-density limit, local correlation functions will not get the right coefficient of the term proportional to (a 0 k F ) 2 , hence CBF corrections to the pairing interaction [30] will also lead to modifications of order (a 0 k F ) 2 .
(ii) The solution of the s-wave gap equation is dominated by the matrix element (3.13 of W(r) at the Fermi surface which leads to the solution (3.15). Only ifW(k) is practically constant for 0 ≤ k ≤ 2k F , we can idenitfy a F with the in-medium scattering length a. The dominant finite-range correction to the pairing interaction comes from the kinetic energy term . For distances much larger than σ but smaller than the average particle distance, this term is dominated by the vacuum solution of the Euler equation, 1 + Γ dd (r) = 1 − a0 r . For distances larger than 1/k F , we obtain from Eqs. (2.7) and (2.11) that Γ dd (r) falls off like (3.21) Consequently, the effective interaction W(k) is quadratic in k for k ≤ k F and has a linear depen-denceW for k > k F . The variation of the pairing interaction between k = 0 and k = k F is interaction dependent and causes, as we shall see, a correction that is larger than the ones due to the effects mentioned above.
All but one of the correction discussed above are of order a 0 k F and higher i.e.
where α is a numerical constant. Among others, the polarization correction discussed by Gorkov [67] has this structure. Inserting the above expansion in (3.15) changes the pre-factor to In other words, to the extent that an expansion in powers of (a 0 k F ) is legitimate, all of these corrections just lead to a modified, but universal, pre-factor in Eq. (3.17). This does not apply to the finite-range correction of the pairing interaction in the relevant regime k ≤ k F . One would, of course, expect that this finite-range correction is of the same order of magnitude but its value depends, as we shall see, on details of the interaction. Hence, we conclude that the exponential behavior of Eq. (3.17) is universal whereas the pre-factor is not.

A. Energetics
We have examined in this paper two model potentials, namely a Lennard-Jones (LJ) potential and an attractive square well (SW) potential Both potentials are parametrized by a typical range σ and the depth of the attractive well ǫ. In both cases, we measure energies in units of 2 /2mσ 2 , and length in units of σ. Thus, the interaction strength ǫ and the density are the only free parameters. Our choice of interactions provides effective potentials designed to avoid the instability against clustering that exists for realistic alkali interactions, but otherwise be close to a realistic situation. The simplest connection to real interactions is provided by the vacuum s-wave scattering length a 0 . The procedure is legitimate in the low-density limit, many observable properties of these gases, such as the energy (2.23), depend indeed only on the s-wave scattering length [28,29]. For higher densities this "universal" behavior ceases. It is the purpose of our calculation to explore that area, and also study the model dependence, by comparing results for the LJ and SW model. To make contact with low-density expansions, as well as to determine the range of "universal behavior", we shall use the s-wave scattering length a 0 instead of the well-depth ǫ to characterize the potential. For the SW potential, a 0 is given by a 0 = 1 κ (κσ − tan κσ) where κ = mǫ/ 2 . For the LJ potential, a 0 must be obtained numerically. We show a 0 in Fig. 2 as a function of the potential well depth ǫ. The attractive SW potential has negative scattering length for an interaction strength below the first resonance, whereas the LJ potential has a positive a 0 below ǫ = 4.336, indicating an effectively repulsive interaction. Thus, when the interaction strength ǫ of the LJ potential is raised starting from zero, we find three regimes. (i) For psilon < 4.336 we have a 0 > 0 and there is no bound state; the many-body ground state is a normal Fermi gas. (ii) For 4.336 < ǫ < 11.18 we have a 0 < 0. There is still no two-body bound state; due to the effectively attractive potential, the many-body ground state in the low density limit is expected to be a BCS state. At ǫ = 11.18 the LJ potential has a resonance at zero scattering energy. (iii) For 11.18 < ǫ, a 0 becomes positive again, and the poten-tial supports one two-body bound state; the many-body ground state in the low density limit is expected to be a BEC state of a molecular Bose gas. All these states of Fermi gases have been studied extensively in experiments with ultracold alkali gases as discussed in the introduction. In a previous paper [54] we have already studied both, a 0 > 0 and a 0 , < 0 for fermions, and a 0 > 0 for bosons. In that work we have also examined more sophisticated versions of the FHNC-EL method and have concluded that these are necessary only at densities comparable to that of liquid helium, the reader is referred to that work to assess the range of densities for which the very simple version of the theory spelled out in Sec. II A is reliable. In short, the accuracy of our energy calculations is expected to be better than 1 percent below a density of ρ = 10 −2 σ, whereas the error of the simple FHNC-EL version is about 10 percent as the density increased to ρ = 0.4 σ which is close to the freezing density of 3 He.
We focus in the present work on the effect of attractive interactions (a 0 < 0 ). We have solved Eqs. (2.7)-(2.11) on a mesh of 2 18 points, with a resolution of 30 points between r = 0 and r = σ, amounting to a box size of 8732 σ. Such a huge box size is necessary to obtain a reasonable momentum space resolution at the very low densities we are considering here: Note that a Fermi wave number of 10 −3 σ −1 corresponds to a wavelength of 6000σ, hence this box size is the bare minimum of what one should take to resolve features of the order of the Fermi wave number. All our calculations are done for the range of interaction strength where there is no two-body bound state, i.e. before the first resonance of a 0 appears (indicated by vertical lines in Fig. 2).
Our equation of state for the two potential models is shown in Figs. 3 and 4. To recover the exact low-density limit (2.23), we have added the second-order CBF correction (2.44). To emphasize the interaction terms, we have subtracted the kinetic energy E kin = 3 5  For both interactions we observe that, dependent on the interaction strength, the equation of state begins to deviate strongly from a simple, smooth power law for a 0 k F > 0.02.
The most interesting feature we observe is that the FHNC-EL equations cease to have solutions at sufficiently large values of the density or of −a 0 . Such a limit is expected for sufficiently attractive interactions: The Fermi gas is, in the low density limit, stabilized by the Pauli pressure. As the density increases, the energy per particle becomes negative and the static incompressibility where c is the hydrodynamic speed of sound, goes to zero. Such an effect has already been reported by Owen [68]. mc 2 → 0 indicates a spinodal instability, where the system separates into a low and a high density phase. On the other hand, it is widely accepted that a low density two-component Fermi gas is subject to dimerization close to the unitary limit, −a 0 → ∞. In the following we will argue that dimerization indeed occurs, but well before the limit −a 0 → ∞. Instead, dimerization is accompanied by the divergence of the in-medium scattering length. Let us examine the question of stability from the point of view of the existence of solutions of the FHNC-EL equations: In general, the FHNC-EL equations cease to have solutions if the assumed wave function is unstable against small perturbations. This is most clearly seen for the case of density fluctuations: The term under the square root of Eq. (2.7) must be positive. In the limit k → 0+ this leads to the condition (4.4) where the right-most expression is the low density limit, a 0 k F ≪ 1, where F s 0 is small and the in-medium scattering length a is well approximated by the vacuum scattering length a 0 , see Eq. (4.6). The limit can be regarded as the low density limit of the particle-hole interaction, where we used that the induced interaction is of second order in the bare interaction. Hence, the system can be driven into an instability by holding the potential fixed and simply increasing the density factor in Eq. (4.5). Note that the above stability limit (4.4) is valid for the local correlation operator (2.2). In an improved calculation that does not rely on a local correlation operators but rather includes CBF corrections to all orders [53] the stability condition would read The condition is immediately recognized as the stability condition of Landau's Fermi Liquid theory, F s 0 > −1. Note that the ground state theory formulated in Eqs. (2.7)-(2.11) does not contain self-energy insertions (called "cyclic-chain" diagrams in the FHNC theory), which means that the effective mass is equal to 1. According to our findings in Ref. 54, this is an acceptable approximation at the low densities under consideration here. Since the interaction contribution to F s 0 is, at low densities, proportional to k F , an instability against density fluctuations could occur for sufficiently attractive interactions. This has the consequence that the FHNC-EL equations cease to have a solution.
We found indeed an instability of the solutions of the FHNC-EL equations, however, this instability does not appear to be caused by F s 0 → −1. In fact, in our calculations we were not able to come close to the point of spinodal instability. We show in Figs. 5 and 6 the value of F s 0 as a function of density for a family of potential strengths. For the SW potential, we have, depending on the coupling strength ǫ, been able to solve the FHNC-EL equations in the regime between ρ = 3.4×10 −11 σ −3 up to a critical density that depends on the coupling strength. At that density, F s 0 begins to drop very rapidly but it does not appear to approach the critical value of −1 for a spinodal instability, see Fig. 5. We also note that an instability signified F s 0 → −1 means a transition to a state with non-uniform density which is not what has been observed.
Whereas a system of fermions interacting with an attractive SW potential exhibits an instability as the density is increased, the LJ model interaction leads, due to its repulsive core, to a richer phase diagram that also features a high-density condensed phase -the liquid phase of 3 He. At low densities, we see the same picture as for the SW potential, namely that the liquid is stabilized by the Pauli pressure, but the FHNC-EL equations cease to have solutions above a certain density where the Landauparameter F s 0 is still far from its critical value F s 0 = −1. The situation is different in the high-density regime: For sufficiently strong interactions, the system also can develop a high-density condensed phase. The noteworthy feature is that one can get much closer to the spinodal instability limit F s 0 = −1 from above than from below, where F s 0 is always significantly higher than -1. To examine the nature of the instability, we show in Fig. 7 the in-medium scattering length a, see Eq. LJ model, because the system is stable for all densities as discussed above. Similar to the vacuum scattering length, the in-medium scattering length exhibits a singularity, however, the location of the singularity depends obviously on both the density and the interaction model. Evidently, medium corrections to the effective interac-tionsW(k) have the effect that the in-medium scattering length exhibits a singularity as k F is increasing. The larger |a 0 |, the smaller the critical k F value where the the divergence happens, see also the next figure. a/a 0 is universal only for very low densities, where all curves merge into a single curve converging towards unity in the zero density limit. For finite density, a/a 0 , and thus k F a, depends on a 0 , k F , and the interaction model, and not just on the dimensionless parameter k F a 0 .
The appearance of such an effect is not surprising because the leading correction to the bare interaction is phonon-exchange which is attractive, leading to an earlier divergence of the in-medium scattering length a than the bare vacuum scattering length a 0 as the density is increased. Thus, we are led to the conclusion that the instability found in our calculations is an indication of phonon-induced dimerization. Further evidence for the appearance of a phonon-exchange induced dimerized phase will be provided in the next section where we discuss distribution functions. Since a dimerized phase is not within the space of wave functions (2.2), the FHNC-EL equations simply cease to have a solution. A more appropriate variational wave function is perhaps the "Pfaffian" form suggested, among others, in Refs. [32,42,69,70]. Similar dimerization effects have been observed in two-dimensional 3 He-4 He mixtures, in double-layers of bosonic dipoles [71], and also in the pro- cess of α-clustering of nuclear matter [72,73]. Figs. 7 show essentially the same scenario for the SW and the LJ model. At a given density, the critical scattering length a 0 where dimerization occurs is modeldependent. But the difference becomes smaller as the density is reduced. In order to determine the critical scattering length, we have extrapolated the in-medium scattering length to the point where it diverges. The result for both interactions is shown in Fig. 8 where we show, as function of k F σ, the inverse of the critical scattering length where the system becomes unstable for both the SW and the LJ interaction. The two curves approach each other as k F σ is well below 0.1, which shows that the critical value of σ/a 0 becomes universal in the low density limit realized in ultracold quantum gases, although, as a consequence of many-body effects, this critical value of σ/a 0 is not zero.

B. Distribution Functions
The pair distribution functions for parallel and antiparallel spins, g ↑↑ (r) (2.20) and g ↑↓ (r) (2.21), contain information about correlations due to Fermi-statistics and, more interestingly, due to interactions. The latter are captured by the direct-direct correlation function Γ dd (r), which we show in Figs. 9 for the LJ and SW model, respectively. The respective potential strengths were chosen such that the scattering length was a 0 = −2.5 σ. Γ dd (r) is shown for three Fermi wave numbers k F σ = 0.001, 0.01, and 0.04. One can discern two regimes: the asymptotic regime where Γ dd (r) falls off as 1/r 2 for rk F 1 due to many-body effects, see Eq. (3.21); and an intermediate regime k F σ rk F 1 where r is smaller than the average particle distance and Γ dd (r) falls off like 1/r as expected from two-body scattering in vacuum. The behavior in the asymptotic regime can be obtained from the k → 0 limit of Eqs. (2.7) and (2.11) which leads to Eq. (3.21). In the low-density limit, the speed of sound is obtained from the equation of state (2.23). Since only the speed of sound enters, the asymptotic form of Γ dd (r) is independent of the interaction model. In Figs. 9 this asymptotic behavior is illustrated by straight lines.
For low densities, the spin-parallel pair distribution function g ↑↑ (r) is dominated by Fermi statistics. The Pauli exclusion principle ensures that g ↑↑ (0) = 0, regardless of interactions. This is guaranteed by the statistical factor 1 − ℓ 2 (rk F ) in Eq. (2.20) which suppresses g ↑↑ (r) for rk F 1, thus screening the interaction. Note also that lim r→0 C(r) + (∆X ee ) 1 (r) = 0. Therefore, the interaction effectively plays no role in a dilute gas of spin-polarized fermions, which has also been established in many experiments. For example for k F σ = 0.01, the curves for g ↑↑ (r) are almost indistinguishable from each other and from the distribution function of noninteracting fermions for a wide sequence of values of the scattering length between a 0 = −0.5 and −4.0. g ↑↑ (r) is essentially identical to the spin-parallel pair distribution of free fermions g F ↑↑ (r), this has also been verified by QMC results [42,50]. Thus, little or no non-trivial information can be obtained from g ↑↑ (r), and we therefore refrain from further discussing or showing this quantity.
The anti-parallel pair distribution function g ↑↓ (r) is not suppressed by statistical correlations.
Instead, g ↑↓ (r) is dominated by correlation effects, which are described by the direct-direct correlation function Γ dd (r), Eq. (2.11). The effectively attractive interaction, a 0 < 0,  leads to an enhancement of g ↑↓ (r) as r is reduced. We show in Fig. 10 g ↑↓ (r) for k F σ = 0.01 and a 0 /σ = −0.5, −1.0, . . . − 4.0, for both the SW and the LJ potential. For small distances r, g ↑↓ (r) is dominated by the interaction potential, leading to a large value at r = 0 for the SW potential: The maximum value of g ↑↓ is attained at r = 0. For the largest |a 0 | considered here, we have obtained solutions of the FHNC-EL equations where g ↑↓ (0) is more than 150 times larger than g ↑↓ (∞). The short range repulsion of the LJ interaction suppresses the pair distribution for r < σ, and therefore g ↑↓ (r) → 0 as r → 0. The maximum value of g ↑↓ (r), which is attained at finite r, is therefore much lower than in the SW case. One might argue that both models are unrealistic at short r, the shape of g ↑↓ (r) for r σ is of only theoretical interest. However, we will show in Section IV C that finite-range effects have a quite visible effect on the superfluid energy gap and, hence, on the condensation energy. We also observe that for r 3σ, g ↑↓ (r) becomes universal in the sense that it is independent of the choice of interaction. This is not a feature special to cold gases: The asymptotic behavior of g ↑↓ (r) is determined by the speed of sound [44], see also Eq. (3.21). g ↑↓ (r) is still strongly enhanced compared to the asymptotic limit g ↑↓ (r → ∞) = 1. Furthermore we find that g ↑↓ (r) is very similar to the zero-energy scattering solution of the two-body Schrödinger equation, this is simply due to the fact that g ↑↓ (r) is dominated by the dynamic correlation function 1 + Γ dd (r), cf. Eq. (2.21). The situation changes again as we approach the dimerization instability: The maximum value of 1 + Γ dd (r) starts to rise rapidly indicating an approaching singularity. The effect becomes stronger and appears at lower density as the coupling strength of the underlying bare interaction is increased.

C. BCS pairing
We now return to the question of the density-and momentum dependence of the pairing interactionW(k). We have discussed already above that the in-medium scattering length can become singular due to phonon-exchange correction.
Before we describe our calculations we go back to the momentum dependence of the pairing interaction. We show in Figs. 12W(k)/W(0+) for two values of k F , k F σ = 0.01 (top panel) and k F σ = 0.04 (bottom panel). For small k F such as k F σ = 0.01 and less, the regime where Γ dd (r) behaves as 1 + Γ dd (r) = 1 − a0 r is rather large, see Figs. 12 and therefore the linear regime inW(k) is well defined and agrees well with the prediction (3.22). For k F σ = 0.04, the regime where the a 0 /r behavior of the correlations is visible is much smaller, hence the linear regime is less well defined andW(k) also becomes model dependent. The important point to be made here is, however, that in neither case, and especially as |a 0 | gets larger, the simple estimateW(k) ≈W(0+) from lowdensity expansions is valid for a significant portion of the integration regime 0 ≤ k ≤ 2k F which is needed for the calculation of the s-wave pairing matrix elements, see Eq. (3.13).
We have calculated the gap in the excitation spectrum at the Fermi surface, ∆ kF (Eq. (3.11)), for the LJ and the SW interaction model, for a wide range of densities, characterized by the Fermi wave number k F , and s-wave scattering lengths a 0 . A first overview is shown in Fig. 13, where the higher values of ∆ F correspond to larger values of |a 0 |. Figs. 13 provide two pieces of information: One is an assessenent of the accuracy of the approximate solution (3.15) which is obviously qute good throughout the whole regime of interaction strengths and densities. The major difference comes, at hight densities, from the deviation of the vaccum scattering length from a F . The general dependence of the gap on both the interaction strength and the density is quite similar, in particular the exponential dependence on the scattering length holds over many orders of magnitude. This is consistent with the general feature spelled out in Sec. III B that mediumand finite-range corrections are manifested in the prefactor in Eq. (3.24). This pre-factor is universal as long as the in-medium scattering length is of the form (3.23) which is the case, amomg others, if the momentum dependence ofW(k) is linear. Deviations from this universal behavior can only be expected from the quadratic behavior ofW(k) for 0 ≤ k k F which is a direct manification of many-body effects. Indeed, this correlation effect is significant: We show in Fig. 14 the ratio a F /a 0 , see Eq. (3.16), in the density/interaction strengt regimes where the gasp is larger than 10 −10 e F . Evidently, the behavior is, in that regime, neither linear, nor a univer-    sal function of a 0 k F . Only at very small values of a 0 k F where the gap is of the order of 10 −8 e F , a linear behavior might be interpreted, but the slope of such a linear behavior depends on the interaction model. Calculations at even lower density where a linear regime might be found require a much bigger mesh than used here to reliably resolve the region 0 ≤ k k F .
In order to disentangle medium and finite-range corrections, we show in the right and left panels of Figs. 15 (LJ model) and 16 (SW model) the ratios ∆ kF /∆ F and ∆ kF /∆ (0) F as defined in Eqs. (3.15) and (3.17), respectively, plotted as function of −1/k F a 0 . The color of the symbols encodes the k F value of the corresponding data point, between red for k F σ = 0.1 and black for k F → 0 -the low density regime which usually is assumed to be universal in the sense that all quantities depend on −1/k F a 0 only. Data for k F σ > 0.1 is not shown in these figures. The two ratios ∆ kF /∆ F and ∆ kF /∆ (0) F give information on two different effects: The ratio ∆ kF /∆ F is an assessment of the accuracy of the approximations leading to Eqs. (3.15) and reflects the importance of the momentum dependence of the pairing interaction. Evidently, this effect is visible and depends little on the density, the ratio does not seem to go to unity with decreasing density. This is simply a consequence of our analysis of Figs. 12: No matter how small the Fermi momentum is, the momentum dependence of the pairing interaction is never given by the naive application (3.22) of low-density arguments but rather reflects genuine manybody physics. For the range of densities and coupling strengths, ∆ kF /∆ (2) F is between 60% and 75%, but there is a rapid drop of ∆ kF /∆ (2) F for −1/k F a 0 < 5, as the system becomes unstable against dimerization. In other words, neglecting the momentum-dependence of the pairing interaction, and thus the finite range of the interaction (including the induced interaction), leads to an overestimation of ∆ kF that becomes more severe as one approaches the dimerization instability, driven by the divergence of the in-medium scattering length a. We also note that, apart for some values at larger k F , ∆ kF /∆ (2) F collapses on a single curve, i.e. it is characterized by a universal dependence on k F a 0 . The curves are identical for the LJ and the SW model.
The difference between our microscopic calculation and the low-density expression (3.17) is much more visible, see left panels of Figs. 15 and 16. This effect has, of course, to do with the divergence of the in-medium scattering length as discussed above. However, even far from the singularity, we see very little evidence of an approach to a "universal" behavior. Evidently one must go to much lower densities. Although we could run our ground state calculations down to densities of 10 −11 σ −3 , the gap itself becomes many orders of magnitude smaller than the Fermi energy, cf. Figs. 13 and are therefore of only methodological interest.

V. SUMMARY
We have in this paper reported ground state calculations for low-density Fermi gases described by two model interactions, an attractive square-well potential and a Lennard-Jones potential. We have used the optimized Fermi-Hypernetted Chain (FHNC-EL) integral equation method which has been proved to provide, in the density regimes of interest here, an accuracy better than one percent [54]. We have examined the low-density expansion of the energy for a local correlation operator, written in the conventional Jastrow-Feenberg form (2.2), our results also apply for fixed-node Monte Carlo calculations for wave functions written in that form. Of course, such a result can only be obtained with optimized correlation functions, using a parametrized Jastrow function instead leads to an unpredictable answer of which one can only say that it should lie above the expansion (2.31). We have demonstrated that a locally correlated wave function does not reproduce the exact low-density limit and have cured the problem by adding the second-order CBF corrections. We have, however, also demonstrated that already at values of a 0 k F ≈ 10 −2 the third term in the expansion (2.23) is overshadowed by model-dependent corrections.
The most interesting result of our work is the appearance of a divergence of solutions of the FHNC-EL equations well before the divergence of the vacuum scattering length a 0 of the interaction potential. We have interpreted this divergence as a phonon-exchange driven dimerization of the system, similar to what one has at zero density when the vacuum scattering length diverges. The appearance of such a divergence is not a surprise: As a function of coupling strength, the vacuum scattering length has a singularity somewhat above the strength where we find the singularity. However, within the medium, the bare interaction is supplemented by the induced interaction w I (r), cf. Eq. (2.10) or (3.19), which describes phonon exchange. This density dependent correction to the bare interaction shifts the interaction strength at which a bound state appears. The closer the bare interaction is to the appearance of the bound state, the smaller the correction needs to be. Since the induced interaction depends on the density, the singularity appears, for stronger coupling of the bare interaction, at lower density.
We have also studied, in the stable regime, the superfluid gap and its dependence on the density and the interaction strength. Two corrections apply to low density expansions: medium corrections and finite-range corrections. We have shown that the most important finiterange corrections are a direct manifestation of the manybody nature of the system. For low densities, the gap can be reasonably well approximated by neglecting finiterange corrections but accounting for medium corrections, but the deviations from the full numerical solution increase on approaching the dimerization instability.
The natural question arises what is the phase "on the other side" of the instability. It is plausible that this is a phase of dimers which could be described by a product of pair wave functions [32,42,70]. Full optimzation of such wave function containing both Jastrow correlations and an antisymmetrized product of pair wave functions requires first to derive the diagrammtic expansion, which goes beyond the scope of this paper.
Appendix A: Low-density expansions for fermions 1. Low density limit for local correlations In the limit of low densities, the exact energy per particle is given by the expansion (2.23) [28,29]. Since the wave function containing a local correlation operator (2.2) is not exact, it is of interest to determine the consequences of this approximation for the low-density expansion.
There are two ways to derive the low-density expansion for the wave function (2.2): One is to assume that the interaction is weak and has a Fourier transform. Then one can proceed in the same way as in Ref. 37, i.e. • Derive the cluster-expansion of the energy keeping all terms that are (a) of second order in the correlation function f 2 (r) − 1 ≡ e u2(r) , or (b) of first order in the potential and first order in the correlation function • Derive the Euler equation, express the correlation function to first order in the potential • Re-insert the optimal correlation function in the energy • Collect all terms that can be written as powers k F a. If we also permit hard-core interactions, we must formulate the expansion in terms of the dynamical correlation function Γ dd (r) but take into account that the solution of the Euler equation is density dependent. It is sufficient to keep pair-correlations only, three-body correlations contribute only to fourth order in the interaction. The analytic expression for the relevant energy contributions is given in Eqs. (2.14), (2.15) and (2.18). The only additional simplification is that we can set, to leading order in Γ dd (r), the factor S 2 Q . The resulting diagrammatic representation of the energy expansion is shown in Fig. 17.
In the low-density limit, v CW (r) is much shorter ranged than both the exchange function ℓ(rk F ) and the dynamic correlation function Γ dd (r),ṽ CW (k) can therefore be replaced by its long-wavelength limitṽ CW (0+) which is, in turn related to the scattering length, see Eq. (2.30). Then, diagrams (b), (d), and (f) are 1/ν times diagrams (a), (c), and (e), respectively. For the correlation function (2.29), all diagrams shown in Fig. 17 can be calculated exactly in terms of the long-wavelength limitṽ CW (0+). The calculation is somewhat tedious but straightforward and leads to the aforementioned result (2.31). However, this is not the exact solution because the exchange diagram (f) has not been treated exactly.
To derive the exact low-density limit, we start from the diagrammatic expansion up to second order is shown in Fig. 17. Carrying out the variation with respect to 1 + Γ dd (r) where (∆X ′ ee ) 1 (r) is represented by the two diagrams shown in Fig. 1, with the dashed line replaced by v CW (r). In leading order in k F , we can replace g F (r) = g F (0) = 1/2 and expand and similarly which lets us combine the last two terms in Eq. (A1) to Multiplying the equation with 1 + Γ dd (r), using and keeping only first order terms in Γ dd (r) reduces Eq.
The approximation leading to Eq. (2.29) is to replace the third term by This is legitimate if the range of Γ dd (r) is small compared to that of ℓ(rk F ) which is not necessarily true. To assess the importance of this term, we writẽ with σ(k) ∝ k as k → 0 and σ(k) → 1 as k → ∞.
In the long wavelength limit v CW (k) ≈ṽ CW (0+), the function σ(k/k F ) becomes universal. We have solved the Euler equation (A4) in that limit, the resulting function σ(k/k F ) is shown in Fig. 18. The deviation from S F (k) is small but visible. Calculating the energy correction to second order in (ak F ) leads to a coefficient 1.5519 instead of 1.5415. The smallness of the effect is plausible because diagram (g) contributes only a few percent to the total energy correction.

Density Dependence of the Correlation Functions
We have above determined the connection between the dynamical correlation function Γ dd (r) and the effective interaction v CW (r). To complete the calculation of the low-density limit of the energy for local correlation operators, we must establish the connection betweenṽ CW (0+) and the scattering length a 0 becauseṽ CW (0+) is calculated with the optimal correlation function (A5) of the many-body problem at finite density, and not with the solution of the vacuum scattering equation (2.24). To calculate the correction, let δψ(r) = ψ(r) − 1 + Γ dd (r) deltaψ(k) =ṽ CW (k)(S F (k) − 1) 2t(k) .
Theñ v CW (0) = ρ d 3 r The last line follows because v(r)|δψ(r)| 2 is of order O(v 3 ), and the mixed term is zero due to the scattering equation (2.24). Since we have eliminated the potentially singular bare potential, we can expand the kinetic energy term to get v CW (0+) = 4πρ 2 m a 0 + 1 2 Collecting all results, one finds ak F π + 1.53517 ak F π 2 + . . . , (A8) see Appendix A 1 for details of the calculation. The numerical factor 1.53517 is to be compared with the factor 4(11 − 2 ln 2)/35 = 1.098 of Eq. (2.23). We emphasize again that the result (2.31) also applies for the "fixednode" approximation because the terms where that approximation deviates from our expansion are of at least fourth order in the potential strength.

Correlated Basis Functions Corrections
Once we have determined that a local correlation operator of the form (2.2) does not lead to the correct lowdensity limit (2.23) of the correlation energy, the question arises what corrections to that wave function are necessary. We show here that second-order CBF perturbation theory (2.44) is sufficient to retrieve the exact low-density expansion. In the low-density (local) approx-imationW(q) andÑ (q) are given by Eqs. (2.43).
There are normally significant numerical cancellations between the two terms in the numerator of Eq. (2.44). However, for the purpose of discussing formal properties, we may expand the numerator. The term containing two matrix-elements pp ′ | W |hh ′ a is just the ordinary expression of second-order perturbation theory [28]. None of the other terms contains an energy denominator, they can therefore be written in terms of the diagrams of the JF method.
We carry out the calculation for just the direct term, notice that the double-exchange is equal to the direct term.
where V F = 4π 3 k 3 F ,n(k) = 1 − n(k), and (A10) Expanding the second term as we identify in the first term the negative diagram (g) of Fig. 17. Thus, these terms cancel exactly all direct, two-line diagrams shone in Fig. 17. The last terms in Eqs. (A11) and (A12) can not be represented by legitimate diagrams of JF-theory, they can be are combined to − 1 4 d 3 q (2π) 3 ρ t(q)Γ 2 dd (q) + 2ṽ CW (q)Γ dd (q) = − 1 4 The first term is seen to cancel the correction (A7), whereas the second term is combined with the first term in Eq. (A9) to give the textbook result [28] − 1 4 The treatment of the exchange terms is done along the same lines and will not be repeated here.
Appendix B: Numerical solution of the gap equation

General strategy
We present in this appendix details of our solution method for the gap equation. Besides describing the sometimes delicate numerics, we obtain information on the accuracy of approximate analytic solutions [65,66] which are often used to discuss superfluidity or superconductivity.
The gap equation is highly non-linear and a simple iteration procedure does not converge. Moreover, in particular for values ∆(k) ≪ e F , the integrand is narrowly peaked at the Fermi momentum and an adaptive mesh is necessary to reach a reliable numerical accuracy. Above, A very rapidly converging algorithm is as follows: Consider the generalized eigenvalue problem λ(ξ)∆ n+1 (k, ξ) = (B3) To shorten the equations, let E n (k, ξ) ≡ (ǫ k − µ) 2 + ξ 2 ∆ 2 n (k) .
The problem can be mapped to a regular symmetric eigenvalue problem by defining Then, the eigenvalue problem is The algorithm is (i) Start with a reasonable estimate ∆ 0 (k ′ ), e.g. a constant (ii) Solve the above eigenvalue problem as a function of the scaling parameter ξ, find the value ξ 1 for which the lowest eigenvalue of the equation is 1. Along with a trial solution, we can calculate the derivative with respect to ξ, see below.
(iv) Go to step (ii) and repeat until convergence which is obtained for ξ = 1.

Adaptive mesh
To deal with the strongly peaked denominator E(k, ξ), we introduce a transformation k = k(x) that generates a dense set of points around the Fermi momentum. Then (B6) The idea is now to find a function k(x) such that dk dx 1 (ǫ k − µ) 2 + ∆ 2 (k) (B7) is almost constant. The function of choice is where δ = ∆(k F )/e F and κ = k/k F . This choice has the property dx dk = k e F (k 2 /k 2 F − 1) 2 + δ 2 (B10) i.e.