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

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, Phys. Rev. 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 , where is the vacuum -wave scattering length and 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 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 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.


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]. 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]. 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 [12], this so-called BCS-BEC cross-over has been widely studied experimentally [8] and theoretically [20], see also reviews Refs. . The observation of quantized vortices on both sides of the BCS-BEC cross-over provided an unambiguous proof of superfluidity by fermionic pairing [25]. Recent work investigated the effect of partial polarization of a two-component Fermi gas on the Fermi liquid parameters [26], the nature of the transition from a BCS state to a state of a molecular BEC [27], and quantifying the superfluid fraction in a Fermi gas by means of second sound [28].

We study in this work the low-density properties of homogeneous Fermi gases at zero temperature. We use for our study a square-well and a Lennard-Jones interaction potential. Changing the interaction strength (coupling constant) of the respective potential changes the scattering length for two-body scattering , which we refer to as vacuum scattering length. The in-medium scattering length will be introduced later. When , i.e. the interaction is effectively attractive, one expects BCS type pairing of particles with opposite spin. As , a low energy resonance of the two-body problem generates bound dimers. is called the unitary limit, since the only relevant length scale is the inverse Fermi wave number . Increasing the attraction further, the Cooper pairs become bound molecules and the fermionic nature of their constituents becomes less visible.

We are in particular interested in structural quantities such as the energetics, distribution functions, the stability of the system against spinodal decomposition, dimerization, and 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 [29] that are routinely applied at low densities.

In the low density limit, many quantities like, for example, the ground state energy, depend only on the dimensionless parameter [30]. We are interested here in the parameter range where this “universal” behavior ceases to persist due to correlation effects. An example for correlation effects is the pair distribution function, . In mean–field approaches, is equal to the distribution function of the non-interacting Fermi gas, . As we will see below, deviates, in particular for spin-antiparallel particles, substantially from when the absolute value of the scattering length becomes large compared to the characteristic length of the interaction potential.

In the limit of weak attractive interactions, the system can be described by a BCS type wave function. When the weak-coupling approximation does not apply (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 [31] is reviewed in Section 3.1, see also Refs. 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 is a full Fermi Hypernetted Chain (FHNC) summation for large gap parameters has been suggested by Fantoni [36], we will comment on this approach further below.

Our paper is organized as follows: In Section 2 we will review briefly 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 [38], from Coupled Cluster theory [39] and from an extension of density functional theory which includes pair correlations. We evaluate in Section 2.2 the low-density limit and show that the exact formula [30] is not reproduced by the Jastrow-Feenberg and/or the “fixed node” approximation. To correct this problem, we apply in section Section 2.3 perturbation theory in a correlated basis. We show that second-order CBF corrections must be added to obtain the correct low–density expansion.

In Section 3.1 we will review the CBCS theory. We will show that the theory can be formulated in exactly the same way as ordinary BCS theory. CBCS simply provides a prescription 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 -matrix used in the low-density expansion of BCS theory [40].

In Section 4, we present our results for the energy, the pair distribution function, the in-medium scattering length, and the gap energy as a function of Fermi wave number and the vacuum scattering length . The dynamical correlations can be characterized by three regimes: For short distances, , these correlations are, of course, determined by the interaction. The intermediate regime is dominated by two-body scattering where correlations decay as . A third, asymptotic regime for larger than the average interparticle distance is dominated by many-body effects where the correlations decay as .

We find an instability of the FHNC-Euler-Lagrange (FHNC-EL) solutions accompanied by a divergence of the in-medium scattering length . This instability occurs well before the divergence of the vacuum scattering length ; it is caused by the induced interaction mediated by phonon exchange. Thus, strongly bound dimers can be formed at finite density even if the bare potential does not have a bound state. Note that this is a property of the normal system, wave functions of the type normally used only in the “BEC” regime [43] would here be more appropriate. We will return to this point further below.

In the regime we solve the CBCS gap equation and show that deviations from the simple BCS approximation can be separated into two contributions. One stems from the density dependence of the in-medium scattering length, the other one from the non-negligible momentum dependence of the pairing interaction in the CBCS gap equation.

2Generic Many-Body Theory

2.1Variational 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

where is a phenomenological, local interaction. The method starts with a variational ansatz for the wave function [46]

is a model state, normally a Slater-determinant for fermions and for bosons, and is the correlation operator written in the form ( ?). There are basically two ways to deal with this type of wave function. In quantum Monte Carlo studies, the wave function (Equation 2) is referred to as “fixed node approximation”. An optimal correlation function is obtained by stochastic means [43]. A decomposition into -body correlations is then, of course, not necessary. Alternatively, one can use diagrammatic methods, specifically the optimized Fermi-hypernetted chain method for the calculation of physically interesting quantities. These diagrammatic methods have been successfully applied to such highly correlated Fermi systems as He and He at [55]. They are naturally expected to work much better in the low density systems of interest here. In fact, we have shown in recent work [56] that even the simplest version of the FHNC-EL theory is accurate within better than one percent at densities less than 25 percent of the ground state density of liquid He Diffusion Monte Carlo calculations typically use a parametrized Jastrow-Feenberg (JF) ansatz for importance sampling, where the parameters are optimized by variational Monte Carlo calculations. JF theory makes explicit use of the form ( ?). It has been shown [38] 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 [57] in both isotopes. They are completely negligible below approximately 25 percent of the respective equilibrium densities. We can identify, at the low densities we are concerned with here, the Jastrow-Feenberg approximation with the fixed-node approximation in quantum Monte Carlo calculations.

The correlations are obtained by minimizing the energy, i.e. by solving the Euler-Lagrange (EL) equations

The evaluation of the energy (Equation 3) for the variational wave function ( ?) and the analysis of the variational problem are carried out by cluster expansion and resummation methods. The procedure has been described at length in review articles [58] and pedagogical material [32]. 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 of approximations is the only systematic approximation scheme that preserves the properties of the exact variational problem [46].

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 [55]. 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 for the wave function Eq. (Equation 2).

In the FHNC-EL//0 approximation [59], which contains both the random phase approximation (RPA) and the Bethe-Goldstone equation in a “collective” approximation [32], the Euler equation ( ?) can be written in the form

where is the static structure factor of the interacting system, is the kinetic energy of a free particle, is the static structure of the non-interacting Fermi system, and

is the so-called “particle-hole interaction”. We have also above defined the “Clark-Westhaus effective interaction” [58]. As usual, we define the Fourier transform with a density factor,

Auxiliary quantities are the “induced interaction”

and the “direct-direct correlation function”

Eqs. (Equation 4)–(Equation 8) form a closed set which can be solved by iteration. Note that the Jastrow correlation function ( ?) has been eliminated entirely.

Figure 1: The two diagrams contributing to (\Delta X_{\rm ee})_1(r). The usual diagrammatic notations of FHNC-EL theory  apply.
Figure 1: The two diagrams contributing to . The usual diagrammatic notations of FHNC-EL theory apply.

The pair distribution function can generally be written as

Roughly speaking, describes dynamic, short-range correlations, is the pair distribution function of the non-interacting Fermi gas, and is the Slater exchange function. describes the combination of statistical and dynamic correlations. In leading order in the dynamic correlations we have

where is represented by the two exchange diagrams shown in Figure 1.

In this approximation, the energy per particle is [55],

The term 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 spin-antiparallel channel:

where indicates the Fourier-transform (Equation 6). To leading order in the density, the term is the only term that reflects Fermi statistics, whereas the factor describes dynamical correlations. , , and are normalized such that they go to unity for large .

2.2Low-density limit

In the limit of low densities, the equation of state and related quantities depend only on the vacuum -wave scattering length and the Fermi wave number . For example, the energy per particle has the expansion [30]

Note that the expansion (Equation 12) is strictly valid only for . For attractive potentials the superfluid condensation energy must be added.

The locally correlated wave function (Equation 2) is not exact, and the question arises whether it recovers the expansion (Equation 12). It is plausible that this is not the case: The calculation of the third term in Eq. (Equation 12) makes explicit use of the form of the energy denominator in second order perturbation theory [29]. 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 (Equation 10) to second order in the vacuum scattering length . 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. for fermions.

We will show the details of the calculation in Appendix 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. (Equation 13) with and using the identity

gives a relationship that will be useful later:

The quantity is structurally identical to as introduced in Eq. (Equation 5), except that is calculated for the zero energy vacuum scattering solution . Integrating Eq. (Equation 16) leads to the relationship

We notice that the induced interaction as defined in Eq. (Equation 7) is of second order in the interaction. To leading order in the density we can also expand Eq. (Equation 4)

and obtain from Eq. (Equation 8) the solution

In addition to calculating the energy contributions (Equation 10) for the correlation function (Equation 18), we must express in terms of the scattering length because is calculated with the optimal correlation function (Equation 8) of the many-body problem at finite density, and not with the solution of the zero-density scattering equation (Equation 13). In Appendix A.1 we will prove the relationship

Collecting all results, one finds

see Appendix A.1 for details of the calculation. The result (Equation 20) is to be compared with the factor of Eq. (Equation 12). 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 described in the next Section 2.3. 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.05690Ry [61] instead of the exact value 0.06218Ry [62]. This deficiency is, for the electron gas, removed by second-order CBF theory [64]. One conclusion of our analysis is that the fixed node approximation in conjunction with Green’s functions [47] or Diffusion [44] should reproduce the expansion (Equation 20) and not (Equation 12).

2.3Elements of Correlated Basis Functions

We have seen above that a locally correlated wave function (Equation 2) does not produce the exact low-density limit (Equation 12) 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 the method only very briefly, details may be found in pedagogical material [32] and review articles [58]. The diagrammatic construction of the relevant ingredients has been derived in Ref. .

CBF theory uses the correlation operator to generate a complete set of correlated and normalized -particle basis states through

where the form a complete basis of model states, normally consisting of Slater determinants of single particle orbitals. Although the are not orthogonal, perturbation theory can be formulated in terms of these states [66].

For economy of notation, we introduce a “second–quantized” formulation of the correlated states. The Jastrow–Feenberg correlation operator in ( ?) depends on the particle number, i.e. (whenever unambiguous, we omit the corresponding subscript). Starting from the conventional that create and annihilate single particle states, new creation and annihilation operators of correlated states are defined by their action on the correlated basis states:

According to these definitions, and obey the same commutation rules as the creation and annihilation operators and of uncorrelated states, but they are not Hermitian conjugates. If is an –particle state, then the state in Eq. (Equation 22) must carry an -particle correlation operator , while that in Eq. ( ?) must be formed with an –particle correlation operator .

In general, we label “hole” states which are occupied in by , , , and unoccupied “particle” states by , , etc.. To display the particle-hole pairs explicitly, we will alternatively to the notation use . A basis state with particle-hole pairs is then

For the off–diagonal elements of an operator , we sort the quantum numbers and such that is mapped onto by

From this we recognize that, to leading order in , any depends only on the difference between the states and , and not on the states as a whole. Consequently, can be written as matrix element of a -body operator

(The index indicates antisymmetrization.)

The key quantities for the execution of the theory are diagonal and off-diagonal matrix elements of unity and ,

Eq. ( ?) defines a natural decomposition [65] of the matrix elements of into the off-diagonal quantities and and diagonal quantities .

To leading order in the particle number, the diagonal matrix elements of become additive, so that for the above -pair state we can define the CBF single particle energies

with where

and is an average field that can be expressed in terms of the compound diagrammatic quantities of FHNC theory. According to (Equation 25), and define particle operators and , e.g.

Diagrammatic representations of and have the same topology [65]. In homogeneous systems, the continuous parts of the are wave numbers ; we abbreviate their difference as .

In principle, the and are non-local -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 Section 2.1. In the same approximation, the operators and are local, and we have [55]

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 [55]. 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 ( ?) corrects the coefficient of the third term in the expansion (Equation 20) and leads to the exact low-density limit (Equation 12).

3BCS Theory with correlated wave functions

3.1General derivation

We show in this section how the variational theory is generalized to a superfluid or superconducting state. We restrict ourselves here to the simplest case of –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. .

The BCS theory of fermion superfluidity generalizes the Hartree–Fock model by introducing a superposition of independent particle wave functions corresponding to different particle numbers [67]

The coefficient functions and are known as Bogoliubov amplitudes. They 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 for small interparticle distances . 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 to . The former is defined in the –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, applying the correlation operator to that state, normalizing the result, and finally summing over all particle numbers. We must therefore distinguish between correlation operators and normalization integrals corresponding to different particle numbers . Thus, the correlated BCS (CBCS) state is

The trial state (Equation 32) superposes the correlated basis states with the same amplitudes with which the model states enter the corresponding expansion of the original BCS vector.

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 (Equation 33) for the superfluid trial state (Equation 32). This has been done successfully for the one– and two–body density matrices corresponding to a slightly different choice of the correlated BCS state [36] 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 , where is the superfluid gap energy. We will demonstrate below that this quantity is indeed small in the regime where the wave function (Equation 32) is appropriate.

In leading order, it is sufficient to retain the terms of first order in the deviation and those of second order in . We refer to this as the “decoupling approximation”. The calculation of for correlated states [31] 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,

With the results ( ?) and ( ?), 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 [41]. The correlation operator serves here to tame the short–range dynamical correlations. The effective interaction is just an energy independent approximation of the -matrix.

We may proceed now in the conventional way to determine the Bogoliubov–amplitudes , , by variation of the condensation energy ( ?) 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 , . The CBCS gap equation becomes

The conventional, i.e. uncorrelated, BCS gap equation [68] is retrieved by replacing the effective interaction by the matrix elements of the bare interaction. Note that our “decoupling approximation” simply means that we assume that the pairing interaction does not depend on the Bogoliubov amplitudes.

To conclude this section, we point to a subtle issue concerning the normalization of the correlated BCS state (Equation 32). Above, we have written the wave function as a specific linear combination of normalized states . In our formulation of a correlated pairing theory [34] as well as in related work [36] the correlated BCS state


is the BCS state (Equation 31) projected into the -particle Hilbert space. This approach leads to a slightly different gap-equation where the gap function under the square root in Eq. (Equation 35) is scaled by a factor , see Eq. (3.16) od Ref. . This factor diverges for optimized or otherwise long-ranged correlations. To avoid this divergence we have used here the method developed in Ref. .

3.2Analysis of the gap equation

In the local approximations appropriate for low densities, the effective interaction is given by Eqs. (Equation 30). The pairing matrix element is expressed in terms of the Fourier-transforms and :

The remaining arguments are standard, cf. Ref. : If the gap at the Fermi surface is small, we can replace the pairing interaction by its -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. . In particular, the second term has the only function to regularize the integral for large . 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 with the vacuum scattering length :

Of course, our equations ( ?), (Equation 35) are much more general: At low densities, the subtraction term – i.e. the second term in the square bracket of Eqs. (Equation 36) and (Equation 38) is important to regularize the integral for large momentum transfers. At higher densities, the finite range of the interaction provides that momentum cutoff and the subtraction term becomes negligible since the energy numerator term is zero at the Fermi momentum.

By comparison with the low-density limit and Eq. (Equation 17) we interpret the constant

as an “in-medium” scattering length. Hence, at finite densities, one expects two types of corrections:

  1. Medium corrections: The effective pairing interaction is related to through Eqs. (Equation 4)- (Equation 8) which leads to

    Because of Eq. (Equation 19) and the fact that the induced interaction is of second order in the interaction, we conclude that that

    In the same order, non-local contributions to the pairing interaction (Equation 30) [65] 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. [40]. Moreover, similar to our analysis of the low-density limit, local correlation functions will not get the right coefficient of the term proportional to , hence CBF corrections to the pairing interaction [31] will also lead to modifications of order .

  2. The solution of the -wave gap equation is dominated by the matrix element (Equation 37) of at the Fermi surface which leads to the solution (Equation 39). Only if is practically constant for , we can identify with the in-medium scattering length . The dominant finite-range correction to the pairing interaction comes from the kinetic energy term . For interparticle distances much larger than the interaction range , but smaller than the average particle distance, this term is dominated by the vacuum solution of the Euler equation, . For interparticle distances larger than , we obtain from Eqs. (Equation 4) and (Equation 8) that falls off like

    Consequently, the effective interaction is quadratic in for and has a linear dependence

    for . The variation of the pairing interaction between and is interaction-dependent and causes, as we shall see in Section 4.3 a correction that is larger than the ones due to the effects mentioned above.

All of the corrections discussed above except the one due to the quadratic momentum dependence of for are order and higher i.e. thy lead to a density dependence of the in-medium scattering length of the form

where is a numerical constant. Among others, the polarization correction discussed by Gorkov [40] has this structure. Inserting the above expansion in (Equation 39) changes the pre-factor to

In other words, to the extent that an expansion in powers of is legitimate, all of these corrections just lead to a modified, but universal, pre-factor in Eq. (Equation 41). This does not apply to the finite-range correction of the pairing interaction in the relevant regime . One would, of course, expect that this finite-range correction is of the same order of magnitude. We shall see in Section 4.3 that its value depends, on details of the interaction. Hence, we conclude that the exponential behavior of Eq. (Equation 41) is universal whereas the pre-factor is not.



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 characteristic length and the depth of the attractive well . In both cases, we measure energies in units of , 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 -wave scattering length . The procedure is legitimate in the low-density limit, many observable properties of these gases, such as the energy (Equation 12), depend indeed only on the -wave scattering length [30]. 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 -wave scattering length instead of the well-depth to characterize the potential.

(color online) The plot shows the scattering length 00 as a function of the interaction strength for the LJ (red, solid) and the SW (blue, dashed) potential. The vertical lines (at \epsilon =
11.18 for LJ and \epsilon = 4.336 for SW) indicate the interaction strength where a two-body bound state appears. The dots on the lines indicate the interaction strengths corresponding to vacuum scattering lengths 00/\sigma = -0.5, -1.0,\ldots\-4.0 for which we highlight the Landau parameter and the in-medium scattering lengths in Figs.  and .

For the SW potential, is given by where . For the LJ potential, must be obtained numerically. We show in Fig. ? as a function of the potential well depth . The attractive SW potential has negative scattering length for an interaction strength below the first resonance. The LJ potential has a positive below , indicating an effectively repulsive interaction. Thus, when the interaction strength of the LJ potential is raised, starting from , we find three regimes. (i) For we have and there is no bound state; the many-body ground state is a normal Fermi gas. (ii) For we have . There is still no two-body bound state. Due to the effectively attractive potential, the many-body ground state is, at low densities, a BCS state. (iii) At the LJ potential has a resonance at zero scattering energy. For , becomes positive again, and the potential supports at least one two-body bound state. 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, [56] we have already studied both, and for fermions, and 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 Section 2.1 is reliable. In short, the accuracy of our energy calculations is expected to be better than 1 percent below a density of , whereas the error of the simple FHNC-EL version is about 10 percent as the density increased to which is close to the freezing density of He We focus in the present work on the effect of attractive interactions (). We have solved Eqs. (Equation 4)-(Equation 8) on a mesh of points, with a resolution of 30 points between and , 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 corresponds to a wavelength of , 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 appears (indicated by vertical lines in Fig. ?).

Our equation of state for the two potential models is shown in Figs. ? and ?. To recover the exact low-density limit (Equation 12), we have added the second-order CBF correction ( ?). To emphasize the interaction terms, we have subtracted the kinetic energy . We have normalized the equation of state to the expansion (Equation 12). Thus, Figs. ? and ? show only the model-dependent correction to the equation of state. Omitting the CBF correction ( ?) and comparing to the low-density expansion (Equation 20) gives practically the same results, they are therefore not shown.

Figs. ? and ? show already a visible dependence of the equation of state on the potential model at approximately where the third term in the expansion (Equation 12) is of the order of . In other words, for both interactions the model-dependent corrections are of the same order of magnitude as the third term in the expansion (Equation 12). 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 .

(color online) The plot shows the interaction contribution to the equation of state, normalized to the low density expansion, i.e. the second and third term in the expansion () for the square-well potential. We show results for scattering lengths 00/\sigma = -1,-2,-3,-4, the symbols indicate the numerical values and the curves a second-order polynomial fit of the form E/E_0 = 1+\alpha (00k_{\rm F}) + \beta
  (00k_{\rm F})^2. E_{\rm HY} is the expansion ().
(color online) The plot shows the interaction contribution to the equation of state, normalized to the low density expansion, i.e. the second and third term in the expansion () for the square-well potential. We show results for scattering lengths , the symbols indicate the numerical values and the curves a second-order polynomial fit of the form . is the expansion ().
(color online) Same as Fig.  for the Lennard-Jones model of the interaction.
(color online) Same as Fig. for the Lennard-Jones model of the interaction.

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 . 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 is the hydrodynamic speed of sound, goes to zero. Such an effect has already been reported by Owen [69]. 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, . In the following we will argue that dimerization indeed occurs, but well before the limit . 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. (Equation 4) must be positive. In the limit this leads to the condition

The right-most expression is the low density limit, , where is small and the in-medium scattering length is well approximated by the vacuum scattering length , see Eq. (Equation 52). The limit can be regarded as the low density limit of the particle-hole interaction,

We have used above the fact 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. (Equation 51). Note that the above stability limit (Equation 50) is valid for the local correlation operator (Equation 2). In an improved calculation that does not rely on a local correlation operators but rather includes CBF corrections to all orders [55] the stability condition would read

The condition is immediately recognized as the stability condition of Landau’s Fermi Liquid theory, . Note that the ground state theory formulated in Eqs. (Equation 4)-(Equation 8) 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. , this is an acceptable approximation at the low densities under consideration here. Since the interaction contribution to is, at low densities, proportional to , 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 . In fact, in our calculations we were not able to come close to the point of spinodal instability. We show in Figs. ? and ? the value of 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 up to a critical density that depends on the coupling strength. At that density, begins to drop very rapidly but it does not appear to approach the critical value of for a spinodal instability, see Fig. ?. We also note that an instability signified means a transition to a state with non-uniform density which is not what has been observed experimentally.

(color online) The plot shows the dependence of the Landau parameter F_0^s for the attractive square-well potential as a function of the density for a sequence of coupling strengths \epsilon = 0.1, 0.2, \ldots, 4.6 (blue, dashed curves) and vacuum scattering lengths 00/\sigma = -0.5,-1.0,\ldots,-4.0 (red, solid curves) that correspond to the dots in Fig. . The curve that ends at the lowest density corresponds to the strongest interaction \epsilon = 4.6,\ 00/\sigma=-11.2\,, whereas the ones corresponding to the weak interactions 0.1\le\epsilon\le 0.7 are stable beyond a density of \rho=0.4\sigma^{-3}.
(color online) The plot shows the dependence of the Landau parameter for the attractive square-well potential as a function of the density for a sequence of coupling strengths (blue, dashed curves) and vacuum scattering lengths (red, solid curves) that correspond to the dots in Fig. . The curve that ends at the lowest density corresponds to the strongest interaction , whereas the ones corresponding to the weak interactions are stable beyond a density of .

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 He At low densities, we see the same picture as for the SW potential: The liquid is stabilized by the Pauli pressure, but the FHNC-EL equations cease to have solutions above a certain density where the Landau-parameter is still far from its critical value . 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 from higher densities than from lower densities.

(color online) Same as Fig.  for the Lennard-Jones model of the interaction. The blue (dashed) lines correspond to coupling strengths \epsilon = 7.0, 9.0, 9.2,
  9.3,\ldots 9.9. Note that the Lennard-Jones model also supports a high-density condensed phase: The curves for the interaction strengths \epsilon = 7.0\ (00/\sigma=-1.12) (blue, dashed), \epsilon = 7.51\ (00/\sigma=-1.5) and \epsilon =
  8.01\ (00/\sigma=-2) (both red, solid) are discontinuous at high density. The curve that ends at the lowest density corresponds to the strongest interaction.
(color online) Same as Fig. for the Lennard-Jones model of the interaction. The blue (dashed) lines correspond to coupling strengths . Note that the Lennard-Jones model also supports a high-density condensed phase: The curves for the interaction strengths (blue, dashed), and (both red, solid) are discontinuous at high density. The curve that ends at the lowest density corresponds to the strongest interaction.

To examine the nature of the instability, we show in Fig. ? the in-medium scattering length , see Eq. (Equation 42). The ratio is shown as function of for increasing values of : for the SW model and for the LJ model. We do not show and for the 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. The location of the singularity depends obviously on both the density and the interaction model. Evidently, medium corrections to the effective interactions have the effect that the in-medium scattering length exhibits a singularity as is increasing. The larger , the smaller the critical value where the the divergence happens, see also Fig. ?. 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, , and thus , depends on , , and the interaction model, and not just on the dimensionless parameter .

The appearance of such an effect is not surprising because the leading correction to the bare interaction is the attractive phonon-exchange. This leads, as the density is increased, to a divergence of the in-medium scattering length . Thus, we conclude 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.