Variational and Parquet-diagram theory for strongly correlated normal and superfluid systems

TL;DR

This paper develops a variational and parquet-diagram theoretical framework for strongly correlated normal and superfluid systems, emphasizing the importance of including particle-hole propagators for physical accuracy and applying the theory to neutron matter and Fermi liquids.

Contribution

It introduces a unified variational and parquet-diagram approach, extending FHNC-EL equations beyond traditional approximations for superfluid systems.

Findings

Inclusion of particle-hole propagators improves physical consistency.

Quantitative superfluid gap changes are small, but stability ranges vary significantly.

Application to neutron matter and Fermi liquids reveals system-specific differences.

Abstract

We develop the variational and correlated basis functions/parquet-diagram theory of strongly interacting normal and superfluid systems. The first part of this contribution is devoted to highlight the connections between the Euler equations for the Jastrow-Feenberg wave function on the one hand side, and the ring, ladder, and self-energy diagrams of parquet-diagram theory on the other side. We will show that these subsets of Feynman diagrams are contained, in a local approximation, in the variational wave function. In the second part of this work, we derive the fully optimized Fermi-Hypernetted Chain (FHNC-EL) equations for a superfluid system. Close examination of the procedure reveals that the na\"ive application of these equations exhibits spurious unphysical properties for even an infinitesimal superfluid gap. We will conclude that it is essential to go {\em beyond\/} the usual Jastrow-Feenberg approximation and to include the exact particle-hole propagator to guarantee a physically meaningful theory and the correct stability range. We will then implement this method and apply it to neutron matter and low density Fermi liquids interacting via the Lennard-Jones model interaction and the P\"oschl-Teller interaction. While the quantitative changes in the magnitude of the superfluid gap are relatively small, we see a significant difference between applications for neutron matter and the Lennard-Jones and P\"oschl-Teller systems. Despite the fact that the gap in neutron matter can be as large as half the Fermi energy, the corrections to the gap are relatively small. In the Lennard-Jones and P\"oschl-Teller models, the most visible consequence of the self-consistent calculation is the change in stability range of the system.

Equations

$$H=-\sum _i\frac{{\hslash }^2}{2m}{\nabla }_i^2+\sum _{i<j}v(i,j).$$

$$\displaystyle\bigl{|}{\Psi_{\bf o}^{(N)}}\bigr{\rangle}$$

$$F_N({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)$$

$$H_{{\bf o}}=\bigl{\langle}{\Psi_{\bf o}}\bigr{|}H\bigl{|}{\Psi_{\bf o}}\bigr{\rangle}$$

$$\rho_{2}({\bf r},{\bf r}^{\prime})=\bigl{\langle}{\Psi_{0}}\bigr{|}\sum_{i\neq j}\delta({\bf r}_{i}-{\bf r})\delta({\bf r}_{j}-{\bf r}^{\prime})\bigl{|}{\Psi_{0}}\bigr{\rangle}\,,$$

$$g(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{{\rho }_2(\mathbf{r},{\mathbf{r}}^{\prime })}{{\rho }^2},$$

$$S(k)=1+\rho \int d^{3}r{e}^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}}\left[g(r)-1\right],$$

$$\frac{\delta H_{\mathbf{o}}}{\delta u_2}({\mathbf{r}}_1,{\mathbf{r}}_2)=0.$$

$$g(r)$$

$$X(r)$$

$$\tilde{N}(k)$$

$$\tilde{f}(k)=\rho \int d^{3}r{e}^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}}f(r).$$

$$F{\nabla }^{2}F=\frac{1}{2}({\nabla }^2{F}^2+F^2{\nabla }^2)+\frac{1}{2}{F}^2\left[\nabla ,\left[\nabla ,\ln F\right]\right]-\frac{1}{4}\left[\nabla ,\left[\nabla ,F^2\right]\right].$$

$$F^2\left[\nabla ,\left[\nabla ,\ln F\right]\right]=\frac{1}{2}{F}^2\sum _{i<j}{\nabla }^2{u}_2(r_{ij}).$$

$$\frac{H_{\mathbf{o}}}{N}$$

$$\equiv$$

$$\frac{\delta }{\delta \sqrt{g(r)}}\frac{H_{\mathbf{o}}}{N}=0.$$

$$S(q)$$

$$V_{\mathrm{p}-\mathrm{h}}(r)$$

$${\tilde{w}}_{\mathrm{I}}(k)$$

$$\Delta V_e(r)=\frac{{\hslash }^2}{4m}{\nabla }^{2}E(r)+\rho \int d^3{r}^{\prime }\frac{\delta E(r^{\prime })}{\delta g(r)}\frac{{\hslash }^2}{4m}{\nabla }^{2}g(r^{\prime })\equiv \frac{{\hslash }^2}{4m}{\nabla }^{2}E(r)+E^{\prime }(r)$$

$$\frac{{\hslash }^2}{m}{\nabla }^2\sqrt{g(r)}=\left[v(r)+\Delta V_e(r)+w_{\mathrm{I}}(r)\right]\sqrt{g(r)}.$$

$$\frac{E}{N}=\frac{T}{N}+\frac{V}{N},$$

$$\frac{V}{N}=\frac{\rho }{2}\int d^{3}rv(r)g(r)$$

$$\frac{d}{d\lambda }\frac{E}{N}=\frac{\rho }{2}\int d^{3}rv(r)g_{\lambda }(r)+\frac{1}{N}\int d^{3}r\frac{\delta E}{\delta g_{\lambda }(r)}\frac{d{g}_{\lambda }(r)}{d\lambda }.$$

$$\frac{E}{N}=\frac{E_0}{N}+\frac{\rho }{2}\int d^{3}rv(r){\int }_0^{1}d\lambda g_{\lambda }(r),$$

$$S(q)=-{\int }_0^{\infty }\frac{d\omega }{\pi }\Im m\chi (q,\omega )$$

$$\chi (q,\omega )=\frac{{\chi }_0(q,\omega )}{1-{\tilde{V}}_{\mathrm{p}-\mathrm{h}}(q){\chi }_0(q,\omega )}$$

$$V_{\mathrm{p}-\mathrm{p}}(r)=v(r)+w(r).$$

$$\frac{{\hslash }^2}{m}{\nabla }^2\sqrt{g(r)}=V_{\mathrm{p}-\mathrm{p}}(r)\sqrt{g(r)}.$$

$${\chi }_0^{\mathrm{Bose}}(q,\omega )=\frac{2t(q)}{(\hslash \omega +\mathrm{i}\eta {)}^2-t^2(q)}.$$

$${\tilde{V}}_{\mathrm{p}-\mathrm{h}}(0+)=m{c}^2,$$

$$S(k)=\frac{\hslash k}{2mc}\mathrm{as}k\to 0+.$$

Full text

In the second part of this work, we derive the fully optimized Fermi-Hypernetted Chain (FHNC-EL) equations for a superfluid system. Close examination of the procedure reveals that the naïve application of these equations exhibits spurious unphysical properties for even an infinitesimal superfluid gap. We will conclude that it is essential to go beyond the usual Jastrow-Feenberg approximation and to include the exact particle-hole propagator to guarantee a physically meaningful theory and the correct stability range.

We will then implement this method and apply it to neutron matter and low density Fermi liquids interacting via the Lennard-Jones model interaction and the Pöschl-Teller interaction. While the quantitative changes in the magnitude of the superfluid gap are relatively small, we see a significant difference between applications for neutron matter and the Lennard-Jones and Pöschl-Teller systems. Despite the fact that the gap in neutron matter can be as large as half the Fermi energy, the corrections to the gap are relatively small. In the Lennard-Jones and Pöschl-Teller models, the most visible consequence of the self-consistent calculation is the change in stability range of the system.

1. 1 Introduction

2. 2 Motivation: Variational and local parquet theory for bosons

3. 2.1 Rationalization of parquet diagram summations

4. 2.2 Hypernetted Chain and Euler equations

5. 2.3 Pair density functional theory: The view from the top

6. 2.4 Stability and Consistency

7. 3 Variational and local parquet diagram theory for fermions

8. 3.1 Generating functional and the generalized distribution functions

9. 3.2 Cluster expansions

10. 3.3 FHNC and Euler equations

11. 3.4 Energy

12. 3.5 Uniform limit approximation

13. 3.6 The low-density limit

14. 3.7 Exchange corrections

15. 3.8 Limitations of local correlation functions

16. 3.9 Elements of correlated basis functions

17. 4 Connections between FHNC and parquet diagrams

18. 4.1 Rings

19. 4.2 Ladders

20. 4.3 Rungs

21. 4.4 Self-energy

22. 5 BCS theory for local correlations

23. 5.1 Weakly coupled systems

24. 5.2 Cluster expansions for a superfluid system

25. 5.3 FHNC and Euler equations

26. 5.4 Energy

27. 5.5 Uniform limit approximation

28. 5.6 Euler equation for the Bogoliubov amplitudes

29. 5.7 Long wavelength analysis

30. 6 Applicatios

31. 6.1 Neutron Matter

32. 6.1.1 Energetics

33. 6.1.2 Effective interactions

34. 6.1.3 Self energy

35. 6.1.4 BCS pairing

36. 6.2 The Lennard-Jones liquid

37. 6.2.1 Energetics and stability

38. 6.2.2 Effective interactions and correlations: A configuration space view

39. 6.2.3 Effective mass

40. 6.2.4 BCS pairing

41. 6.3 Pöschl-Teller interaction

42. 6.3.1 Energetics and stability

43. 6.3.2 BCS pairing

44. 7 Discussion

45. A Cluster expansions for the generating function. Diagonal terms

46. B Cluster expansions for the generating function. Off-diagonal terms

47. C Cluster expansions for the energy numerator terms

48. D Calculation of exchange diagrams

1 Introduction

The repertoire of methods for the quantitative microscopic description of normal quantum many-body systems has condensed, over the past few decades, to a relatively small number of techniques. These can be roughly classified on the one hand side as various shades of large-scale numerical simulation methods and, on the other hand, diagrammatic approaches summing, in some approximation, the parquet class of Feynman diagrams. Numerical simulations are capable of high precision but are computationally demanding and limited to relatively simple Hamiltonians and mostly ground state properties. They also need the physical intuition of the user about the possible state of the system. Semi-analytic, diagrammatic approaches lead to a better understanding of the underlying physical mechanisms, but have limited accuracy due to some necessary approximations. These diagrammatic approaches are versions of Green’s function methods [1], Coupled Cluster theory [2], and variational methods [3]. The interconnections between the various methods are well understood for Bose systems [4, 5]; the level and the details of implementation for different systems varies, however, vastly.

For normal systems, Jackson et al. make compelling arguments [4] that the summation of the so-called parquet diagrams is a minimum requirement for a microscopic treatment of the many-body problem that that treats both the short-ranged structure and the long-wavelengths excitations on equal footings. We will review these arguments further below. There are presently basically two theoretical approaches that have the diagrammatic completeness of the parquet diagrams, these are the Jastrow-Feenberg variational method [3] and the local parquet-diagram summation of Refs. 4 and 6. These methods have led, for boson systems, to exactly the same equations.

The situation is also intuitively clear for fermions, although technically more complicated due to the multitude of exchange diagrams generated by the antisymmetry of the fermion wave function, Unfortunately the fermion version [7] has so far not led to practical applications.

Among the many-body methods that sum, in some approximation, the parquet class of diagrams, the Jastrow-Feenberg method has been developed farthest. Both, local parquet theory and the Jastrow-Feenberg method are “robust” in the sense that exactly the same equations can be used for very different interactions like electrons, nucleons, and quantum fluids. Coupled Cluster theory [2] has also been very successful for electrons [8, 9] and nuclear systems [10, 11, 12] but it requires different truncation schemes for these two classes of many-body system. It lacks, therefore, the robustness of the Jastrow-Feenberg method. It was also less successful in predicting the ground state properties of the helium fluids. For bosons, a version of coupled cluster theory - the so-called “super-SUB-2” approximation has been developed [5] that is equivalent to the local parquet or Jastrow-Feenberg theory.

All of the above statements refer to normal systems. However, pairing phenomena are ubiquitous in the physics of many-body systems. Sixty years ago, the proof of the Cooper theorem [13, 14] provided the key to understanding the pairing phenomenon in interacting many-fermion systems, triggering the creation of the BCS theory of electronic superconductivity. It was quickly understood [15] that this work also has significant implications for nuclear phenomena, basically explaining the energy gap between the ground state and the first excited state for certain classes of nuclei. Since that landmark development, the basic paradigm of BCS theory has been extended to diverse fermionic systems with considerable success. We refer to two recent excellent review articles [16, 17] for a very complete account of the present situation and a comprehensive survey of the relevant literature. The work to be presented in this paper is meant to be complementary to these papers. We will, as far as justifiable, avoid overlaps. Rather, we shall focus on the technical aspects of many-body theory and spend very little space reviewing and comparing specific calculations.

A theory for superfluid many-body systems of the same diagrammatic completeness that was achieved for normal systems is presently unavailable, although specific partial summations of the perturbation series have been carried out [18, 19, 20, 21]. A version of Coupled Cluster theory for BCS-type wave functions has also been developed [22].

In systems where the pairing is due to the underlying many-body Hamiltonian, one often relies on effective interaction approximations which have the useful feature to permit the examination of mechanism and dependencies, but come with all the uncertainties involved in constructing effective interactions.

One of the intentions of our paper is therefore to develop a theory superfluid systems that is diagrammatically, as far as justified by the problems at hand, equivalent to the Jastrow-Feenberg or parquet theory for normal systems.

Our paper is organized as follows: In the following section, we will first give a pedagogical review of the motivations behind the correlated wave function method and the parquet-diagram theory. We will begin with the optimized hypernetted chain (HNC-EL) method for bosons which has proven to be the preferred systematic method of summing infinite classes of cluster diagrams.

As mentioned above, the boson version of the theory is identical to the summation of local parquet diagrams or to a specific version of coupled cluster theory. A priori, all two-particle vertices in a Feynman-diagram based theory are functions of four energy/momentum variables. Energy and momentum conservation and isotropy reduce that to ten variables which is still too much for diagram summation methods. Local parquet theory then introduces specific approximations to make all vertices functions of the momentum transfer only, and derives a procedure to replace the energy dependence by an average energy.

The situation is much more complicated for fermions for two reasons: One is that the antisymmetry requirement for the wave function leads to a multitude of exchange diagrams, the other is that the Fermi sea provides a natural frame of reference, whereas the Jastrow-Feenberg theory makes the approximation that all correlations depend only on the distance between particles. That requires futher approximations.

We will here examine what approximations must be made to go from a specific set of Feynman diagrams to a corresponding set of Jastrow-Feenberg diagrams. It will turn out that exactly the same procedure of defining an average energy that has been established for Bose systems [4] can be carried over to Fermi systems. The existence of a preferred frame of reference will require an additional Fermi-sea averaging procedure to generate vertices that depend only on the momentum transfer. We will see that exactly the same procedures for defining an average energy and an average momentum apply in the three channels, particle-hole, particle-particle, and single-particle propagators.

Realizing the relevant correspondences we will be led to formulate a hybrid theory that has the same diagrammatic content as the variational theory but avoids some important approximations.

We then develop the generalization of the Jastrow-Feenberg method for superfluids. We derive the diagrammatic expansions, carry out the FHNC summations for a correlated superfluid state, and derive the Euler equation for optimizing the correlations. By examining the Euler equation, we will unveil a severe problem of the Jastrow-Feenberg wave function: We will demonstrate that the Euler equation for the pair distribution function displays spurious instabilities for net-attractive interactions, characterized by an attractive Landau parameter $F_{0}^{s}<0$. We will demonstrate that these problems are caused by the so-called “collective” or “single-pole” approximation [26, 27] for the particle-hole propagator, which is implicit to the Jastrow-Feenberg wave function and will be discussed at length in Section 4. We are therefore lead to conclude that the Jastrow-Feenberg wave function for a superfluid system does not permit a sensible optimization of the pair correlations. Hence, one must go beyond the simple Jastrow-Feenberg method and to implement the correlated-basis functions (CBF) or parquet-diagram theory for superfluid states. We will show that the instabilities are then removed.

The need for developing CBF/parquet-diagram methods for the superfluid state has also a quantitative reason. We have argued  – and demonstrated – many years ago that local correlations of the kind (2.3) are inadequate to deal with pairing phenomena. The qualitative explanation for that is quite simple: local correlation functions treat all particles within the Fermi sea in the same way. BCS-pairing occurs at the Fermi surface; correlations that are independent of the location of the particle within the Fermi seas should therefore be particularly poor to describe pairing phenomena. Quantitative evidence for this was provided in our neutron matter calculations of Refs. 28 and 29. This is another reason that one must go beyond the Jastrow-Feenberg theory to deal with pairing phenomena reliably.

In Secs. 6.1, 6.2, and 6.3, we apply our theoretical methods to a few physically interesting cases: model Fermi gases at low densities, and neutron matter. With that, we follow up on previous works [30, 31] which was partly motivated by the interest in the BCS-BEC crossover in cold Fermi gases (see Refs. 32 and 33 for review articles) and superfluidity in neutron matter (see Ref. 34 for a collection of review papers). We have recently examined pairing phenomena in both model Fermi systems [35] and neutron matter [36]. We extend these calculations, which have assumed a small superfluid gap and treated the BCS correlations perturbatively, within the much more advanced theory to be developed in this paper.

In applications to neutron matter, we demonstrate that the inclusion of the full superfluid propagators in both the density and the spin-channel have a rather visible consequence for the superfluid gap.

The second case to be discussed are many-particle systems interacting via a family of Lennard-Jones model interaction. The attractive Lennard-Jones liquid has a more interesting phase diagram than neutron matter since it can have two spinodal points at which the speed of sound vanishes. One spinodal point appears at negative pressure at about 60 percent of the equilibrium density. This point can be reached by gradually lowering the density; it is characterized by the fact that the Fermi liquid Landau parameter $F_{0}^{s}\searrow-1$. A second spinodal point appears at very low density when the attractive interaction begins to dominate over the Pauli pressure. The appearance of this instability is obvious from the equation of state, it has already been observed by Owen [37].

We have studied BCS pairing for the Lennard-Jones interactions extensively in Ref. 35 in an approximation that assumed that the BCS correlations are weak. Going beyond this approximation, we face the aforementioned spurious instabilities of locally correlated wave functions which have a drastic effect in the Lennard-Jones liquid: In the whole density regime where the weakly coupled theory predicted a superfluid transition, FHNC-EL equations for local correlations have no physically acceptable solutions. We solve this problem by including the correct superfluid particle-hole propagator which has a quite visible quantitative effect on both the phase diagram and pairing properties.

Computations in the vicinity of the spinodal points become very demanding since the correlations become very long-ranged. In Ref. 35, we have observed that it is rather easy to come close to the upper spinodal point. However, getting close to the lower spinodal point turned out to be impossible since the solutions to the FHNC-EL equations diverge already a distance from the limit $F_{0}^{s}\searrow-1$. This divergence was identified as a divergence of the in-medium scattering length. We find exactly the same property in our much more advanced calculations to be presented in this paper. In fact, our inclusion of exchange diagrams, which improves the predictions of the Fermi-liquid parameter $F_{0}^{s}$ significantly, hardly changes the location of the instability.

Computations for the purely attractive Pöschl-Teller potential, which has due to the absence of a repulsive hard core no stable high-density phase and only the lower spinodal point confirm our conclusions.

We conclude this paper with a brief summary of our results.

2 Motivation: Variational and local parquet theory for bosons

2.1 Rationalization of parquet diagram summations

To describe the physics in the interaction-dominated short-distance region within diagrammatic perturbation theory, short-ranged correlations must be dealt with properly. These are treated, in perturbation theory, by summing the ladder diagrams which determine, among others, the pair distribution function $g(r)$ at small distances. The description of generally long-ranged effects, in particular phonons or plasmons and the behavior of the static structure function $S(q)$ at long wavelengths $q$, requires the summation of chain diagrams. The simultaneously correct treatment of both short- and long-ranged effects requires, therefore, the self-consistent summation of ring- and ladder-diagrams which defines the set of parquet diagrams [4, 6, 38].

The resulting two-body vertices are still functions of three independent four-momenta. To make execution of the theory practical, approximations must be made. Local parquet theory localizes the vertices by choosing an average energy such that the contribution to the pair distribution of the full vertex is the same as the that of the localized one. That way, the iterative procedure is, in every step, connected to a physical observable. Moreover, we shall see below that the pair distribution function can indeed be considered the only necessary independent variable that determines the properties of the system.

Carrying out this procedure for bosons, it turns out that one arrives at a set of equations that had been known for many years [3, 39, 40, 41], namely the HNC-EL equations [6] (dubbed “Paired-Phonon Analysis at that time) and the HNC-EL energy functional [38]. In fact, the analogy goes farther in the sense that the inclusion of the leading non-parquet diagrams [42] is the same as the inclusion of optimized three-body correlations in the wave function (2.3) [43].

The situation is more complicated for fermions, mostly due to the multitude of additional diagrams. Of course, since the correspondence between parquet diagrams and the HNC-EL method is clear for bosons, a similar correspondence should be expected for fermions.

In the nexy section we shall demonstrate the analogy between fermion parquet theory and FHNC-EL for important classes of diagrams, namely the ring- and ladder-diagrams. This is the essence of the simplest version of FHNC-EL, referred to FHNC-EL//0 [44]. More complete implementations [27] also include RPA-exchange diagrams, self-energy corrections, and mixture of all of these. These versions are referred to as FHNC//n where n is the level of higher-order exchange diagrams retained. As we shall see, the prescription of making all vertices energy-independent by choosing a well-defined average energy can be carried over from the boson parquet theory. The existence of a preferred reference frame, the Fermi sea, still causes the simplest vertices depend on three momenta; turning these into functions of momentum transfer will require the introduction of an additional specific averaging procedure of single-particle energies over the occupied states in the Fermi sea.

Both of these localizing procedures might seem ad-hoc from the point of view of conventional perturbation theory, and without further consideration other procedures might look equally well justified. They are rationalized by the fact that these localization prescriptions lead to the Jastrow-Feenberg wave function (2.3). The optimization prescription (2.8) makes sure that one has, that way, constructed the best wave function that can be represented in terms of local functions. In fact, a generalization of the Hohenberg-Kohn theorem to two-body functions [45, 46] shows that the pair distribution $g(r)$ is indeed quite generally defined by a variational problem.

The Jastrow-Feenberg wave function is known to reproduce the properties of both helium fluids with better than 90 percent accuracy [47]; below about 25 percent of the helium saturation density the accuracy of the FHNC//0 approximation is better than 1 percent, and Coulomb systems, both bosons [43] and fermions [48] are generally reproduced at the percent level or better.

2.2 Hypernetted Chain and Euler equations

We choose the Jastrow-Feenberg approach here to derive what we shall refer to as “generic” many-body method because it requires relatively little formal input. It is suitable for a non-relativistic many-body Hamiltonian

[TABLE]

The method starts with an ansatz for the $N$-body wave function,

[TABLE]

where {I_{\bf o}}=\bigl{\langle}{{\bf o}}\bigr{|}F_{N}^{\dagger}F_{N}\bigl{|}{{\bf o}}\bigr{\rangle} is the normalization constant. Here \bigl{|}{{\bf o}}\bigr{\rangle} is a model state, which is normally a Slater determinant for normal Fermi systems, and $F_{N}$ is an $N$-body correlation operator. The explicit reference to the particle number $N$ will be necessary later when we generalize the method to BCS states, we shall omit it from now on for the normal system which has a fixed particle number.

When truncated at the two-body term $u_{2}({\bf r}_{i},{\bf r}_{j})$, Eq. (2.3) defines the standard Jastrow theory. Historically [49] the theory was developed as a “quick and dirty” way to deal with the strong, short-ranged forces prevalent in nuclei. The task of the function $f(r)=\exp\left(\frac{1}{2}u_{2}(r)\right)$ is to bend the wave function to zero inside the regime of the repulsive hard core. The energy expectation value

[TABLE]

and other physically interesting quantities like the pair density

[TABLE]

the pair distribution function

[TABLE]

and the static structure function

[TABLE]

are then calculated by cluster expansion and resummation techniques. We will deal in this paper exclusively with translationally invariant and isotropic systems, hence $g({\bf r},{\bf r}^{\prime})=g(|{\bf r}-{\bf r}^{\prime}|)$.

A typical situation is seen in Fig. 1 where we examine the example of an interaction with a strong repulsive core and an attractive well. We show correlation and distribution functions for that interaction. The specific example is for one of the Lennard-Jones models to be discussed in Section 6.2. At higher densities, the pair distribution function begins to develop oscillations typical for strong short-ranged order. The correlation function $f^{2}(r)$ and the “dressed” correlation function $1+\Gamma_{\!\rm dd}(r)$ to be introduced in Section 3.3 describe the dynamic correlations induced by the interaction whereas the pair distribution function $g(r)$ is determined by both dynamic and statistical correlations.

The evaluation of physical quantities for the wave function (2.3) requires approximations. It was quickly realized [3] that the hypernetted chain summation (HNC) and truncation scheme has the advantage over other integral equation techniques like Born-Green-Yvon [50] or Percus Yevick[51] that it facilitates the unconstrained functional optimization of the pair correlations by solving the Euler equation

[TABLE]

To set the scene for the further discussions, and to demonstrate both the simplicity and the power of the method, as well as its physical content, let us first discuss the simpler case of a Bose liquid. In that case, the model state is \bigl{|}{{\bf o}}\bigr{\rangle}=1. The HNC scheme for bosons is known from the theory of imperfect gases where the Jastrow correlation function $u_{2}(r)$ is replaced by $-\beta v(r)$, $\beta$ being the inverse temperature [52, 54, 53]. The equations are

[TABLE]

Above, $X(r)$ and $N(r)$ are the sums of “non-nodal” and “nodal” diagrams, and $E(r)$ is the sum of “elementary” diagrams which can be expressed in terms of the pair distribution function $g(r)$. We define, as usual in this field, the dimensionless Fourier transform by including a density factor $\rho$:

[TABLE]

The elementary diagram contributions $E(r)$ have to be included term by term; they change the numerical values of the results, but not the analytic structure of the equations. Three-body correlations also lead only to a modification of that term [55].

The pair correlation function can be eliminated entirely from the theory by utilizing the Jackson-Feenberg identity

[TABLE]

For a Jastrow wave function, the second term can be written as

[TABLE]

and, eliminating $u_{2}(r)$ via the HNC equations (2.9), leads after a few manipulations to

[TABLE]

where $t(k)=\hbar^{2}k^{2}/2m$ is the kinetic energy.

It is then straightforward[39, 3, 56] to derive the Euler equation

[TABLE]

Skipping the technical details we display the resulting equations:

[TABLE]

Above,

[TABLE]

is the contribution from elementary diagrams and, if applicable, multiparticle correlations.

A few algebraic manipulations show that the pair distribution function satisfies the coordinate-space equation [56]

[TABLE]

Eq. (2.18) recognized as a Bogoliubov equation in terms of an effective “particle-hole” interaction $\tilde{V}_{\rm p-h}(k)$. Likewise, Eq. (2.22) is recognized as the Bethe-Goldstone equation in terms of the interaction $v(r)+\Delta V_{e}(r)+w_{\rm I}(r)$. This observation led Sim, Woo, and Buchler [40] to the conclusion that “it appears that the optimized Jastrow function is capable of summing all rings and ladders, and partially all other diagrams, to infinite order”.

The results (2.18) and (2.22) also substantiates the assertion made above that the HNC summation scheme is the method of choice over alternatives because it facilitates the optimization of the correlations: No matter which approximation we choose for the elementary diagrams, and whether we include higher order correlation functions, the only thing that changes is the correction term $\Delta V_{e}(r)$, but the structure of the equation remains the same.

2.3 Pair density functional theory: The view from the top

We have already commented in Section 2.2 that the pair correlation function $f(r)$ can be eliminated entirely from energy expression which is then formulated entirely in terms of the physical observable $g(r)$. It is therefore natural to ask whether a general minimum principle exists for the pair distribution function. Effectively, we are looking for a two-body version of the Hohenberg-Kohn [57, 58] theorem.

Let us write the energy per particle as

[TABLE]

where

[TABLE]

is the potential energy, and $T$ the kinetic energy whose form is yet unspecified.

Following the line of arguments leading to the Hohenberg-Kohn theorem for the one-body density, three statements can be made:

• (1)

The kinetic energy $T$ depends only on $g(r)$ and not on $v(r)$.

• (2)

Assuming that the interaction goes to zero at large distances, there is a bijective mapping between $v(r)$ and $g(r)$.

• (3)

The total energy has a minimum equal to the ground state energy at the physical ground state distribution function, in other words the ground state distribution function can be obtained by functionally minimizing the energy (2.23) with respect to the distribution function $g(r)$.

The proof parallels exactly the proof of the original Kohn-Hohenberg theorem and does not need to be repeated here.

Let us assume now that we have a variational problem of the form (2.17) with an energy functional (2.23), (2.24). We then can calculate the pair distribution function (or the static structure function) for any potential $\lambda v(r)$ with $0<\lambda<1$.

Replacing, in Eq. (2.24) $v(r)$ by $\lambda v(r)$ and differentiating with respect to $\lambda$ gives

[TABLE]

The second term in Eq. (2.25) vanishes, we can therefore recover the energy by coupling constant integration,

[TABLE]

where $g_{\lambda}(r)$ is the pair distribution function calculated for a potential strength $\lambda v(r)$, and $E_{0}$ is the energy of the non-interacting system which is zero for bosons, and equal to the energy $T_{\rm F}$ of the non-interacting Fermi gas for fermions.

In Eq. (2.26) we recover, of course, the Hellmann-Feynman theorem [59, 60] which was originally proven for the exact ground state. The above derivation [38] shows that the theorem is true not only for the exact ground state, but also for any approximate energy functional, as long as the pair distribution function is obtained by minimizing that functional.

Evidently, the above statement is much more general than the optimization condition (2.17) for the Jastrow–Feenberg wave function because it defines a whole class of many–body theories which can be characterized by the central role of the pair distribution function. To summarize, the above consideration shows that the many-body theory of strongly interacting systems can indeed quite generally be formulated in terms of a local two-body quantity. This provides, similar to the tremendously successful density functional theory of inhomogeneous electron systems, a significant simplification compared to Greens’s function theories.

The above analysis is evidently independent of the statistics, it applies equally well for bosons and fermions.

So far, our considerations were entirely parallel to conventional density functional theory; the energy functional is still unspecified. The next step is therefore the construction of a “pair density functional”. Unlike the conventional density functional of inhomogeneous electron systems, some exact features of the pair distribution function are known that can be used for the construction of a pair-density functional (2.23).

The static structure function $S(q)$ is related to the density–density response function $\chi(q,\omega)$ through

[TABLE]

Assuming, for example, an RPA form of the density-density response function defines a local “particle–hole interaction” $\tilde{V}_{\rm p-h}(q)$ by an RPA formula

[TABLE]

where $\chi_{0}(q,\omega)$ is the Lindhard function. 2. 2.

The short-ranged structure of the pair distribution function $g(r)$ is determined by a Bethe-Goldstone equation in terms of a yet unspecified particle–particle interaction $V_{\rm p-p}(r)$ We can assume that, for short, distances, $V_{\rm p-p}(r)$ is dominated by the bare interaction $v(r)$, i.e. we can write

[TABLE] 3. 3.

We require that $g(r)$ and $S(k)$ are consistent in the sense that they are related by Eq. (2.7).

For bosons, all the calculations can be carried out analytically. In the absence of Pauli operators, the Bethe Gladstone equation is simply a zero-energy Schrödinger equation

[TABLE]

Just as we can think of Eqs. (2.27) and (2.28) as a definition of $\tilde{V}_{\rm p-h}(q)$ in terms of $S(q)$, we can think of Eq. (2.29) as a definition of $V_{\rm p-p}(r)$ in terms of $g(r)$.

For bosons we have

[TABLE]

Then, the frequency integration (2.27) can be carried out analytically and leads to the familiar Bogoliubov formula (2.18). Simply manipulating Eqs. (2.27)-(2.30) then leads to the HNC-EL equations (2.18)-(2.20), (2.22) [45]. The only undetermined quantity is the correction $V_{e}(r)$ which is, in HNC-EL, determined by the elementary diagrams and multiparticle correlation functions retained; in local parquet theory it is the set of diagrams that are both particle-particle and particle-hole irreducible [42].

To summarize, the HNC-EL theory supplements the variational prescription following from the Hohenberg-Kohn theorem for the pair distribution function by the requirement that $g(r)$ function satisfies both, a Bethe-Goldstone equation and an RPA equation. Hence, we shall refer to the equations as generic equations because they can be obtained without ever mentioning a Jastrow-Feenberg function.

2.4 Stability and Consistency

A condition for the existence of solutions of the Euler equation is that the term under the square-root in Eq. (2.18) is positive. We must identify the long-wavelength limit with the hydrodynamic speed of sound,

[TABLE]

to obtain the correct long wavelength limit

[TABLE]

An immediate consequence is that the HNC-EL or local parquet equations have no solution of the system is unstable against infinitesimal density fluctuations. This is a very desirable feature and unique to theories that have the diagrammatic completeness of the parquet theory.

Alternatively we can calculate the hydrodynamic speed of sound from the equation of state

[TABLE]

The definitions (2.31) and (2.33) will, in any approximate theory, not be identical. In fact, it can be shown in both Jastrow-Feenberg theory [61] and in parquet-diagram theory [62] that they agree only when all diagrams and correlations to all orders are included. Turning this feature in an advantage, the comparison between the definitions (2.31) and (2.33) can serve as a convergence test of approximate evaluations. We will utilize this feature in our numerical studies below.

3 Variational and local parquet diagram theory for fermions

The Jastrow-Feenberg method for fermions has been successfully applied to the relatively simple electron liquid [63, 64, 48], nuclear systems [65] as well as highly correlated Fermi systems like 4He [55, 46] and 3He [37, 61, 27] at $T=0$ and finite temperatures [66, 67, 68]. The full fermion HNC equations [69, 70, 71] are significantly more complicated than the bosons equations (2.9); instead of one set of nocal, non-nodal, and elementary diagrams we have four sets. Moreover, very specific truncation schemes of exchange diagrams are necessary to permit an unconstrained optimization of the correlations [61], the above-mentioned hierarchy of FHNC//n approximations.

We have shown in recent work [72] that even the simplest version of the FHNC-EL theory reproduces the equation of state within better than one percent at densities less than 25% of the saturation density of liquid 3He. This statement applies to the energy, other quantities are, as we shall see, more senesitive to level at which the FHNC are implemented. A similar statement applies for nuclear systems [36]. It is not much more complicated to solve the full set of FHNC-EL equations, including elementary diagrams and triplet correlations [27]. The version to be presented here permits, however, a clearer identification of sets of FHNC diagrams with Feynman diagrams.

3.1 Generating functional and the generalized distribution functions

We describe in this and the following sections the basic techniques of cluster expansions and resummations for Fermi systems. The manipulations of Section 2.2 relied on the simplicity of the HNC equations (2.9) for bosons. We must now be more systematic, this is also necessary in view of the generalization to superfluid systems to be described below.

The central quantity for all derivations is the “generating function” $G(\beta)$ defined as follows: Substitute in the variational wave function (2.3)

[TABLE]

where

[TABLE]

is the “Jackson-Feenberg effective interaction”. With this, all quantities depend on the dummy parameter $\beta$. Define then the generalized normalization integral

[TABLE]

and the generating function

[TABLE]

The value $G(0)$ is just the logarithm of the norm of the wave function, but evidently it is technically no more complicated to calculate $G(\beta)$ than it is to calculate $G(0)$.

The energy expectation value may be obtained from the generating functional through

[TABLE]

Here $T_{\rm F}$ is the kinetic energy of the non-interacting Fermi gas, and

[TABLE]

is an energy correction in Fermi systems that will be dealt with below. Further, it is convenient to define generalized densities

[TABLE]

The (generalized) two-body density may be derived as a variational derivative of the generating function with respect to the (generalized) pair correlation function

[TABLE]

a relation that will be particularly useful for the derivation of the Euler equations. The last expression shows that the procedure generates a $g({\bf r},{\bf r}^{\prime})$ which has an overall factor $f^{2}({\bf r},{\bf r}^{\prime})=\exp(u_{2}({\bf r},{\bf r}^{\prime})$.

For the case $\beta=0$, the definition (3.7) reduces to the two-body density introduced above (cf.  Eq. (2.5)),

[TABLE]

The construction of both the energy and the distribution functions via derivatives of one common generating functional is a welcome economy since it is sufficient to develop an algorithm for the calculation of $G(\beta)$; one does not need to start over for each individual quantity of interest. For example, we can immediately write down the exact Euler equation for the pair distribution function:

[TABLE]

where

[TABLE]

3.2 Cluster expansions

A very straightforward technique to derive cluster expansions for the generating function and, or course, any other quantity of interest is the power-series method [71]. The method introduces a dummy parameter $\alpha$ in the correlation operators

[TABLE]

where $h(r)\equiv f^{2}(r)-1$ and we have suppressed the parameter $\beta$. One then expands the quantity of interest in a power series in $\alpha$, evaluated at $\alpha=1$. For example,

[TABLE]

That way, a series of cluster contributions $(\Delta G)^{(n)}(\beta)$ of increasing number of $h(r_{ij})$ factors is generated. These are best represented diagrammatically [44]:

Each point (open dot) represents a particle coordinate ${\bf r}_{i}$. 2. 2.

Each filled point (solid dot) implies the integration over the coordinate ${\bf r}_{i}$, and multiplication with a density factor $\rho=N/\Omega$, where $\Omega$ is the normalization volume. 3. 3.

Correlation line connecting the points ${\bf r}_{i}$ and ${\bf r}_{j}$ represent a function $h(r_{ij})$. These are depicted as dashed lines connecting the two points.

[TABLE] 4. 4.

Exchange lines represent the function

[TABLE]

These are depicted by oriented solid line connecting the point ${\bf r}_{i}$ to ${\bf r}_{j}$,

[TABLE]

We call a diagram linked if each point is connected to every other point by at least one continuous path of graphical elements. Linked diagram contributions to the generating functional are proportional to the particle number. We call diagram irreducible if it cannot be calculated as a product of two or more simpler diagrams.

The cluster expansion of the generating function is then represented in terms of all topologically distinct irreducible diagrams without external points constructed by the following rules:

Each point is attached by at least one correlation line. 2. 2.

Two points can be joined by at most one correlation line. 3. 3.

Any point of a contributing diagram is joined by at most one incoming exchange line which must be accompanied by a single outgoing exchange line. Hence, exchange lines come in closed loops. For each closed loop of $L$ exchange lines there is a factor $(-1/\nu)^{L-1}$ in the corresponding analytic contribution, as well as the exchange function. (Generally, the sign $(-)^{L-1}$ is displayed with the diagram, but the numerical factor $1/\nu^{L-1}$ is left implicit.) Here, $\nu$ is the degree of degeneracy of the single particle states.

The first terms in the diagrammatic expansion of the generating function are shown in Fig. 2.

The generalized pair distribution function can then be derived by the variational prescription (3.8) and the energy expectation is obtained from the general expression (3.5) as follows

Replace, in turn, each correlation line $h(r_{ij})$ by a line $f^{2}(r_{ij})v_{\rm JF}(r_{ij})$. This is what the $\beta$-derivative does. 2. 2.

To calculate the kinetic energy corrections $T_{\rm JF}$ replace any pair of exchange lines $\ell(r_{ij}k_{\rm F})\ell(r_{ik}k_{\rm F})$ by a pair $\frac{\hbar^{2}}{8m}\nabla_{i}^{2}\ell(r_{ij}k_{\rm F})\ell(r_{ik}k_{\rm F})$.

Cluster contributions for the (generalized) pair distribution function can then be generated by the prescription (3.8). Diagrammatically, the construction of $g(r,\beta)$ amounts to the following operation on the generating function $G(\beta)$:

Remove, in turn, each correlation line and turn its endpoints into open (“reference”) points ${\bf r}$ and ${\bf r}^{\prime}$ 2. 2.

Multiply this by $\exp(u_{2}({\bf r},{\bf r}^{\prime}))=1+h({\bf r},{\bf r}^{\prime})$. This factor is needed to model the short-ranged structure of the wave function for hard-core potentials.

A few low-order diagrams contributing to the pair distribution function are shown in Fig. 3.

The figure shows the expansion of $g(r)$ in terms of correlation and exchange lines. Two adjacent diagrams (a) and (b), (c) and (d) etc. can always be combined to obtain the form (3.9), in other words the function $C(r)$ is represented by diagrams (a), (c), (e), etc.. This is exactly the intention of Jastrow-Feenberg theory [49, 3], namely to have the correlation function model the short-ranged structure of the wave function.

Close inspection of the individual contributions to $g(r)$ reveals, however, a dilemma: For the stability of the system as well as for meaningful solutions of the Euler equation, it is important long ranged correlations are important, cf. Eq. (2.32). To get the correct behavior of $S(q)$ for $q\rightarrow 0+$, the diagrams must be grouped differently. For example, diagrams (b), (h), and (k) combine to

[TABLE]

where

[TABLE]

is the static structure function of the non-interacting Fermi gas. Likewise, it is straightforward to show [61] that the sum of diagrams (g) and (o) as well as the sum of diagrams (d), (i), and (m) go as $q^{2}$ in the limit $q\rightarrow 0+$. The analysis can be easily extended to other, more complicated cases.

To summarize, there is no finite truncation of the expansion of the pair distribution function that is exact in both, the short-distance and the long-wavelength limit. One can deal with this situation in three ways:

One can ignore the problem entirely and use, consistent with the original idea of Jastrow-Feenberg theory, an approximation for $g(r)$ of the form (3.9). This is, among others, the idea of the FHNC summations of Ref. 71. One must live with the fact that the static structure function has the incorrect long-wavelength limit and the correlation functions are limited to simple parameterized forms.

• 2.

One can use different approximations for $g(r)$ and $S(q)$ depending on which quantity is of interest. Choosing a form for $S(q)$ that has the exact long-wavelength behavior permits the unconstrained optimization of the pair correlations. Of course, one must then construct a pair distribution $g(r)$ of the form (3.9).

• 3.

One must sum infinite sets of exchange diagrams or approximations thereof as was done in Ref. 27.

For the purpose of this work, we shall use the second approach because it is sufficiently accurate for all of our purposes [72, 36] and allows the most direct identification of JF diagrams with parquet diagrams.

To conclude this section, we mention another specific set of diagrams, the so–called “cyclic chain” (cc-) diagrams. These are two body diagrams that have a continuous exchange path connecting the external points. A few examples are shown in Fig. 4. We can again identify “chain” diagrams (diagrams (c), (d), and (e)), “parallel connections” (diagrams (a) and (b)) and one “elementary” diagram (diagram (f)).

The summation of these diagrams suggests the introduction of a “dressed” exchange line $L(r)$ which is the sum of all 2-point diagrams having an exchange path connecting the two external points. The “cc” diagrams are related to self-energy corrections in perturbation theory, see Sections 3.9 and 4.4 below.

3.3 FHNC and Euler equations

We discuss here the simplest implementation of the FHNC theory that is compatible with the variational problem, called FHNC//0 approximation. The approximation keeps in $g(r)$ only those diagrams that contain exchange loops of the form $\ell^{2}(r_{ij}k_{\rm F})$. These are, for example, the diagrams (a-g), (k) and (l) shown in Fig. 3. We must also keep in mind that a different approximation has to be used in $S(q)$. The implementation and relevance of higher order exchange corrections will be discussed below in Section 3.7.

In the FHNC//0 approximation, the FHNC equations are no more complicated than the Bose HNC equations (2.9):

[TABLE]

Above, $\Gamma_{\!\rm dd}(r)$ is the set of all diagrams contributing to $g(r)$ that have no exchange lines attached to the external points. Examples are diagrams (b), (e), (f), (k) and (l) shown in Fig. 3. We need these equations to eliminate $u_{2}(r)$ from the $v_{\rm JF}(r)$ and to formulate the theory entirely in terms of the physically observable static structure function $S(k)$ and derived quantities.

In this approximation, the energy correction $T_{\rm JF}$ has the form

[TABLE]

Similar to the Bose case, the Euler equation is best formulated in momentum space. Formally, the optimization condition‘ can be written as [3]

[TABLE]

The $S^{\prime}(k)$ is generated by the $\beta$-derivative technique outlined above. This amounts to

Replace, in turn, each correlation line by $f^{2}(r)v_{\rm JF}(r)$. This can be done with the $\beta$-derivative trick. 2. 2.

To add the correction from $T_{\rm JF}$, replace, in turn, each $S_{\rm F}(k)$ by $-\frac{t(k)}{2}\left[S_{\rm F}(k)-1\right]$.

Using the version (3.20), (3.22) of the FHNC00 equations as well as the “priming” operation outlined above we obtain

[TABLE]

and, inserting this in the Euler equation (3.24) and solving for $S(k)$:

[TABLE]

In the FHNC//0 approximation, the effective interaction $\tilde{V}_{\rm p-h}(k)$ is approximated by the “direct” part of the particle–hole interaction, $\tilde{V}_{\rm p-h}(k)\approx\tilde{V}_{\rm dd}(k)$; we will discuss the importance to exchange corrections further below. The quantity is structurally identical to the one for bosons:

[TABLE]

where

[TABLE]

is the “Clark-Westhaus” effective interaction [44] and $w_{\rm I}(r)$ is the “induced interaction”.

[TABLE]

The second line is obtained by using Eq. (3.26) to eliminate $\tilde{V}_{\rm p-h}(k)$. The Bose limit is obtained by replacing $S_{\rm F}(k)\rightarrow 1$.

To derive the equation determining the short-ranged structure of the correlations, begin with Eq. (3.30) which we can write, using the Euler equation (3.26), as (cf. Eq. (2.62) of Ref. 73)

[TABLE]

Using the identity

[TABLE]

Eq. (3.31) becomes

[TABLE]

The right-hand side is evidently zero for bosons, and the Euler equation is a simple zero-energy Schrödinger equation where the bare interaction is supplemented by the induced potential which guarantees that the scattering length of the effective interaction $v(r)+w_{\rm I}(r)$ is zero. This is the well-known result of Refs. 56 and 6. For fermions, the right hand side changes the short-ranged behavior of the correlation function $\Gamma_{\!\rm dd}(r)$ and, hence, the short-ranged behavior of the pair distribution function $g(r)$. This is consistent with the prediction of the Bethe-Goldstone equation that the Pauli principle changes the short-ranged behavior of the wave function [74], see Eq. (4.22) below.

3.4 Energy

For $E_{\rm FHNC}$ we use Eq. (2.14) of Ref. 35 which implies the form (3.23) for the $T_{\rm JF}$. Summarizing, we use

[TABLE]

where $T_{\rm F}$ is the kinetic energy of the non-interacting Fermi gas, theory, and the pair distribution function is given by [35]

[TABLE]

Above,

[TABLE]

The two terms $(\Delta\tilde{\tilde{X}}_{\rm ee})_{1}^{(3)}(k)$ and $(\Delta\tilde{X}_{\rm ee})_{1}^{(4)}(k)$ are diagrammatically represented bu the 3-body and 4-body diagram shown in Fig. 5

That term is omitted the FHNC//0 approximation, it is for energy calculations in our case very small and only needed to establish the exact low density expansion of the energy in powers of $a_{0}k_{\rm F}$ to second order [75]. We found that non-universal contributions to the equation of state are overwhelming well below the densities where the second-order terms become visible. However, we will see that the correction to the particle-hole interaction originating from these diagrams is substantial even in the low-density limit.

3.5 Uniform limit approximation

The so-called “uniform limit” approximation [3] is obtained by assuming a weak, but possibly long–ranged interaction. We study this because it will permit direct contact to be made to the random phase approximation. It will also highlight a special feature for the superfluid system, see 5.5 below.

Specifically, the uniform limit approximation is valid when

[TABLE]

note that we do not assume that $\tilde{v}(k)\tilde{\Gamma}_{\!\rm dd}(k)$ or $\Gamma^{2}_{\!\rm dd}(k)$ are negligible.

In the energy term $E_{\rm R}$ in Eq. (3.34) we can use

[TABLE]

This term is not negligible because we can write the energy contribution as

[TABLE]

Then, the energy becomes

[TABLE]

where $E_{\rm HF}$ is the energy expectation value in Hartree-Fock approximation. The Euler equation (3.38) becomes

[TABLE]

follows from this expression by minimization with respect to $S(k)$. We will see below that both (3.37) and (3.38) are is indeed approximations for of the RPA energy and structure function. Since the $S(k)$ is obtained from a variational principle, it follows from our analysis of Section 2.3 that the energy (3.37) can be obtained by coupling constant integration which is now, of course, best performed in momentum space.

3.6 The low-density limit

Cold gas applications focus on very low density systems where information about the interaction can be reduced to a single parameter, the $S$-wave scattering length $a_{0}$. We have studied this area very carefully in Ref. 35, see also sec. 3.8. We have identified three areas where the correlations are determined by different effects:

The short-distance regime of the order of the interaction range is, of course, dominated by the shape of the interaction. 2. 2.

The intermediate regime is the range between the interaction regime and the average particle distance. If this regime is large, which is the case for very low–density systems, the correlations are determined by the vacuum scattering Knight $a_{0}$,

[TABLE] 3. 3.

For interparticle distances larger that $1/k_{\rm F}$, the correlations fall off as

[TABLE]

where $mc_{\rm F}^{2}=\hbar^{2}k_{\rm F}^{2}/3m$ is the incompressibility of a non-interacting Fermi gas. Evidently this is a many-body effect, the fall-off of the long-range correlations is needed to have the wave function normalized.

A manifestly microscopic calculation begins, of course, with the bare interaction which determines the vacuum scattering length. Since all diagrammatic calculations imply some approximations, we need to make sure that any approximations we are making will not destroy this property. This will turn out to be an important consideration in the next section.

To derive the low-density limit of the Euler equation, begin with the diagrammatic expansion of the pair distribution function shown in Fig. 3. In that limit, only the first line of diagrams contributes; moreover we can identify $\Gamma_{\rm dd}(r)=f^{2}(r)-1$. Then the Euler equation reads

[TABLE]

where $g_{\rm F}(r)=1-\frac{1}{\nu}\ell^{2}(rk_{\rm F})$ is the pair distribution function of the non-interacting system, and $\Gamma_{\!\rm dd}^{\prime}(r)$ is constructed in analogy to $S^{\prime}(k)$ above which yields

[TABLE]

We can cancel the term $\frac{\hbar^{2}}{4m}\nabla^{2}g_{\rm F}(r)$ on the left and the right; the term $\nabla g_{\rm F}(r)\cdot\nabla\Gamma_{\rm dd}(r)$ goes as $k_{\rm F}^{2}$ and can be ignored; hence we end up with

[TABLE]

which is exactly the bosons form. Going through the same manipulations as in Section 3.3 we finally can write the low-density limit of the Euler equation as

[TABLE]

Evidently, this is identical to the low-density limit of Eq. (3.33) saying that the FHNC//0 approximation as spelled out in the previous section describes both the short– and the long– ranged correlations consistently. Of course, the simplicity of this equation is caused by the fact that the Pauli principle plays no role in the low-density limit, see Eq. (3.33).

3.7 Exchange corrections

We have above formulated a version of FHNC-EL that contains the simplest versions of the RPA and the Bethe-Goldstone equation. These describe the qualitatively correct physics, but have, even at very low densities, some quantitative inconsistencies which we have to address and handle.

Let us go back to the energy expression (3.34). The dominating term at low densities is

[TABLE]

Taking this Hartree-Fock like expression and ignoring the density dependence of the correlation functions, we get

[TABLE]

which evidently differs from

[TABLE]

by approximately a factor of $1-1/\nu$.

That means, exchange diagrams must be included to obtain consistency between the hydrodynamic derivative (2.33) and $\tilde{V}_{\rm p-h}(0+)$ in the low-density limit. The simplest version of the FHNC hierarchy that corrects for this deficiency is FHNC//1 which includes the exchange diagrams shown in Fig. 5. We can extract the relevant modification from the full FHNC-EL equations as formulated in Ref. 27 by keeping only the exchange term $V_{\rm ee}(k)$. The Euler equation remains practically the same, except that the static structure function of the non-interacting system becomes

[TABLE]

where $S_{\sigma}(k)$ is the spin-structure function obtained from the wave function (2.3), and the combination $S_{\rm F}(k)+(\Delta\tilde{X}_{\rm ee})^{(1)}(k)$ is the two-body approximation for that function. Note that $(\Delta\tilde{X}_{\rm ee})^{(1)}(k)\propto k^{2}$ as $k\rightarrow 0+$ whereas $S_{\rm F}(k)\propto k$.

The particle-hole interaction is modified by

[TABLE]

where $\tilde{V}_{\rm ee}(k)$ is represented by the sum of the three diagrams shown in Fig. 5. The red wavy line must then be identified with

[TABLE]

which is, of course, in the low density limit equal to $V_{\rm CW}(r)$. The Euler equation becomes

[TABLE]

The induced interaction is also modified and has, in the form of Eq. (3.29) an additional term

[TABLE]

Skipping the technical details (see also D), the long-wavelength expansion of $\tilde{V}_{\rm ee}(k)$ is found to be

[TABLE]

The factor $4/3$ compared to Eq. (3.45) is an incorrect consequence of local correlation functions, but that term vanishes in the low-density limit. The leading term in the density expansion comes out correctly if the first order exchange diagrams are included. Moreover, $W(r)\approx V_{\rm p-h}(r)$, and we see that the first order exchange corrections lead to the desired factor $1-1/\nu$ as they should.

However, the naïve addition of exchange diagrams is problematic in the limit of low densities. In that limit, the three and four-body diagram in $\tilde{V}_{\rm ex}(k)$ can be neglected, and we have

[TABLE]

Following the derivations of Section 3.3 we and up with a coordinate space equation of the form

[TABLE]

The last line in Eq. (3.54) goes to zero in the low density limit, but the term in the second line does not. This leads to solutions that are very different from the vacuum solution derived in the previous section. The only way to rectify this situation (short of solving the full FHNC-EL) equations is to use a slightly modified relationship between $S(k)$ and the effective interactions, namely

[TABLE]

This relationship has a number of very interesting and very desirable features: First, the square-root term in the numerator may be identified with a “collective RPA” expression for the spin-structure function,

[TABLE]

we refer the reader to Section 4.1 for a justification. The expression (3.50) is then obtained by expanding $S_{\sigma}(k)$ to first order in the interactions and identifying

[TABLE]

As mentioned above, the FHNC//1 approximation only leads to the two-body approximation which is unsatisfactory for another reason: As stated above, the spin-static structure function in that approximation goes, for small $k$, as $S_{\rm F}(k)$ which disagrees with experiments. The expression (3.56) does not have this problem. In other words, the thorough examination of the variational problem and the demand for consistent treatment of short– and long ranged correlations as well as the proper low-density limit provides crucial information on adequate approximation schemes for the Euler equation.

We have commented above on the fact that the positivity of the term under the square root in the denominator is, with the qualification that the Jastrow-Feenberg wave function is not exact, related to the stability against density fluctuations. Likewise, the positivity of the numerator is connected with the stability against spin-density fluctuations.

In perturbation theory, the diagrams shown in Fig. 5 correspond to the particle-hole ladder diagrams, driven by the exchange term of the particle-hole interaction

[TABLE]

This non-local term leads to a rather complicated addition to the summation of the ring diagrams in the sense that it would supplement the RPA sum by the RPA-exchange (or particle-hole ladder) summation. We can again make the connection to the (local) FHNC expression (3.48) by realizing that this expression is obtained from the exact expression (3.57) by exactly the same hole-state averaging process that was discussed in Section 4.1, Eq. (4.3):

[TABLE]

The discussion of the preceding section shows that this localization procedure of the exchange term maintains the dominant part of the long-wavelength limit.

3.8 Limitations of local correlation functions

One often learns most about a theory by examining situations where it fails. The wave function (2.3) is in principle exact for bosons when correlation functions to all orders are included. It is not exact for fermions since the nodes of the wave function are identical to those of the non-interacting system.

A first consequence of the limitations of local correlation functions was pointed out by Zabolitzky [64]: The exact high-density expansion of the equation of state of a homogeneous electron gas is, in units of the Wigner-Seitz radius $r_{s}$ [76, 77].

[TABLE]

A wave function of the form (2.3) leads, in this case, to 0.05690$\,\ln r_{s}\,$Ry for the logarithmic term. This deficiency is corrected [78] by second order correlated basis functions (CBF) theory which will be reviewed in the next Section 3.9

More recently [35], we have examined the low-density limit of equation of state of the weakly interacting Fermi gas. The exact limit is expressed as a power series expansion in terms of the parameter $a_{0}k_{\rm F}$ [75, 79]

[TABLE]

With the wave function (2.3) one obtains for the coefficient of the third term the result of $1.5415$ instead of the exact value $4(11-2\ln 2)/35=1.098$. Again, second order CBF theory corrects this limit [35].

A third issue is the stability of the system against infinitesimal density fluctuations. We have commented about the connection between the long-wavelength limit $V_{\rm p-h}(0+)$ and the hydrodynamic compressibility for bosons. In a Fermi fluid, we also have the Pauli repulsion, i.e. we should identify

[TABLE]

where $c_{\rm F}^{*2}=\frac{\hbar^{2}k_{\rm F}^{2}}{3mm^{*}}$ is the speed of sound of the non-interacting Fermi gas with the effective mass $m^{*}$, and $F_{0}^{s}$ is Landau’s Fermi liquid parameter. Requiring a positive compressibility leads to Landau’s stability condition $F_{0}^{s}>-1$.

Solutions of the FHNC-EL equation exist if the expression under the square root of Eq. (3.26) is positive, which leads to the stability condition $F_{0}^{s}>-4/3$. This result is again a manifestation of the fact that the wave function (2.3) is not exact. We shall return to this issue when we discuss the generalization of the Euler equation to superfluid systems.

A final word is in order on the identification of $\tilde{V}_{\!\rm p-h}(0+)$ with the Fermi-liquid parameter $F_{0}^{s}$. We have already seen above that exchange corrections are important and improve the agreement in leading order in the exchange corrections. More generally, the identification is, of course, only approximate and, just like for bosons, exact agreement can be achieved only in an exact theory [62]. A diagrammatic analysis that makes no assumptions on the level of FHNC theory used is found in Refs. 61, 80. The FHNC//0 and FHNC//1 versions makes more approximations, we have discussed these in Section 3.7.

3.9 Elements of correlated basis functions

We have seen above that a locally correlated wave function (2.3) fails to reproduce several known exact features of many-particle systems. A way to cure the problem is provided by perturbation theory with correlated basis functions (CBF theory) [3, 81]. Some of the basic ingredients of CBF theory are also required for the formulating a theory for weakly coupled superfluid systems [30, 31]. We review CBF theory here only very briefly for the purpose of defining the essential ingredients, details may be found in pedagogical material [25] and review articles [44, 27]. The diagrammatic construction of the relevant ingredients has been derived in Ref. 23.

CBF theory uses the correlation operator $F$ to generate a complete set of correlated and normalized $N$-particle basis states through where the correlated states $N$-body \bigl{|}{\Psi_{\bf m}^{(N)}}\bigr{\rangle} are

[TABLE]

where the \{\bigl{|}{{\bf m}^{(N)}}\bigr{\rangle}\} form complete sets of $N$-particle states, normally Slater determinants of single particle orbitals. When unambiguous, we will omit the superscript $N$ referring to the particle number. Although the \bigl{|}{\Psi_{\bf m}^{(N)}}\bigr{\rangle} are not orthogonal, perturbation theory can be formulated in terms of these states [82, 3].

In general, we label “hole” states which are occupied in \bigl{|}{{\bf o}}\bigr{\rangle} by $h$, $h^{\prime}$, $h_{i}\;,\ldots\,$, and unoccupied “particle” states by $p$, $p^{\prime}$, $p_{i}\;,$ etc.. To display the particle-hole pairs explicitly, we will alternatively to the notation \bigl{|}{\Psi_{\bf m}}\bigr{\rangle} use \bigl{|}{\Psi_{p_{1}\ldots p_{d}\,h_{1}\ldots h_{d}}}\bigr{\rangle}. A basis state with $d$ particle-hole pairs is then

[TABLE]

For the off-diagonal elements $O_{\bf m,n}$ of an operator $O$, we sort the quantum numbers $m_{i}$ and $n_{i}$ such that \bigl{|}{{\bf m}}\bigr{\rangle} is mapped onto \bigl{|}{{\bf n}}\bigr{\rangle} by

[TABLE]

From this we recognize that, to leading order in the particle number $N$, any matrix element of an operator $\hat{O}$

[TABLE]

depends only on the difference between the states \bigl{|}{{\bf m}}\bigr{\rangle} and \bigl{|}{{\bf n}}\bigr{\rangle}, and not on the states as a whole. Consequently, $O_{\bf m,n}$ can be written as matrix element of a $d$-body operator

[TABLE]

(The index $a$ indicates antisymmetrization.)

The key quantities for the execution of the theory are diagonal and off-diagonal matrix elements of unity and $H\!-\!H_{\bf o}$,

[TABLE]

Eq. (3.68) defines a natural decomposition [23, 25] of the matrix elements of $H$ into the off-diagonal quantities $W_{\bf m,n}$ and $N_{\bf m,n}$ and diagonal quantities $H_{\bf m}$.

To leading order in the particle number, the diagonal matrix elements of $H\!-\!H_{\bf o}$ become additive, so that for the above $d$-pair state we can define the CBF single particle energies

[TABLE]

with $e_{ph}=e_{p}-e_{h}$ where

[TABLE]

and $u(p)$ is an average field that can be expressed in terms of the compound diagrammatic quantities of FHNC theory [23]

According to (3.66), $W_{{\bf m},{\bf n}}$ and $N_{{\bf m},{\bf n}}$ define $d-$particle operators ${\cal N}$ and ${\cal W}$, e.g.

[TABLE]

Diagrammatic representations of ${\cal N}(1,2,\ldots,d)$ and ${\cal W}(1,2,\ldots,d)$ have the same topology [23]. In homogeneous systems, the continuous parts of the $p_{i},h_{i}$ are wave numbers ${\bf p}_{i},{\bf h}_{i}$; we abbreviate their difference as ${\bf q}_{i}$.

In principle, the ${\cal N}(1,2,\ldots,d)$ and ${\cal W}(1,2,\ldots,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. 5.3. In the same approximation, the operators ${\cal N}(1,2)$ and ${\cal W}(1,2)$ are local, and we have [27]

[TABLE]

Explicit formulas for the single-particle spectrum $e_{k}$ may be found in Ref. 27. Since we are only interested in relatively weakly correlated systems, and also want to make the connection to perturbation theory as transparent as possible, we spell out the simplest form:

[TABLE]

where

[TABLE]

Note that we have above approximated, among others, $L(r)\approx\ell(rk_{\rm F})$.

One of the most straightforward application of CBF theory is to calculate corrections to the ground state energy. In second order we have, for example,

[TABLE]

The magnitude of the CBF correction is normally comparable to the correction from three-body correlations [27]. It is also important to note that there are significant cancellations between the two terms in the numerator. We have shown in previous work that including the second order CBF correction (3.76) leads in the electron gas [78] to the correct expansion (3.59) as well to Huang-Yang expansion (3.60) for a low-density Fermi gas [35].

4 Connections between FHNC and parquet diagrams

4.1 Rings

The expression (3.26) reduces to the Bogoliubov equation for the case of bosons, $S_{\rm F}(k)=1$. It is therefore expected that it can also be derived, for fermions, from the random phase approximation for the dynamic structure function

[TABLE]

where $\chi_{0}(k,\omega)$ is the Lindhard function, and $\tilde{V}_{\rm p-h}(k)$ is a local quasiparticle interaction or “pseudopotential” [83, 84]. Consistent with the convention (2.10) according to which $\tilde{V}_{\rm p-h}(k)$ has the dimension of an energy, we have defined the density-density response function slightly different than usual [1], namely such that has the dimension of an inverse energy.

Eq. (3.26) can be obtained by approximating the Lindhard function $\chi_{0}(k,\omega)$ by a “collective” Lindhard function (occasionally also referred to as “one-pole approximation” or “mean spherical approximation”) $\chi_{0}^{\rm coll}(k,\omega)$.

This approximation may be justified in several ways. For further reference, define for any function $f({\bf p},{\bf h})$ depending on a “hole momentum” $|{\bf h}|<k_{\rm F}$ and a “particle momentum” ${\bf p}={\bf h}+{\bf q}$ with $|{\bf p}|>k_{\rm F}$ the Fermi-sea average

[TABLE]

where $V_{\rm F}$ is the volume of the Fermi sphere, $n(k)=\theta(k_{\rm F}-k)$ is the Fermi distribution, and $\bar{n}(k)=1-n(k)$. Replacing the particle-hole energies in the Lindhard function $\chi_{0}(q,\omega)$ by the above “collective” or “Fermi-sea averaged” energies

[TABLE]

leads to a “collective approximation” of the Lindhard function

[TABLE]

and the RPA response function

[TABLE]

The frequency integration (4.2) can then be carried out analytically and leads to equation (3.26).

Another way to justify the collective approximation (4.5) is to approximate the particle-hole band by an effective single pole such that the $m_{0}$ and $m_{1}$ sum rules are satisfied [26, 27]

[TABLE]

This leads to the same form (4.5). A third way to understand the approximation (4.5) is to realize that the Lindhard function can be thought of as the Fermi-sea average of the particle-hole excitations.

[TABLE]

Approximating this average value of the inverse of the particle-hole spectrum by the inverse of the average value,

[TABLE]

also leads to the collective Lindhard function (4.5).

The uniform limit approximation discussed in Section 3.5 amounts to replacing the particle–hole interaction $\tilde{V}_{\!\rm p-h}(k)$ by $\tilde{v}(k)$. We can then use the familiar coupling constant integration (2.26),

[TABLE]

where $\chi_{\lambda}(k,\omega)$ is the RPA density–density response function (4.1) for an interaction $\lambda v(r)$. Using the collective approximation (4.5) in the RPA energy correction, we can carry out both the coupling constant and the frequency integration exactly and end up with the “uniform limit” expression (3.37) [73].

The connection between the chain diagrams of FHNC-EL theory and the RPA ring diagrams can also be made more technical by summing all ring diagrams in CBF perturbation theory. This is a rather tedious procedure which involves calculating all contributions to both the energy and the effective interaction ${\cal W}$ and ${\cal N}$ introduced in Section 3.9 with the momentum flux of a ring diagram in FHNC-EL and infinite order CBF perturbation theory. The plausible result of the calculation is that the contribution of all FHNC ring diagrams which involve the “collective” particle-hole propagator (4.5) are canceled and replaced by the RPA ring diagrams involving the exact Lindhard function. For details, see Refs. 23, 24 and also Ref. 25 for pedagogical material.

4.2 Ladders

We now turn to discuss the connection between the Euler equations of FHNC-EL, specifically Eq. (3.33), and the Bethe-Goldstone equation. A treatment that is similarly rigorous as the one of the ring diagrams is not available for the ladder diagrams that describe short-ranged correlations [85, 86, 87, 88, 89]. We can nevertheless highlight the relationship between the coordinate-space formulation (3.33) of the Euler equation and the Bethe Goldstone equation.

We begin with the Bethe-Goldstone equation as formulated in Eqs. (2.1), (2.2) of Ref. 87. As above, it is understood that ${\bf p},{\bf p}^{\prime}$ are particle states and ${\bf h},{\bf h}^{\prime}$ are hole states; it is then unnecessary to spell out the projection operators.

For the present purpose it is not convenient to go to the center of mass frame. Then, the Bethe Goldstone equation reads

[TABLE]

Following Ref. 87, introduce the pair wave function

[TABLE]

Comparing this with Eq. (4.2) we see that

[TABLE]

This gives the relationship

[TABLE]

for the pair wave function $\psi$. This is obviously still an quantity that depends on three meomenta.

Bethe and Goldstone set the center of mass momentum zero which leaves us still with a function of two variables. This is customary in nuclear physics applications, what follows is what local parquet diagram theory or FHNC-EL would suggest.

Making the connection to FHNC-EL we can proceed again in two different ways. One is to approximate the energy denominator in Eq. (4.15) by its Fermi-sea average, note that the $\bar{n}(k_{i})$ factors in Eq. (4.15) restrict the integration range to unoccupied states.

[TABLE]

Since the bare interaction is, per assumption, a function of momentum transfer, the pair wave function $\psi$ also becomes a function of momentum transfer only. This gives a local equation

[TABLE]

The Fermi-sea integrals of the energy denominator defining $E_{1}(q)$ can be carried out analytically, see problem 1.5 in Ref. 1.

Another way to deal with this, which is more in the spirit of Bethe and Goldstone, is to write Eq. (4.15) as

[TABLE]

Approximating now

[TABLE]

gives

[TABLE]

or, in coordinate space

[TABLE]

Making the connection to the collective response functions of Section 4.1 we note that we can write (4.16) as

[TABLE]

Replacing here $\chi_{0}(q,\omega)$ by $\chi_{0}^{\rm coll}(q,\omega)$ leads to the same expression,

[TABLE]

Thus, the different localization prescriptions discussed in Section 4.1 all lead to the localized Bethe-Goldstone equation (4.19). The direct calculation of the localized energy denominator (4.16) offers an alternative which we have not further pursued.

Both procedures are, of course, approximations; in particular the second form permits to rewrite the Bethe-Goldstone equation in the form of a non-local Schrödinger equation. The legitimacy can be tested by comparing $\left\langle E({\bf q},{\bf k},{\bf k}^{\prime})\right\rangle(q)$ with $E_{1}(q)$ as defined in Eq. (4.16). The maximum deviation of the ratio $\left\langle E({\bf q},{\bf k},{\bf k}^{\prime})\right\rangle(q)/E_{1}(q)$ from 1 is about 18 percent around $k\approx 2k_{\rm F}$, in the relevant regime $0\leq k\leq k_{\rm F}$ the agreement is better than 10 percent.

Comparing Eq. (4.20) with Eq. (3.33), it makes sense to identify $\psi(r)\approx\sqrt{1+\Gamma_{\!\rm dd}(r)}$. Eq. (4.20) is then obtained by the further assumption $\psi^{2}(r)-1\ll 1$. More importantly, the bare interaction of the Bethe Goldstone equation is supplemented by the induced interaction. This has the very important feature that the scattering length of the effective interaction $v(r)+w_{\rm I}(r)$ is zero and that, hence, the pair wave function falls off faster than $1/r$. Historically, a rapid ”healing” of the wave function was often accomplished by introducing a gap in the single particle spectrum [90]. The literature on this subject is vast; the reader is directed to review articles [91, 92] for details.

The identification between the two expressions (4.20) and (3.33) is not as precise as in the case of ring diagrams, but note that FHNC-EL//0 contains more than just particle-particle ladders but also particle-hole and hole-hole ladders [26].

One can also apply the same procedure to define a localized $G$ matrix which is then

[TABLE]

The answer is, for $S_{\rm F}(q)=1$, the well-known Bose limit which should come out. The right-hand side of Eqs. (4.20) and (3.33) manifests the fact that the short–ranged behavior of the pair wave function is modified by the Pauli principle. The effect has already been noted by Gomes, Walecka, and Weisskopf [74].

4.3 Rungs

Let us now turn to the localization procedure of parquet diagram theory. There are two issues to clarify: One is to determine the approximations that are made such that the procedure leads, similar to the Bose case, to the Euler equations for the Jastrow correlations. The second is to generalize the procedure to fermions without such approximations.

We have commented above about the importance of the induced interaction correction $w_{\rm I}(r)$ in the Bethe-Goldstone equation or the coordinate-space Euler equation. From diagrammatic analysis, we should identify $w_{\rm I}(r)$ with a local, energy independent approximation for the energy-dependent particle-hole reducible vertex. Assuming a local particle-hole interaction $\tilde{V}_{\rm p-h}(q)$ we can write the chain approximation for the full, energy dependent vertex as

[TABLE]

and the particle-hole reducible part of that as

[TABLE]

Following Refs. 4, 6, we now define an energy-independent vertex by taking $W(q,\omega)$ at an average frequency $\bar{\omega}(q)$ defined by

[TABLE]

The frequency integral in the second term can be carried out analytically [1]. In the “collective approximation” (4.5) for $\chi_{0}(q,\omega)$, we obtain

[TABLE]

where we recover the $\tilde{W}(q)$ of Eq. (3.49) and the induced interaction (3.30). We can also calculate the averaged frequency

[TABLE]

This is the same as Eq. (8) of Ref. 6 for bosons when we replace $S_{\rm F}(q)\rightarrow 1$.

In fact, the collective approximation is not necessary; the frequency integral can be carried out exactly for the full Lindhard function. Of course, both versions are approximations for the fully energy-dependent induced interaction, their comparison gives information on the robustness of this approximation. We have commented about this already in the section on the ladder diagrams. The above analysis clarifies the relationships between the ring diagrams in parquet theory and those of FHNC-EL.

The above-mentioned experience with the accuracy of the collective approximation for energy calculations must of course be qualified: Approximating the response function by a “collective” response function can, of course, be expected to be a good approximation only if the system indeed has a strong collective mode. This is indeed the case for the density channel in 3He which has, at zero pressure, a Fermi-liquid parameter $F_{0}^{s}=9.15$[93] or for electrons where the plasmon contains all the strength. If $F_{0}^{s}$ is small or even negative, the collective mode is Landau damped. The collective approximation (4.5) in the RPA expression (4.2) is still reasonably accurate, but not at the percent level.

Another issue that needs to be addressed when moving from the Jastrow-Feenberg description to parquet diagrams is the definition of $\tilde{\Gamma}_{\!\rm dd}(k)$. In FHNC//0 we can obtain this quantity from $S(k)$ via Eq. (3.22)

[TABLE]

To construct the equivalent of this relationship in parquet theory, we go back to Eq. (4.25). We should identify

[TABLE]

Here, we encounter again the integral (4.25) which can be carried out either exactly of in the collective approximation. The latter leads to the connection (4.27) as a definition for $\tilde{\Gamma}_{\!\rm dd}(q)$. The FHNC approximation $\tilde{\Gamma}_{\!\rm dd}^{\rm FHNC}(q)$ is then obtained by replacing $\chi_{0}(q,\bar{\omega}(q))$ by $\chi_{0}^{\rm coll}(q,\bar{\omega}(q))$.

4.4 Self-energy

Single-particle properties are in perturbation theory normally discussed in terms of the Dyson-Schwinger equation [94, 95, 96]. While in principle exact, approximations are necessary to make the theory useful. A popular approximation is the so-called G0W approximation [97, 98, 1, 99] for the self–energy,

[TABLE]

where

[TABLE]

is the free single-particle Green’s function. $U(k)$ is a static field, the Fock term in Hartree-Fock approximation, possibly the Brueckner-Hartree-Fock term in Brueckner theory [91].

In the term $\tilde{V}_{\rm p-h}^{2}(q)\chi(q,\omega)$ we recover the RPA for the energy dependent effective interactions $\tilde{w}_{\rm I}(q,\omega)$ (4.24). Setting, as above, $\tilde{w}_{\rm I}(q,\omega)\approx\tilde{w}_{\rm I}(q,\bar{\omega}(q))$ we can carry out the frequency integration and obtain a static self-energy simply in the form of a Fock term in terms of the induced interaction $w_{I}(q)$. Since this sums all chain (or particle-hole reducible) diagrams, the static field $U(k)$ should be the Fock term generated by the particle–hole irreducible interaction, i.e. $\tilde{V}_{\rm p-h}(q)$. To summarize, we obtain in this simplest approximation of the FHNC-EL equations that the numerator $\tilde{X}_{\rm cc}^{\prime}(k)$ of Eq. (3.73) is the same as static approximation for the G0W self-energy when the effective interaction is taken at the average frequency $\bar{\omega}(q)$. Of course, the summation of the full set of the “cc” equations and the associated linear equation for $\tilde{X}_{\rm cc}^{\prime}(k)$ as well as retaining the denominator $1-\tilde{X}_{\rm cc}(k)$ implies that much larger classes of diagrams are routinely calculated in FHNC.

The study of the self–energy highlights, at a somewhat more technical level, another limitation of locally correlated wave functions: Strictly speaking, $\Sigma(k,\omega)$ depends on energy and momentum. Approximating this function by an “average” energy-independent function misses the important non-analytic structure of the self-energy around the Fermi surface. This has the well-known consequence of an enhancement of the effective mass in nuclei around the Fermi surface [100]. The effect is quite dramatic in 3He [101, 102] where a Jastrow-Feenberg wave function predicts an effective mass ratio $m^{*}/m<1$ in massive contrast to the experimental value around $m^{*}/m\approx 3$ [93, 103] which are well reproduced when the full self-energy is calculated [104, 105].

To conclude this discussion, a few more remarks are in order:

The strong enhancement of the effective mass in 3He is due to an interplay between density-fluctuations describing basically hydrodynamic backflow, and spin-fluctuations describing the emission and re-absorption of a low-energy magnon [105]. We have spelled out in Eq. 4.30 only the density channel which is the most important one for out purposes.

• 2.

One can, of course, also use the collective density–density response function (4.6). This would describe the coupling to density fluctuations but treat the coupling to particle-hole excitations only very approximately; in particular it would wipe out the structure of the self-energy around the Fermi momentum. The approximation can be useful for describing the motion of impurities in a Fermi liquid.

5 BCS theory for local correlations

Let us now turn to the generalization of the correlated wave functions method to superfluid systems. Since we have reviewed the FHNC-EL theory and its relation to parquet diagrams above, we can restrict ourselves to the discussion of what changes for a superfluid system. Previous work has either assumed that the superfluid state deviates little from the normal state [36, 30, 31, 28, 29, 35] and/or adopted low-order cluster expansions [106, 107, 108, 109]. The Jastrow-Feenberg variational approach has never been developed to a level comparable to the normal system which made the identification with parquet-diagrams possible. This is one of the tasks of our work.

We construct a correlated wave function for a superfluid system by combining the BCS wave function of a weakly interacting system,

[TABLE]

where $u_{\bf k}$, $v_{\bf k}$ are Bogoliubov amplitudes satisfying $u_{\bf k}^{2}+v_{\bf k}^{2}=1$, with the Jastrow-Feenberg wave function (2.2), (2.3), to the form

[TABLE]

We have commented on alternative choices [110, 107] of the correlated BCS wave function in Ref. 35.

In what follows, we will refer to expectation values with respect to the uncorrelated state (5.1) as $\left\langle\ldots\right\rangle_{0}$ and those with respect to the correlated state (5.2) as $\left\langle\ldots\right\rangle_{c}$. Physically interesting quantities like (zero temperature) Landau potential of the superfluid system

[TABLE]

are then calculated by cluster expansion and resummation techniques; the correlation functions are determined by the variational principle

[TABLE]

5.1 Weakly coupled systems

We have simplified in Refs. 35 and 36 the problem by expanding $\left\langle H^{\prime}\right\rangle_{c}$ (5.3) in the deviation of the Bogoliubov amplitudes $u_{{\bf k}}$, $v_{{\bf k}}$ from their normal state values $u^{(0)}_{{\bf k}}=\bar{n}(k)$, $v^{(0)}_{{\bf k}}=n(k)$. This approach adopts a rather different concept than the original BCS theory: A wave function of the form (5.1) begins by creating Cooper pairs out of the vacuum. Instead, the approach (5.2) begins with the normal, correlated ground state and generates one Cooper pair at a time out of the normal system as suggested recently by Leggett [111]. Adopting an expansion in the number of Cooper pairs, the correlation functions $u_{n}({\bf r}_{1},\ldots{\bf r}_{n})$ can be optimized for the normal system.

Carrying out this expansion in the number of Cooper pairs, we have arrived at the energy expression of the superfluid state

[TABLE]

Above, $H_{\bf o}^{(N)}$ is the energy expectation value of the normal $N$-particle system, and $\mu$ is the chemical potential. The $e_{{\bf k}}$ are the single particle energies derived in correlated basis function (CBF) theory [23], see also Section 4.4. The paring interaction has the form

[TABLE]

The effective interaction ${\cal W}(1,2)$ and the correlation corrections ${\cal N}(1,2)$ are given by the compound-diagrammatic ingredients of the FHNC-EL method for off-diagonal quantities in CBF theory [23], see Section 3.9.

The Bogoliubov amplitudes $u_{\bf k}$, $v_{\bf k}$ are obtained in the standard way by variation of the energy expectation (5.5). This leads to the familiar gap equation

[TABLE]

The conventional (i.e. “uncorrelated” or “mean-field”) BCS gap equation [1] is retrieved by replacing the effective interaction ${\cal P}_{{\bf k}{\bf k}^{\prime}}$ by the pairing matrix of the bare interaction.

5.2 Cluster expansions for a superfluid system

The expansion in the number of Cooper pairs created out of the normal ground state described above is legitimate as long as the superfluid gap function $\Delta_{\bf k}$ is small. When that assumption is not satisfied, one must evaluate all physical quantities of interest for the fully correlated BCS state (5.2). This is the purpose of this section.

The central quantity in the development of our method is the zero-temperature grand (or Landau) potential (5.3) of the superfluid system which we can write as

[TABLE]

For the development of the formal theory, we utilize the methods developed for the cluster expansions of the normal system [23] and outlined in Section 3.9 and modify them for the present case. We can write $\left\langle H\right\rangle_{c}$ in terms of these quantities as

[TABLE]

where $E_{\rm diag}$ is the term in the first line.

The decomposition (5.11) has the following purpose: In what follows, we will see that one can obtain the combination $E_{\rm G}$ by appropriate generalization of the diagrammatic expansions for the normal system, in particular the cluster expansion of this term is irreducible and one can utilize the technique of a “generating function” to obtain it. The “energy numerator” term $E_{\rm enum}$ requires separate treatment; it leads to the second term in the definition (5.6) of the pairing matrix elements.

Carrying out the cluster expansions is a rather tedious matter [112]; their main purpose is to verify, to a convincingly high order, the diagrammatic rules according to which the basic quantities of the theory are constructed. We display here only the characteristic steps and the most relevant result. A few simple examples on how the cluster expansions are carried out will be shown in A and B.

We utilize the “generating function” technique by replacing $u_{2}(r)$ by $u_{2}(r;\beta)$ in the correlated state (5.2). Apart from the “Jackson-Feenberg” kinetic energy terms which originate from the last term in Eq. (2.11) and will be spelled out below, the diagonal and off-diagonal energy terms are

[TABLE]

The purpose of the representation (5.14) is that cluster expansion and resummation techniques for $G$ can be constructed directly from the corresponding expression for the normal system. We have discussed the rules for the normal system in Section 3.2 above. Our diagrammatic analysis for the superfluid system has resulted in the assertion that the corresponding expansion is derived from that by the following rule:

Interpret the density factor as

[TABLE]

where $\Omega$ is the normalization volume. Note that this is the density of the model system described by the uncorrelated BCS state (5.1) and not the density corresponding to the correlate state (5.2). See Eq. (C.20) for the correlation correction to the density. 2. 2.

Re-interpret all exchange lines $\ell(rk_{\rm F})$ as

[TABLE] 3. 3.

In each exchange loop, replace, in turn, each pair of exchange lines $\ell_{v}(r_{ij})\ell_{v}(r_{kl})$ by a pair $-\ell_{u}(r_{ij})\ell_{u}(r_{kl})$, where

[TABLE]

The terms containing only $\ell_{v}(r_{ij})$ lines come from the “diagonal” term $G_{\rm diag}$ and those containing at least one pair of $\ell_{u}(r)$ lines originate from $G_{\rm offd}(\beta)$. The diagrammatic expansion of the generating function $G(\beta)$ can therefore be read off from the diagrammatic expansions for the normal system discussed in Section 4. As an example, we show in Fig. 6 the diagrammatic representation of some leading cluster contributions to $G(\beta)$.

From the cluster expansion for the generating function $G(\beta)$ we can construct the expansion for $E_{\rm G}$

[TABLE]

Note, in particular, that the expansion of $E_{\rm G}$ is irreducible. Diagrammatically, the corresponding expansion of the energy $E_{\rm G}$ can be generated from the generating functional by generalizing the construction rules spelled out in Section 3.2 by

To calculate the kinetic energy corrections $T_{{\rm JF},{\rm G}}$ replace any pair of exchange lines with one common point, $\ell_{\{u,v\}}(r_{ij})\ell_{\{u,v\}}(r_{ik})$, by a pair $\frac{\hbar^{2}}{8m}\nabla_{i}^{2}\ell_{\{u,v\}}(r_{ij})\ell_{\{u,v\}}(r_{ik})$.

The last term, $E_{\rm enum}$ gives rise to the “energy numerator corrections” shown in (5.6). A discussion of the significance of these terms is found in Refs. 35 and 36. Basically, the formulation (5.6)-(5.9) amounts to a reformulation of the gap equation in terms of the $T$-matrix as carried out, for example, in Ref. 113.

Cluster expansions for the energy numerator terms must be treated separately. These corrections consist of a series of products of diagonal matrix elements of the Hamiltonian, and off-diagonal matrix elements of the correlation operator. We classify the individual terms according to the number of common states contained in these two sets of matrix elements. The simplest, and dominating, energy numerator terms are those that have only one common state; in particular these are the only ones that survive in the weakly coupled limit discussed in Section 5.1. The derivation of these terms will be outlined in C. To express the construction rules, we need the following two generalized exchange lines

[TABLE]

where the $e_{{\bf k}}^{(v)}$ are single particle energies generated from the ordinary CBF single particle energies $e_{{\bf k}}$ by replacing all exchange lines by $\ell_{v}(r)$. The leading term in the $e_{{\bf k}}$ is simply the kinetic energy $t(k)$.

The construction rule to obtain these diagrams from the generating functional is then:

Replace, in turn, each $\ell_{u}(r_{ij})$ (red-arrow line) by $\ell_{u}^{\prime}(r_{ij})$ (double red-arrow line) as defined in Eq. (5.19). 2. 2.

Replace, in turn, $\ell_{v}(r_{ij})$ (black-arrow line) by $\ell_{v}^{\prime}(r_{ij})$ (double black-arrow line) as defined in Eq. (5.20). For this, we must interpret the single dot as an $\ell_{v}(r_{ii})$ line returning into itself.

5.3 FHNC and Euler equations

We formulate in this section the simplest version of the FHNC-EL equations, or FHNC//0 approximation. It is straightforward to formulate and implement the full FHNC equations [27]; in fact we have used the FHNC//1 approximation in our numerical applications. As above, the FHNC//0 form displays the physical content of the theory more clearly. Moreover, we shall show that corrections beyond the local correlation function are much more important than adding more complicated FHNC diagrams.

For further reference, we need the following quantities

[TABLE]

as well as their Fourier transforms $\tilde{\sigma}_{v}(k)$, $\tilde{\sigma}_{u}(k)$, $\tilde{\sigma}_{v}^{\prime}(k)$, and $\tilde{\sigma}_{u}^{\prime}(k)$. We begin with the expansion of the “dressed” correlation line $\Gamma_{\!\rm dd}(r)$. In FHNC//0, one keeps only the diagrams shown in Fig. 8 plus, of course, all longer chain diagrams and all parallel connections.

The FHNC//0 equations for a superfluid system are then identical to those for the normal system, except that the structure function of the normal system is replaced by the one of the BCS state,

[TABLE]

where $\hat{\rho}_{{\bf k}}$ is the density operator. It has the form

[TABLE]

From the definitions (5.26), (5.21), and 5.22) it is seen that the long-wavelength limit of $S_{\rm F}(k)$ is

[TABLE]

Hence, $S_{\rm F}(0+)\neq 0$ for the superfluid system. The FHNC equations are, in this approximation, the same as Eqs. (3.20)-(3.22). Note, however, that we can in general not directly construct the pair distribution function $g(r)$ from Eqs. (3.35), (3.36) if energy numerator terms beyond the kinetic energy are included.

The derivation of the Euler equation also is very similar to that for the normal system; the only change is that we must take the energy numerator terms into account. The $S^{\prime}(k)$ then consists of three instead of two terms, the rule spelled out in Section 3.3 has to be augmented by

Replace, in turn, each $S_{\rm F}(k)$ by $\tilde{\sigma}_{u}^{\prime}(k)$. This may be thought of a correction to the second item $S_{\rm F}^{\prime}(k)$.

The remaining derivation is then identical to that for the normal system if we re-interpret

[TABLE]

Inserting this in Eq. (3.25) and the Euler equation (3.24) and solving for $S(k)$:

[TABLE]

where we have defined

[TABLE]

5.4 Energy

Since the $E_{\rm G}$ is, with replacement of the exchange lines by $\ell_{v}(r_{ij})$ and $\ell_{u}(r_{ij})$ described above, identical to that of the normal system, we only need to discuss the correction $E_{\rm enum}$. We include in the energy numerator all those terms that have the same topology as those retained in the above energy expressions (3.34). There are the ones derived in B containing $\tilde{\sigma}_{u}^{\prime}(k)$ and $\tilde{\sigma}_{v}^{\prime}(k)$ as defined in Eqs. (5.24) and (5.23).

[TABLE]

where $\tilde{\Gamma}_{\!\rm dd}^{(v)}(k)$ is represented by the subset of diagrams containing only $\ell_{v}(r_{ij})$ exchange lines. In our numerical calculations we found that the term $\left[\tilde{\Gamma}_{\!\rm dd}(k)-\tilde{\Gamma}_{\!\rm dd}^{(v)}(k)\right]$ is tiny, the term is henceforth omitted.

5.5 Uniform limit approximation

Eq. (5.29) has an interesting consequence: Note that the correction term $R(k)$ does not contain correlation contributions. In other words it as also present even if the interaction is set to zero. The consequence of this is most easily seen in the “uniform limit approximation” [3] discussed in Section 3.5 which amounts to identifying $V_{\rm p-h}(r)=v(r)$ in Eq. (5.29). The energy is, in this approximation [114]

[TABLE]

The first terms in this expression is simply the (Hartree-Fock) energy expectation value with respect to the uncorrelated BCS state

[TABLE]

and the remaining terms have been manipulated, using the Euler equation (5.29) in the uniform limit to eliminate $\tilde{v}(k)$.

Our result (5.32) demonstrates that the energy correction due to correlations is always negative. This result is not entirely surprising, it simply says that adding correlations will lower the energy expectation value. The more interesting statement is that the second term is negative even in a non-interacting system. The reason for this is that the BCS wave function (5.1), when projected to a fixed particle number, can be written in the form of an independent pair wave function [115]

[TABLE]

where ${\cal N}$ is the normalization integral, ${\cal A}$ stands for antisymmetrization, and $\phi(ij)$ is a pair wave function given by the Fourier transform of $v_{\bf k}/u_{\bf k}$. Thus, if we begin with an incorrect assumption about the Bogoliubov amplitudes (they are, of course, equal to normal state values $u^{(0)}_{{\bf k}}=\bar{n}(k)$, and $v^{(0)}_{{\bf k}}=n(k)$ for a non-interacting system) the Jastrow-Feenberg correlations try to compensate for that and lower the energy expectation value.

5.6 Euler equation for the Bogoliubov amplitudes

The study of the Euler equation has so far kept the Bogoliubov amplitudes $u_{{\bf k}},v_{{\bf k}}$ fixed and only dealt with the optimization of pair correlations for such a fixed model state. In particular the analysis of the preceding section showed that the Jastrow correlations and the correlations implicit to the BCS wave function are not completely independent. It is clear that, in an independent step, the Bogoliubov amplitudes must also be optimized.

The derivation of the Euler-equation of the Bogoliubov amplitudes $u_{\bf k}$, $v_{\bf k}$ follows closely the derivation of the conventional gap equation by minimization of the Hartree-Fock approximation [116]. One difference is that one has, in Hartree-Fock-Bogoliubov approximation, only one pair of $\ell_{v}(r)$ and one pair of $\ell_{u}(r)$ lines, see Eq. (5.33). The full variational energy expression $\left\langle H^{\prime}\right\rangle_{c}$ in Eq. (5.11) can, on the other hand, have any number of these exchange lines. This means, in practice, that the effective pairing interaction and the single-particle energies depend implicitly on the $u_{\bf k}$, $v_{\bf k}$. The second difference is the appearance of the “energy numerator” terms $E_{\rm enum}$.

As usual, we guarantee the normalization condition $u_{\bf k}^{2}+v_{\bf k}^{2}=1$ by setting

[TABLE]

The minimization of the energy with respect to the $\chi_{\bf k}$ can be done in both momentum and coordinate space. In coordinate space we have

[TABLE]

or, in momentum space with

[TABLE]

we have

[TABLE]

Then, the optimization of $\left\langle H^{\prime}\right\rangle_{c}$ with respect to $\chi_{\bf k}$ results in

[TABLE]

There are two contributions to the variational derivatives appearing in Eq. (5.39): One comes from $E_{\rm G}$ and the other from the energy numerator terms $E_{\rm enum}$. The variations of $E_{\rm G}$ and $E_{\rm enum}$ with respect to $\ell_{v}(r_{ij})$ or $\ell_{u}(r_{ij})$ are best done diagrammatically by removing, in turn, from the expansion shown in Fig. 7 one exchange line and opening its external points.

The variation of $E_{\rm G}$ with respect to the $\tilde{\ell}_{v}(k)$ gives

[TABLE]

where the $e_{\bf k}$ are the generalization of the single particle energies of CBF theory [23], performed according to the rules formulated in Section 5.2.

• 2.

The variation of $E_{\rm G}$ with respect to the $\tilde{\ell}_{u}(k)$ is constructed according to the same rules. It generally leads to a non-local pairing interaction of the general CBF form [23]. Since we keep, in the FHNC//0 version, only diagrams with $\ell_{v}^{2}(r)$ or $\ell_{u}^{2}(r)$ loops, we can write the result as

[TABLE]

where ${\cal W}(r)$ is diagrammatically obtained by (a) taking all diagrams in $E_{\rm G}$ than contain only $\ell_{v}^{2}(r)$ or $\ell_{u}^{2}(r)$ loops, and (b) removing, in turn, each $\ell_{u}^{2}(r)$ loop and opening its external points.

• 3.

We have restricted the discussion of the energy numerator terms to the simple approximation (5.31), omitting the second term. Then the variation with respect to $\chi_{\bf k}$ yields two terms:

[TABLE]

The second term, which is proportional to $u_{\bf k}^{2}v_{\bf k}^{2}$, contributes a correction $\delta e_{k}$ to the single particle spectrum whereas the first term gives a correction to the pairing matrix element. There is a third term, which we have not spelled out in Eq. (5.42), which originates from the implicit dependence of the single-particle energies $e_{{\bf k}}^{(v)}$ on the $\ell_{v}(r)$. This term also contribute, if kept, to the shift $\delta e_{k}$ of the single particle spectrum.

We can now write the minimization condition for the Bogoliubov amplitudes as

[TABLE]

where

[TABLE]

With that, we have brought the minimization condition for the Bogoliubov amplitudes in exactly the same form as the one for the weakly interacting system [116, 117], the remaining manipulations can therefore be skipped. The only change compared to the weakly coupled limit (5.9) is that all ingredients depend implicitly on the $u_{{\bf k}},v_{{\bf k}}$ and the equation must be solved iteratively. Of course, in a more complete evaluation of the energy numerator terms as outlined in B, the kinetic energy terms in Eq. (5.44) are also replaced by the CBF single particle energies plus corrections from higher order diagrams. Also, note that the single particle energies and the pairing matrix elements depend implicitly on the Bogoliubov amplitudes. We also recover the weak pairing limit (5.6) by setting $u^{(0)}_{{\bf k}}=\theta(k-k_{\rm F})$, $v^{(0)}_{{\bf k}}=\theta(k_{\rm F}-k)$, then $(t(k)-\mu)(u_{{\bf k}}^{2}-v_{{\bf k}}^{2})\rightarrow|t(k)-\mu|$.

5.7 Long wavelength analysis

The result (5.29) points to one of the major problems of local correlation functions and is, as such, one of the key messages of our paper. At the first glance, the result looks innocuous. Recall that we have discussed in Section 3.8 how Landau’s stability condition $F_{0}^{s}>-1$ is a condition for existence of the parquet equations. In particular, within that limit, a negative value of $\tilde{V}_{\rm p-h}(0+)$ is permitted.

If, on the other hand, the system is superfluid, we have $S_{\rm F}(0+)>0$, and therefore $\tilde{V}_{\rm p-h}(0+)+R(0+)$ must be positive to have a solution. $R(0+)$ vanishes in the limit of a normal system whereas $\tilde{V}_{\rm p-h}(0+)$ remains finite. Thus, if  $\tilde{V}_{\rm p-h}(0+)<0$, the Euler equation ceases to have a solution even if the gap is infinitesimally small. We have shown this here only for the case of the FHNC//0 approximation, but it is also true in the general case that all FHNC diagrams, and possibly also higher-order correlation functions, are included since the Euler equation remains structurally the same, see Section 3.7. Our observation applies, of course, equally to “fixed-node” Monte Carlo calculations which may see this instability only in large stochastic fluctuations.

On a less drastic level, assume that $\tilde{V}_{\!\rm p-h}(0+)+\tilde{R}(0+)>0$. Eq. (5.29) then predicts the long-wavelength limit

[TABLE]

which is obviously the wrong behavior since one should have $S(k)\sim\hbar k/2mc$. In other words, even if the interaction is repulsive, $\tilde{V}_{\rm p-h}(0+)>0$, Eq. (5.29) predicts the incorrect behavior of the static structure function at long wave lengths since the contribution of the free kinetic energy is missing.

This is evidently not a statement about the physics, but rather a statement on the approximations implicit to the wave function (5.2), specifically the collective approximation (4.5) for the Lindhard function. One can expect that the correct Lindhard function removes this spurious instability.

There have been several suggestions for a Lindhard function for a superfluid system [118, 18, 19, 21], the most frequently used form for $T=0$ is given below. In view of the need for the spin-spin response function in Section 6.1 we cite here both the spin and the density channels. In terms of the usual relationships of BCS theory,

[TABLE]

with $\xi_{{\bf k}}=t(k)-\mu$ and $E_{{\bf k}}=\sqrt{\xi_{{\bf k}}^{2}+\Delta_{{\bf k}}^{2}}$ we have [115, 119, 120, 118]

[TABLE]

where

[TABLE]

where the upper sign applies to the density channel, and the lower to the spin channel, respectively.

In the limit of a normal system, the coefficient $b_{{\bf k},{\bf q}}^{\rho,\sigma}$ become

[TABLE]

as it should come out.

This superfluid Lindhard function is consistent with the $S_{\rm F}(k)$ as defined in Eq. (5.25).

[TABLE]

We can now return to the frequency integration (4.2). All we need to show is that this expression exists, for small gaps, for $-1<\frac{\tilde{V}_{\rm p-h}(0+)}{mc_{\rm F}^{2}}$. For that, it is sufficient to look at the limit $k\rightarrow 0$ of the static response function and the static Lindhard function.

For $\omega=0$ we get for the static Lindhard function

[TABLE]

This is identical to the same limit of the Lindhard function of the normal system and leads to the correct stability condition.

Thus, the conclusion of our analysis is a more precise version of the statement made in Sections 3.8 and 4 about the validity of locally correlated wave functions for physical phenomena that involve mostly particles close to the Fermi surface. Previous work [28, 29] pointed out quantitative deficiencies; we have demonstrated here a much more profound problem, namely that there are indeed very serious qualitative difficulties in the sense that the minimization problem (2.8) has no solution. We stress again that this feature is not a consequence of the specific level of FHNC approximations. It is a general problem of the Jastrow-Feenberg form of the wave function. More elaborate versions of the FHNC summations and/or the inclusion of triplet correlations can change the numerical values but not the general features.

We conclude this section by recalling that the $\chi_{0}^{(\rho)}({\bf k},\omega)$ does not satisfy the $f$-sumrule. Improved versions have been suggested [18, 21]. It would be very interesting to examine the wave function that corresponds to that work. Recall [121] that, to satisfy the $f$-sumrule, the wave function must be an eigenfunction of the real-world Hamiltonian (2.1), which would be simply the Hamiltonian of a non-interacting Fermi system. Since that Hamiltonian commutes with the particle number operator, the wave function would also be, unlike the state (5.1), an eigenstate of the particle number operator.

6 Applicatios

Most of the calculations reported in this paper use the “weak coupling approximation” spelled out in Section 5.1 with the additional qualification that the “collective” Lindhard function (4.5) that was used throughout our earlier work was replaced by the exact Lindhard function in the rings and rungs. Only very few of the ground state properties changes visibly by including the superfluid Lindhard functions (5.46), we will mention that where appropriate. Only for the calculations of the gap function $\Delta(k)$ we have used also used the superfluid Lindhard function.

We have also used the simple RPA formulas (4.2) and (4.24) to sum the rings and rungs, for that it is appropriate to set $\tilde{R}(k)=0$. A more elaborate version of time-dependent Hartree-Fock is available [20] which should provide an interpretation of the correction $\tilde{R}(k)$; we have not included this improvement because the corrections due to including the functions (5.46) are throughout very small, but computationally very time consuming.

6.1 Neutron Matter

Neutron matter is the first natural application of our method. It is of astrophysical interest because the magnitude of the superfluid gap is of critical importance for the cooling rate of neutron stars [122, 123, 124]. The system has been the subject of two extensive recent reviews [16, 17] which give a very complete account of the current literature. We can therefore restrict ourselves here to those aspects that are specific to the high-level many-body treatment of our paper.

From a quantitative point of view one should question the validity of the “weakly paired” approximation described in Section 5.1 because the gap is, at low densities, of the order of $0.5e_{\rm F}$ [125, 126, 36]. This concern is one of the reasons for developing the methods described in Section 5.

A few additions are necessary to deal with the dependence of the neutron-neutron interaction on the relative spin of the interacting particles. The Jastrow-Feenberg correlation function (2.3) does not contain spin-dependent correlations. It is, of course, straightforward to calculate the energy expectation value for an interaction of the type (6.1). Indeed, the spin-channel of the interaction gives an important correction to the energy; we have included it in Ref. 36. In the Euler equation, the spin-channel operator only contributes to the exchange terms discussed in Section 3.7; this correction is omitted in the FHNC//0 approximation.

Within CBF, it is also straightforward to include the simplest spin-correlations, either in low order perturbation theory [28], or by summing the CBF ring diagrams [24]. The latter provides an RPA-type energy correction from spin-correlations; the procedure then corresponds to an optimized “single operator chain” approximation [127].

6.1.1 Energetics

Let us now go through the relevant computations step-by-step. All intermediate results that are shown are for the Reid $V_{4}$ soft-core potential which we write in the operator basis [11]

[TABLE]

We have carried out the identical calculations for the Argonne $V_{4}^{\prime}$ nucleon-nucleon interaction [128]. The results are sufficiently similar to those of the Reid potential; we will therefore display only the final results for that case.

The essential difference compared with what we have discussed in previous work [35, 129] is the operator form (6.1) of the interaction. Instead of a local ${\cal W}(r)$, we have the spin–dependent interaction

[TABLE]

where

[TABLE]

in $V_{\rm ee}(r)$. Since the wave function (2.3) contains no spin correlations, there are no chain diagrams contributions to $\tilde{W}(k)$ in the spin channel. In the CBF or parquet calculation, the appropriate interaction is then

[TABLE]

where

[TABLE]

We have commented above about the fact that the two ways to calculate $F_{0}^{s}$, Eqs. (3.61) and (2.33) agree only in an exact theory. To assess the consistency between the two definitions, we have fitted, in the regime $0\leq k_{\rm F}\leq 0.5\,{\rm fm}^{-1}$ the equation of state by a function

[TABLE]

and calculated $F_{0}^{s}$ from both the fit (2.33) and the long-wavelength limit (3.61). Here, and throughout the rest of this paper, we will set the effective mass ratio $m^{*}/m=1$. Justification for this is derived from calculations of the self-energy as described in Sections 6.1.3, 6.2.3 and 6.3.1. Fig. 9 shows the situation for the low-density equation of state of neutron matter, interacting via the Reid $V_{4}$ potential.

We find in particular the FHNC//0 approximation for the particle–hole interaction leads to an incorrect positive value of $F_{0}^{s}$. The FHNC-EL//1 approximation improves the agreements by including the exchange term in the particle-hole interaction, but it still disregards the chain-diagrams in the spin-channel who seem to be insignificant. These are included in the parquet calculations. The agreement between $F_{0}^{s}$ calculated in these two ways is, for low densities, evidently quite satisfactory.

6.1.2 Effective interactions

The energetics of neutron matter is relatively insensitive to the quality of the many-body wave function. We have already seen that this is not the case for the quasiparticle interaction. Similarly, one should expect visible corrections to the pairing interaction where medium polarization effects, especially due to spin fluctuations, are expected to have an impact on the superfluid gap [130, 131, 132, 133].

A very thorough examination of polarization effects has been performed by Schulze et al. [134]; our calculations go beyond that work in the sense that we determine the effective interactions by summing the parquet-diagrams, and take the superfluid particle-hole propagators. We have included the exchange terms that have been included in Ref. 134 in the localized form described in Section 3.7. Following the work of Ref. 134, we have also taken the dynamic interactions $\tilde{W}^{(\rho,\sigma)}(q,\omega)$ and $\tilde{W}^{(S)}(q,\omega)$ at $\omega=0$. The interactions derived from the FHNC approximation directly are, of course, energy independent.

A calculation for the superfluid phase that goes beyond the “weakly coupled” approximation requires the implementation of the full parquet-level theory because the Euler equations have, due to $F_{0}^{s}<0$, no solution. Among others, we need the rather time-consuming computation of the superfluid Lindhard function which is an essential input for the theory. The fact that this is indeed necessary is documented in Fig. 10. There we show, for $\omega=0$, the Lindhard function of the normal system, as well as the density- and the spin-channel of the Lindhard function of the superfluid system, Eqs. (5.46).

The discrepancy between the density and the spin-channel superfluid Lindhard functions is quite interesting and has, to our knowledge, not been noted in neutron matter calculations. It is caused by the fact that the superfluid Lindhard function in the spin-channel, Eq. (5.46) goes to zero as $q\rightarrow 0$ for all values of the gap. Of course, the regime where the spin-channel function deviates from the density-channel function becomes smaller with decreasing gap. The density-channel function is expectedly close to the one of the normal system even when the gap is relatively large.

This result has significant consequences for the effective interactions (6.5). These interactions are shown in figures 11 and 12

The conclusion to be drawn from Fig. 11 is that the use of the dynamic effective interaction in the density channel potential at $\omega=0$ makes the interaction somewhat less attractive, but there is practically no difference between using the normal system Lindhard function or the one for the superfluid one. The situation is quite different in the spin-channel: Using the superfluid Lindhard function changes the effective interaction in that channel significantly, note that the gap at the density considered here is only about 20 percent of the Fermi energy, see Fig. 15 below. This is to a large extent due to the fact that superfluid $\chi_{0}^{(\sigma)}(q,0)$ goes to zero at zero momentum transfer. The particle-hole interaction and the paring interaction are therefore identical in that limit. This is not the case if one uses the normal system Lindhard function to calculate $W^{(\sigma)}(q)$.

In the pairing matrix elements, we must couple the interactions (6.5) into the singlet and triplet channel:

[TABLE]

The most relevant quantity for our purposes is the singlet pairing interaction $W^{(S)}(q,\omega=0)$ which is shown in Fig. 12. Similar to the spin-channel, using the superfluid Lindhard function changes the interaction visibly, in particular for long wavelengths $0\leq q\leq 2k_{\rm F}$ which is the important regime for pairing phenomena. Interestingly, the choice of the superfluid Lindhard function suppresses the correction from spin-fluctuations; the $S$-wave interaction is closer to the one where the spin-channel is omitted altogether than to the one where the spin part is taken from the normal system.

6.1.3 Self energy

We have spelled out in Eqs. (3.73), (3.74) and (3.75) working formulas for the single particle spectrum based on a local correlation operator. In Section 4.4 we have then made the connection to the G0W approximation and have shown how these two procedures are related. The situation is similar to the one for other quantities: FHNC offers a much more extensive summation of diagrams. Note, for example, that the denominator $1-\tilde{X}_{\rm cc}(r)$ can be related to high-order exchange diagrams but has no equivalence in the G(0)W approximation. The price for that is that the evaluation of these diagrams is less accurate. We shall address these issues now.

We are examining the self-energy here for a number of reasons. One is to justify our choice of keeping only the kinetic energy in $\ell_{v}^{\prime}(r)$ and $\ell_{i}^{\prime}(r)$, Eqs. (5.19) and (5.20), defining $E_{\rm enum}$. The second reason is to justify the use of a free spectrum in the gap equation (5.9). We have studiend the sensitivity of the superfluid gap agaist changes of the effecive mass in Ref. 36.

The self energy (4.30) is most conveniently evaluated by Wick rotation in the complex $\omega$- plane. That way, the RPA sum for the self-energy is decomposed into two terms, a smooth “background–” or “line–term”, which consists of the frequency integral along the imaginary $\omega$ axis, and a second, “pole term” from the residue of the single-particle Green’s function, i.e.

[TABLE]

with

[TABLE]

and

[TABLE]

Since the effective mass can have an visible effect on pairing properties, we have calculated the ratio

[TABLE]

from the on-shell self-energy $\Sigma(k,t(k))$ for both the Reid and the Argonne potential; results are shown in Fig. 13. The effective masses are all close to 1, suggesting that the single particle spectrum of non-interacting fermions is an adequate approximation in this particular system. It also justifies the on-shell approximation. While the effective mass ratio per-se is close to 1, we see that the corrections from going beyond local correlations are substantial on that scale.

It is also of interest to look at the individual contributions to the self-energy as a function of wave number. We show these in Fig. 14.

The FHNC expression for the static field (3.74) and the “Fock” differ only by the interaction entering the expression: Whereas tilde $X^{\prime}_{\rm cc}(k)$ is defined in terms of the full static interaction ${\cal W}(r)$, the static field $U(k)$ contains only the particle-hole irreducible interaction $V_{\!\rm p-h}(r)$. Evidently these two terms are numerically rather similar. The remaining two terms, $\Sigma_{\rm line}(k,t(k))$ and $\Re e\Sigma_{\rm pole}(k,t(k))$, sum to practically a constant, hence the real part of the spectrum is, at this density, indeed rather similar to the FHNC approximation. We also see that there is, in neutron matter, no visible enhancement of the effective mass around the Fermi surface as observed in nuclear matter [100]. This is consistent with the fact that imaginary part of the self-energy is rather small.

These results essentially confirm that our choice of $m^{*}/m=1$ is justified.

6.1.4 BCS pairing

We have generally solved the full gap equation (5.9), using the eigenvalue solver with an adaptive mesh described in Ref. 35. Results will be reported for the value of the gap function at the Fermi surface, $\Delta(k_{\rm F})$.

Several approximate formulas for the superfluid gap at the Fermi surface are available: At very low densities, the gap can be expressed in terms of the vacuum scattering length $a_{0}$ [113] as

[TABLE]

The exponential $\exp\left(\frac{\pi}{2a_{0}k_{\rm F}}\right)$ is universal in the sense that it depends only on the product $a_{0}k_{\rm F}$, but not on details of the interaction. The pre-factor is modified by polarization corrections [135] which are also universal, and by many-body and finite-range corrections some of which are non-universal [35].

Eq. (6.13) gives, apart from very low densities, quite poor predictions. More useful is the estimate

[TABLE]

where

[TABLE]

is the $S$-wave matrix element of the pairing interaction at the Fermi surface. Eq. (6.14) shows clearly the influence of the effective mass correction; we have studied this in Ref. 36. Since the effective mass is always close to the bare mass, see Fig. 13, we have throughout this paper used $m^{*}/m=1$.

Our results for the superfluid gap are shown in Figs. 15, they basically reflect our findings on the pairing interaction. Including exchanges does not change the gap significantly; the most significant change comes from using the dynamic interaction $\tilde{W}^{(S)}(q,\omega=0)$, Eq. (6.7). Using the superfluid Lindhard function moves the result a little up but is somewhat unexpectedly not very significant, despite the suppression of the long-wavelength components in $\chi_{s}^{(\sigma)}(q,\omega=0)$ as demonstrated in Fig. 10. The approximation (6.14) underestimates $\Delta(k_{\rm F})$ by about a factor of 2. We also show the results obtained from the estimate (6.13) using the experimental $S$-wave scattering length a_{0}\approx-18.7\fm [136].

We have, in the neutron matter calculations, only solved the parquet- equation in the density channel because our primary interest is the pairing interaction; in the spin–channel we have only kept the polarization diagrams. One can extend the calculation to include spin- and tensor- operators [137] in the full parquet summation.

6.2 The Lennard-Jones liquid

Fermi fluids interacting via a family of Lennard-Jones (LJ) interactions

[TABLE]

have been studied extensively in Refs. 72 and 35. The system is interesting from the viewpoint of many-body theory since the equation of state can have two spinodal points where the hydrodynamic speed of sound (2.33) vanishes. For the sake of discussion, we show in Fig. 16 a schematic equation of state for a typical self-bound Fermi liquid.

The dominant energy contribution at very low densities is the Pauli repulsion, $E/N\propto k_{\rm F}^{2}$. As the density increases, interaction terms, which grow as $k_{\rm F}^{3}$, begin to dominate and bend, if the interaction is sufficiently attractive, the equation of state downwards. This leads to a local maximum. If that is the case, the equation of state also has a local minimum at some finite density before the short-ranged repulsion begins to dominate. Hence, there must be two spinodal points. The Lennard-Jones liquid shows, for sufficiently large values of the interaction strength $\epsilon$, both of these points. In contrast to that, an attractive square-well potential or the attractive Pöschl-Teller interaction studied in Ref. 125 and in Section 6.3, has only the lower spinodal point. These systems collapse, due to the absence of a short-ranged repulsion, when the density is increased to a level where the interactions begin to dominate over the Pauli pressure. Neutron matter is not self-bound and does not show spinodal points.

6.2.1 Energetics and stability

We discuss in this section our results for the above family of Lennard-Jones 6-12 interactions. The strength $\epsilon$ or, better, the dimensionless strength parameter $V_{0}=2m\epsilon\sigma^{2}/\hbar^{2}$ of the interaction can be adjusted such that the interaction has the desired vacuum scattering length; a bound state and a corresponding singularity of the vacuum scattering length appears at $V_{0}=11.18$. Fig. 17 shows the scattering length for both the Lennard-Jones and the Pöschl-Teller (PT) interaction discussed in the next section as a function of the potential well depth. We also show, for comparison, the scattering length for the attractive square-well interaction discussed in Ref. 35.

Before we discuss our findings, we comment on the expected accuracy of our results. This was examined in Ref. 72 where we found that the FHNC//0 approximation is accurate within a percent at low densities up to approximately 25 percent of the saturation density. This corresponds to a Fermi wave number of $k_{\rm F}\sigma\leq 0.7$ or a density $\rho=0.0116\sigma^{-3}$. This is the regime of interest here. Around saturation, the FHNC//0 approximation still displaus the correct physics, but leads to an equation of state that is somewhat too soft. The numerical values in that regime should be therefore be considered only qualitative. These observations apply, of course, predominantly to the energy which is relatively insensitive to the quality of the wave function. The stability range is indeed affected, even at low densities, by including the exchange diagrams discussed in Section 3.7.

We have studied here a density regime of $0.01\,\leq k_{\rm F}\sigma\,\leq 2\sigma$, and interactions strengths $1.0\leq V_{0}\leq 9.0$. In Ref. 35 we went to the much smaller values of $0.001\leq k_{\rm F}\sigma$. A proper treatment of many-body effects requires that the simulation box is much larger than $(k_{\rm F}\sigma)^{-1}$. On the other hand, the core region must be properly resolved, hence decreasing the minimum value of $k_{\rm F}$ by an order of magnitude increases the necessary number of mesh points on an equidistant mesh (as required by the fast Fourier transformation algorithm) by an order of magnitude. In Ref. 35 we have worked with $2^{18}$ mesh points; since the implementation of Eq. (4.2) required by the parquet-level calculations implies an energy integration for each mesh point, we have reduced in this work the number of mesh points to $2^{15}$. This limits the minimum density to the above value. Since the lower spinodal points move, with interaction strength, to lower densities, this constraint also limits us to interaction strengths $V_{0}<9.0$.

The equation of state of for a sequence of interaction strengths $1\leq V_{0}\leq 8.01$ is shown in Fig. 18. Results for stronger interactions $8.01<V_{0}<\leq 9$ are not visible in the figure and have been omitted. The figure shows clearly the gap in the equation of state between the two spinodal points which appears when the interaction is sufficiently strong, $V_{0}>6.0$.

The most interesting result of Ref. 35 is that the FHNC-EL equations diverge, at low density and as a function of density, well before the lower spinodal point is reached. To examine this effect more closely, we have included the two improvements discussed in Sections 3.7 and 4.1. As mentioned above, the inclusion of the exchange diagrams shown in Fig. 5 can indeed change the long-wavelength limit of $V_{\rm p-h}(0+)$ by a factor of $1-1/\nu=1/2$ and, hence, can in principle increase the stability regime significantly.

The second correction is due to the transition from the variational wave function to parquet diagrams. A summary of our results for the Fermi liquid parameter $F_{0}^{s}$ is, for the sequence of interaction strengths used in Fig. 18, shown in Figs. 19.

The most visible (and expected) quantitative correction is caused, at low densities, by the inclusion of exchange diagrams. Fig. 20 shows, for four selected densities close to and above the density where spinodal points exist, the $F_{0}^{s}$ as calculated on the FHNC//0 approximation and the parquet//1 approximation. A number of observations apply: First we notice that the results with and without including exchange diagrams differ, in particular in the low density regime, by roughly a factor of 2. This is consistent with the argument of Section 3.7 that the leading term in the density expansion is reduced by exchange diagrams. The factor of $1-1/\nu=1/2$ is visible in the slope of $F_{0}^{s}$ at $k_{\rm F}\rightarrow 0$ as a function of $k_{\rm F}$. (The argument does not apply to the nuclear case since the interaction in the exchange channel is different from the one in the direct channel, see Fig. 9 and Eq. (6.3).) The second, and much more interesting, observation is that the density where the FHNC-EL equations diverge is almost independent of the inclusion of exchange diagrams.

The FHNC and the parquet results are, at the same level of implementation, almost indistinguishable. The inclusion of the correct particle-hole propagator only limits the regime of existence of solutions to the regime $F_{0}^{s}>-1$ but has otherwise little effect. Evidently, the existence of the spinodal points does not change qualitatively with the parquet theory developed in Section 4, however, the location of the divergence changes somewhat.

To assess the accuracy of our predictions we have again carried out the consistency test between the between $F_{0}^{s}$ obtained from the equation of state, (2.33) and from the long-wavelength limit (3.61), $\tilde{V}_{\rm p-h}(0+)$. At the strongest couplings, $2m\epsilon\sigma^{2}/\hbar^{2}>7.75$, the coefficient $b$ of the $k_{\rm F}^{5}$ term cannot be determined reliably, we have therefore kept only the linear in Eq. (6.6). A comparison of the two procedures is shown in Fig. 21; evidently the agreement is in the regime of interest not quite as good as in the nuclear case, this is because the density dependence of the correlations cannot be ignored.

As in previous work [35] we found that it is impossible to find solutions of the FHNC-EL or parquet equations close to the lower spinodal point. The fact that the equations of state all come to an endpoint has been identified in Ref. 35 as due to a divergence of the in-medium scattering length which is, in the local approximation used here

[TABLE]

This divergence is the reason that the Landau stability limit $F_{0}^{s}\rightarrow-1$ could not be reached. The same situation occurs, expectedly, in the present case, see Fig. 22. One might have argued that the divergence found previously is indeed just the spinodal instability and one could not get closer due to some numerical problem. The comparison of the FHNC//0 and the parquet//1 results shown in Fig. 20 shows that the divergence found here is indeed unrelated to the Fermi liquid parameter: It appears at practically the same density, although the $F_{0}^{s}$ differ by a factor of 2.

Thus, our conclusions on this topic have not changed, moreover, we have eliminated the most relevant concern about our previous work, namely that the Landau parameter $F_{0}^{s}$ has not been calculated accurately enough.

6.2.2 Effective interactions and correlations:

A configuration space view

We have throughout the preceding sections emphasized the momentum space properties on effective interactions, in particular the connection between the long wavelength limits, Fermi liquid parameters, the in-medium scattering length, and the pairing interaction. Of course, the structure of correlations is mostly determined by the short–ranged properties of the interactions. Given the knowledge about the bare interactions, the short–ranged structure of the correlation functions is given by the Bethe–Goldstone equation which is, in its essence, a zero-energy Schrödinger equation where the short-ranged structure of the pair wave function is modified by the Pauli principle [88, 87, 74].

The relationship between the pair correlation function $f(r)$, the “dressed” correlation function $\Gamma_{\rm dd}(r)$ and the pair distribution function has been discussed in many places, we therefore focus here on the low-density expansions (3.35), (3.36). Fig. 23 shows, at the relatively low density of $k_{\rm F}=0.3\sigma^{-1}$ a sequence of pair correlation and distribution functions. The “dressed” correlation function $1+\Gamma_{\!\rm dd}(r)$ describes the dynamic correlations and is determined by the bare interaction. The pair distribution function contains, in addition, statistical correction manifested predominantly in the factor $g_{F}(r)$ in Eq. (3.35), but also in the corrections $C(r)$ and $\tilde{X}_{\rm ee}(r)$, cf. Eq. (3.36). Fig. 23 also shows the simple approximation $\left[1+\Gamma_{\!\rm dd}(r)\right]g_{\rm F}(r)$ for the pair distribution function. Evidently, the correction from chain diagrams, $(S_{\rm F}^{2}(k)-1)\tilde{\Gamma}_{\!\rm dd}(k)$ is not negligible; we assert, however, that the contribution from the exchange parts $(\Delta\tilde{X}_{\rm ee})_{1}(r)$ are invisible on the scale of the plot.

The dynamic correlations determine the effective interactions $V_{\!\rm p-h}(r)$ via Eqs. (3.27) and (3.28) and $W(r)$. Aldrich and Pines [83, 84] give a physically intuitive description of the effects contributing to $V_{\!\rm p-h}(r)$ which is called “pseudopotential” in that work:

At short distances, the interaction is screened by short-ranged correlations; 2. 2.

The fact that the wave function is bent at short distances – downwards for repulsive interactions, upwards for the attractive interactions as in the attractive square well potential and the Pöschl-Teller potential – has a price in kinetic energy; 3. 3.

Finally, the presence of other particles should lead to an enhancement of the interaction in the attractive regime.

Eqs. (3.27) and (3.28) give a quantitative meaning to these effects: The short-ranged screening (1) is described by the factor $1+\Gamma_{\!\rm dd}(r)$, the cost in kinetic energy is described by the term $\frac{\hbar^{2}}{m}\left|\nabla\sqrt{1+\Gamma_{\rm dd}(r)}\right|^{2}$, and many-body effects are described by the last term $\Gamma_{\!\rm dd}(r)w_{\rm I}(r)$. Fig. 24 shows these three parts of the particle-hole interaction at $k_{\rm F}\sigma=0.3$ and an interaction strength $V_{0}=7.0$. We had to scale the induced interaction term by a factor of 100 to make it visible at the scale of the plot. This does not mean that this term is negligible; it guarantees that the dynamic correlation function $\Gamma_{\!\rm dd}(r)$ falls off as $r^{-2}$ as $r\rightarrow\infty$. cf. Eq. (3.39).

Fig. 25 shows, finally, the dependence of $V_{\!\rm p-h}(r)$ at $k_{\rm F}\sigma=0.3$ on the interaction strength for the sequence $V_{0}=1,2,\ldots 7$. We have not included the exchange correction, it satisfies to an extremely good accuracy our estimate $V_{\rm ex}(r)=(g_{\rm F}(r)-1)W(r)$.

6.2.3 Effective mass

The Lennard-Jones $6-12$ interaction may be thought of a model for the interaction between helium atoms. More modern interactions [138, 139, 140] are perhaps more accurate; but the Lennard-Jones model reflects the correct physics of 3He. In our parametrization that uses the core size $\sigma$ as unit of length, and $\hbar^{2}/2m\sigma^{2}$ as unit for the energy, 3He is characterized by $V_{0}=8.26$ and the equilibrium density has a Fermi momentum $k_{\rm F}\sigma=2.006$, in other words the equilibrium properties of 3He are beyond the regime where out approximations are reliable. See Ref. 72 for a discussion and full FHNC-EL calculations.

Liquid 3He is known to undergo below 3 mK superfluid phase transitions [141, 142] to two different phases of $P$-wave superfluidity. At the experimental saturation density of 3He, the soft spin-fluctuation mode is very important and is the essential mechanism for the strong effective mass enhancement around the Fermi momentum [93, 103, 104, 105]. If both backflow and spin-fluctuations are omitted, one ends up with an effectiv mass less than 1 [143]. We are in this work interested in low-density systems and potential $S$-wave pairing. We will see below that $S$-wave superfluidity can occur only at very low density; in particular below the lower spinodal point which is around $k_{\rm F}\sigma\approx 0.13$ for $V_{0}=8.26$, see Fig. 19. Therefore, we expect that the much simpler G0W approximation spelled out in Eqs. (4.30), (4.31) is adequate.

Fig. 26 shows the effective mass obtained as described above way for a sequence of interaction strengths $V_{0}$. Throughout, the effective mass never differs from the bare mass by more than 10 percent. Interestingly, we find for the most repulsive interactions a slight reduction of the effective mass ratio which turns into an enhancement around $V_{0}=4.0$ and then above $V_{0}=7.25$ again into a reduction.

The explanation for that is that the value of the effective mass ratio is an effect of both Fermi statistics and hydrodynamic backflow. Fermi statistics has mostly (but not always) the effect of reducing the effective mass, this is a well-known effect in nuclei [100] but also comes out for 3He in the naïve formulation (3.73) [127]. Hydrodynamic backflow [144] increases the effective mass. That effect is included in the G0W approximation which includes hydrodynamic backflow; see Ref. 27 how Feynman’s backflow operator comes out already in second order CBF perturbation theory. In this work, we have not included spin fluctuations, hence the results shown are only due to Fermi statistics and backflow. Apparently, backflow effects become stronger with increasing attractiveness of the interaction. This makes sense because more particles are relatively close to each other; see Fig. 23 how the nearest neighbor peak increases. For that to be effective, of course the density must be high enough.

6.2.4 BCS pairing

Let us now turn to the calculation of the superfluid gap. By construction, our pairing interaction ${\cal W}(1,2)$ contains polarization corrections, see the discussion of the connection between parquet theory and FHNC-EL in Section 4, cf. Eqs. (4.23). These polarization corrections are important modifications of the pairing interaction and can lead to significant changes of the superfluid gap, see, for example, Refs. 131, 132, 133, 134

The details of the choice of interactions (6.5), i.e. if exchange diagrams are included and whether we take the collective Lindhard function or the energy dependent one, define different implementations of the pairing interaction in FHNC or parquet diagram theory. Using local correlations means that the Lindhard function $\chi_{0}(q,\omega)$ is replaced by the collective Lindhard function $\chi_{0}^{\rm coll}(q,\bar{\omega}(q))$ defined in Eq. (4.5), and taken at the average frequency $\bar{\omega}(q)$, Eq. (4.27). In FHNC//0 approximation, this leads to the expression (4.26) which was used in previous work [35, 36]. The equivalence between FHNC and parquet-diagrams permits us, however, to be more flexible. Since the pairing gap is expected to be relatively small one is led, in the parquet version the theory, to take the pairing interaction at zero frequency [134].

Our results are shown in Fig. 27. Indeed, both modifications, exchange diagrams and using the sero-energy Lindhard function for the particle-hole propagator can cause a change of up to a factor of two in the superfluid gap; the effect is about half due to exchange diagrams, and half due the propagator correction. The right pane of Fig. 27 shows a comparison of the solution of the full gap equation with the two approximations (6.13) and (6.14).

A somewhat surprising result is that fully self-consistent theory where all particle-hole propagators (5.46) are taken for the superfluid system causes, despite the smallness of the gap, a visible effect in the stability regime as indicated by the end of the respective lines. Apparenty, being close to an instability, makes every quantity sensitive to small details.

In the model we have considered here it seems impossible to get close to the “BCS-BEC” crossover regime. The reason for lies both in physics and in formalism: From Fig. 19 it is clear that one is, whenever the superfluid gap is sizable, close to the lower spinodal point, even if one cannot reach it. For systems like the one studied here, and, even more, for the attractive Pöschl-Teller interaction studied in the next section, one has two competing effects: To get a stronger gap, one needs a more attractive interaction. Such an attractive interaction has, on the other hand, the effect that the Pauli-pressure is more easily overwhelmed by the interaction and, hence, the spinodal point moved to lower densities.

The formal reason is also quite plausible: The implementation of FHNC-EL or parquet theory involves two sets of diagrams, rings and ladders. A third set, the “cyclic chain” diagrams that would indicate a divergence of the single particle propagator which might be present in spin-polatrized 3He [145] has not been considered here. A divergence of the ring-diagrams is related to the existence of spinodal points or, more generally, to Landau’s stability conditions $F_{\ell}^{s,a}>-(2\ell+1)$. What we find, however, is a divergence in ladder equation (3.33). It is harder to associate that with a specific quantity because the induced interaction $w_{\rm I}$ always tries to adjust itself such that the $v(r)+w_{\rm I}(r)$ has zero scattering length and that, therefore, Eq. (3.33) has a solution with the property $\Gamma_{\!\rm dd}(r)\propto r^{-2}$ as $r\rightarrow\infty$, see Eq. (3.39). But it is easy to envision that the potential becomes attractive enough such that the equation develops a bound state. This has the consequence that $1+\Gamma_{\!\rm dd}(r)$ falls off exponentially as $r\rightarrow\infty$. Then the basic assumptions of the equations of the theory are violated and it is no surprise that they have no solutions. This is apparently what happens here, where the effect is caused by the attractiveness of the induced interaction which describes predomilantly phonon exchange. Hence, we have concluded in Ref. [35] that the divergce is due to a dimerization caused by phonon exchange, in other words it is a genuine many-body effect.

This is an interesting statement per-se: Note that in a weakly interacting system, the gap equation can describe the transition between a “BCS” state where the Cooper pairs are weakly coupled, to a “BEC” phase where the pairs are strongly bound [146]. The gap equation is a proper subset of our diagram summation; but we find a divergence of a set of diagrams that are not included in weakly interacting systems described by the gap equation alone. The issue deserves further investigation.

6.3 Pöschl-Teller interaction

The purely attractive Pöschl-Teller (PT) interaction [147]

[TABLE]

is characterized by the strength $V_{0}$ and the range $\sigma$. For this interaction, the scattering length can be calculated analytically [148].

[TABLE]

where $\gamma$ is Euler‘s constant. The interaction is somewhat more convenient to use in Monte Carlo calculations than the square-well interaction [149]; it has been used in Ref. 125 to study pairing in the unitary limit. We have recently examined the ground state and pairing properties of the Pöschl-Teller gas [129]; in this work we have improved these calculations by applying the improvements outlined above. Since many techical details have been discussed in the preceding sections, we can restrict ourselves to display just the new aspects of our calculations.

6.3.1 Energetics and stability

Since the Pöschl-Teller interaction does not have a repulsive core, there is no self-bound, high-density phase. The low density properties are basically determined by the balance between the Pauli repulsion and the attractive interparticle attraction. Our findings for the PT interaction are sufficiently similar to those of the LJ interaction sich that we can be very brief..

Fig. 28 shows the equation of state, normalized to the kinetic energy of the free Fermi gas. In the limit of low densities, the equation of state should be give by the Huang-Yang expansion (3.60), we have studied this in Ref. 35 for the Lennard-Jones and the attractive square-well model and found that the term proportional to $(a_{0}k_{\rm F})^{2}$ is already overshadowed by non-universal many-body effects. We compare therefore only the linear term in the expansion (3.60) and restrict, of course, the expansion to small values of $a_{0}k_{\rm F}$. Similar to our previous work we find that the expansion (3.60) is reasonably accurate for $a_{0}k_{\rm F}<0.05$, see Fig. 29.

Fig. 30 shows the Fermi liquid parameter $F_{0}^{s}$ as obtained from Eq. (3.61) for the same sequence of interaction strengths as used in Fig. 28 as a function of the density. In general, there is little qualitative change; the stability regime is changed visibly by including exchanges, but the result that there is no way to get close to Landau instability regime $F_{0}^{s}\rightarrow-1$ is not changed; in fact, the most advanced calculation leads, similar to the Lennard-Jones model, to a much smaller value of $F_{0}^{s}$ where the divergence occurs. Due to this instability we have not been able to reach the rather large values of $-k_{\rm F}a_{0}$ reported in Ref. 125 before the system became unstable. We have again checked the consistency between Eqs. (2.33) and (3.61) with the procedure outlined in the previous section. We found that the numerical values are practically identical for weak couplings when exchange diagrams are included.

For completeness, we show in Fig. 31 the effective mass ratio for the above sequence of coupling strengths. As in our previous calculations, these ratios are close to 1 which is certainly a consequence of the low density of the system.

6.3.2 BCS pairing

Fig. 32 shows the calculated energy gap in FHNC//0, parquet//1 and fully self–consistent approximation. Evidently, inclusion of the energy-dependent effective interaction can change the value of the gap visibly. This is, of course, not a statement on the specific FHNC approximation, but more generally on the quality of the locally correlated (or “fixed-node”) wave function. The fully self-consistent calculation brings us here back into the vicinity of FHNC//0, we hasten to point out that this observation is circumstantial and should not be generalized.

In conclusion, we note again that the most important change duu to going beyond the “weak coupling” approximation (5.5) is that the stability regime is further reduced.

7 Discussion

We have formulated in this paper the variational and parquet approach for strongly interacting, both normal and superfluid, systems. We have then studied the Euler equation for the local correlations in a version that exhibits all the exact properties of variational and parquet-diagram theories but does not require the formulation of the full FHNC theory. We have derived the interesting, yet disturbing, result that the Euler equation for locally correlated wave functions of the form (2.3) has, for a superfluid system, no physically meaningful solution for net attractive interactions $\tilde{V}_{\rm p-h}(0+)<0$, and displays pathological features even for repulsive interactions.

The problem is caused by the very form of the Jastrow-Feenberg wave function. The plausible way solve this problem is by comparison with parquet-diagrams. This has been carried out in Section 4, and that formulation can be naturally carried over to the superfluid system. That has the effect that the “collective approximation” for $S(k)$, Eq. (5.29) is replaced by the RPA expression (4.2). The formulation of Coupled Cluster theory with correlated wave functions [31] would perhaps be another way to carry this through, but, on the other hand, our result is sufficiently convincing.

We have then carried out numerical calculations for different model systems of physical interest: neutron matter and model systems interacting via a family of Lennard-Jones and Pöschl-Teller interactions. With the recent interest in pairing in cold gases and the BEC-BCS crossover, it has become fashionable to stress the similarity between cold gases and neutron matter. In fact, the neutron scattering length -18.7\pm 0.6\fm [136] which is much larger than the range of the interaction. We found in this work that these systems are actually rather different. Whereas there was no problem calculating the properties of neutron matter at all densities, the parquet-equations showed inevitable divergences for both the Pöschl-Teller and the Lennard-Jones interactions at scattering lengths of the order of $a_{0}=5\dots 8\sigma$. The divergence is driven by the approaching of a “soft mode” at the spinodal points(s) $F_{0}^{2}\rightarrow\rightarrow-1$, but not caused by it. Spinodal points appear in both the Pöschl-Teller and the Lennard-Jones interactions, see Figs. 19 and 30, but neutron matter is far from such an instability, see Fig. 9. The cause of the instability rather seems to be a divergence in the ladder equation.

Monte Carlo simulations for a square-well interaction seem to indicate that the equation of state comes close, but do not reach, a spinodal decomposition before the BEC-BCS crossover is reached [150]; the data of Refs. 151, 125 are too sparse to make a conclusive statement.

Two points should, of course, be emphasized: One is that the “fixed-node” wave function displays the spurious instability discussed extensively above. Simulations based on such a wave function will at their best converge to the infimum of the energy that can be reached by wave functions of the type (5.2). Another problem arises since an exact simulation should converge towards the exact ground state, which could also be a droplet or a film covering the walls of the simulation box. However, we are interested here in the features of the uniform, but metastable state that can be present close to the spinodal point. We believe that this problem deserves a very careful examination.

The most serious problem in the application to nuclear systems is that there is evidence that parquet-class of diagrams is not sufficient for reliable calculations [152]. Unfortunately, this is insufficient when the core or the interaction is very different in the spin-singlet and the spin-triplet channel. In that case, the so-called “twisted chain” diagrams can be very important. We must keep this problem in mind but leave it for future work. In the language parquet theory, the “rungs” cannot be considered to be local interactions, but all time-orderings must be kept. In the language of the Jastrow-Feenberg method, the so-called commutator diagrams are important. The effect can be very dramatic. We will examine this in future work on neutron matter.

Appendix A Cluster expansions for the generating function. Diagonal terms

We show in the following series of appendices some details of the calculation of the generating functions given in Eqs. (5.12) and (5.13). First, we calculate the cluster expansion of the $G_{\rm diag}(\beta)$: Define

[TABLE]

where $(\Delta G_{{\bf m},{\bf m}})_{n}^{(s)}$ is the sum of all cluster contribution containing $n$ points and $s$ correlation lines $h(r_{ij})$. The individual terms of the cluster expansion are identical to the expansion of the normal ground state if we re-interpret the exchange line $\ell(rk_{\rm F})$ as

[TABLE]

Next, we formulate the resulting matrix elements as matrix elements of a diagonal second quantized operator,

[TABLE]

where

[TABLE]

where ${\cal P}(m_{1},\ldots,m_{n})$ stands for some permutation of the quantum numbers $(m_{1},\ldots,m_{n})$, and $G_{n}^{(s)}({\bf r}_{1},\ldots,{\bf r}_{n})$ is a combination of $s$ correlation lines connecting $n$ points. The leading order terms are

[TABLE]

Finally, we calculate the expectation value of these operators with respect to \bigl{|}{\mathrm{BCS}}\bigr{\rangle} by the usual contraction technique [153, 110]

[TABLE]

where $\bar{m}$, $\bar{n}$ refers the states with opposite momentum and spin of the states $m$, $n$, respectively, and $\operatorname{sgn}(m)=1$ ($-1$) when the spin of $m$ is up (down).

[TABLE]

and

[TABLE]

The diagrammatic representation of the two-body contribution (A.10) is shown in the diagrams (a) and (b) in Fig. 6 whereas $\left(\Delta G_{\rm diag}\right)_{3}^{(2)}$ corresponds to twice diagram (f) and diagram (d) of that figure.

We have to address a somewhat subtle issue connected with the appearance of “equivalent” diagrams which have a different topological structure, but the same value. Diagrams (d) and (e) of Fig. 2 are the simplest example. In the normal system it is easily seen that these diagrams have the same value because the exchange function $\ell(rk_{\rm F})$ satisfies the “convolution property:

[TABLE]

In contrast, as pointed already out in Ref. 110, the convolution property does not hold for the $\ell_{v}(r_{ij})$ exchange lines. Superficially one might therefore expect to get two different expressions for diagrams 4 and 5. An immediate consequence would also be that the expansion becomes point-reducible [110]. We therefore examine the issue here carefully.

We go back to the expansion (A.1), represented diagrammatically in Fig. 2. Diagram (e) can be written as

[TABLE]

The last equation holds because $\ell_{{\bf m}}(r_{ij})$ also satisfies the convolution property (A.12). The latter representation has then to be written in the second quantized form (A.3) which gives the correct weight factor of diagram (f) in Fig. 6.

The reason why these diagrams must be collected before calculating the expectation value with respect to the \bigl{|}{\mathrm{BCS}}\bigr{\rangle} is that we have two disconnected correlation lines $h(r_{12})h(r_{34})$ in $(\Delta G_{{\bf m},{\bf m}})_{4}^{(2)}$ leads additional momentum conservation:

[TABLE]

which would give a zero contribution to the second-quantized form. A very similar consideration was necessary for the calculation of the CBF single particle energies (3.69), see Ref. 23.

Appendix B Cluster expansions for the generating function. Off-diagonal terms

To derive cluster expansions for the generating function $G_{\rm offd}-\ln I_{\rm CBCS}$ we can use again the power-series method explained in Section 3.2 and define the normalization of $|\mathrm{CBCS}\rangle$ as

[TABLE]

where $N_{{\bf m},{\bf n}}^{(N)}(\alpha)$ is given by

[TABLE]

with

[TABLE]

see Eq. 3.67. We now use the power-series method [71] described in Section 3.2 to obtain a cluster expansion

[TABLE]

We can then write a term with $n$ correlation lines as

[TABLE]

where the $(\Delta N)_{m}^{(n)}$ are $m$-body clusters with $n$ lines. We then write these as off-diagonal matrix elements of second quantized operators

[TABLE]

A few leading order terms, where the states \bigl{|}{{\bf m}}\bigr{\rangle} and \bigl{|}{{\bf n}}\bigr{\rangle} differ by exactly two states, are

[TABLE]

In this formulation, we can calculate the expectation values of these operators with respect to \bigl{|}{\mathrm{BCS}}\bigr{\rangle}. The diagrams contributing to $G_{\rm offd}$ containing one or two correlation lines are shown in Fig. 33.

As a matter of course we have to pay again attention to the appearance of equivalent diagrams which must be combined before the BCS expectation value is calculated. Since the operations are identical to those outlined for the diagonal components in A, we skip this step.

Appendix C Cluster expansions for the energy numerator terms

In this appendix, we show details of the calculation of the energy numerator corrections $E_{\rm enum}$ defined in Eq. (5.11). We begin with the definition

[TABLE]

We can again use the power-series method to generate a cluster-expansion for $E_{\rm enum}$. Unlike above, we expand the diagonal matrix elements $H_{{\bf m},{\bf m}}^{(N)}$ and the off-diagonal matrix elements $N_{{\bf m},{\bf n}}^{(N)}$ individually. Analogously to A, we can write the diagonal matrix elements of the Hamiltonian in terms of a sum of diagonal, second-quantized $n$-body operators

[TABLE]

Similarly, we write the matrix elements $N_{{\bf m},{\bf n}}^{(N)}$ in the form (B.6). Then, each single contribution to the sum in the numerator of Eq. (C.1) has the form

[TABLE]

The expectation values of the creation and destruction operators with respect to \bigl{|}{\mathrm{BCS}}\bigr{\rangle} are then done with the usual contraction technique [153]. Unlike above, we now have two groups of creation and annihilation operators: The group $\{k_{1},\ldots,k_{n}\}$ belonging to the diagonal matrix elements, and the group $\{i_{1},\ldots,i_{m},j_{1},\ldots,j_{m}\}$ belonging to the $(\Delta\hat{N})_{m}^{(t)}$. We now classify the expansion in terms of the number of contractions between the two groups. The term where only operators within each group are contracted with each other cancels the second term in the above equation.

Hence, the simplest contribution is the set of diagrams where one pair of creation/destruction operators of the set $\{k_{1},\ldots,k_{n}\}$ is contracted with one pair of the group $\{i_{1},\ldots,i_{m},j_{1},\ldots,j_{m}\}$. All other contractions of operators of the set $\{k_{1},\ldots,k_{n}\}$ are done just as for the $G_{\rm diag\/}$ and the resulting $E_{\rm diag\/}$. It is clear that the resulting quantity can depend only on the quantum number of the state that has been contracted with the group $\{i_{1},\ldots,i_{m},j_{1},\ldots,j_{m}\}$. Hence, we can think of these contributions as a one-body operator

[TABLE]

where

[TABLE]

and $[\ldots]^{\cal F}$ denotes the Fourier transform defined as (2.10), and $k$ is the momentum carried by the common state $n_{k}$.

The leading term in $\epsilon^{(v)}(k)$ is kinetic energy $t(k)$, higher other correction terms to $t(k)$ are shown in Fig. 34. They reduce to the CBF single-particle spectrum in the normal system if one replace $\ell_{v}(r_{ij})$ by $\ell(r_{ij})$.

We can now rewrite $E_{\rm enum}^{(1)}$ as

[TABLE]

where

[TABLE]

and $\cal N$ can again be expanded in a cluster expansion as in B:

[TABLE]

A few leading order terms in $\hat{\cal N}$ are given by

[TABLE]

Thus, we can define

[TABLE]

We give details of the calculation here only for two typical terms: The two-body term $(\Delta E_{\rm enum})_{2}^{(1)}$ is

[TABLE]

The rest of the calculation is performed by the usual contraction technique which leads to

[TABLE]

This leads to diagrams (a), (b), and (c) shown in Fig. 35, where the black and red double-arrow respectively stand for a kinetic exchange lines,

[TABLE]

A new structure comes it at the next order, $(\Delta E_{\rm enum}^{(1)})_{3}^{(2)}$:

[TABLE]

We look only at the first term in $(\Delta\hat{N})_{3}^{(2)}$ in Eq. (C.10) which is evaluated as

[TABLE]

The first two terms are represented by diagrams (d), (e), and (f) in Fig. 35 and the last term is represented by diagram (g) in Fig. 35.

The diagrams we obtained in Fig. 35 can be further simplified by introducing the modified exchange lines $\ell_{v}^{\prime}(r_{ij})$ and $\ell_{u}^{\prime}(r_{ij})$ defined in Eqs. (5.20) and (5.19). We represent $\ell_{v}^{\prime}(r_{ij})$ and $\ell_{u}^{\prime}(r_{ij})$ by black and red double arrowed lines as shown in Fig. 36. For instance, the diagrams (a), (b), and (c) in Fig. 35 can be combined as the first diagram shown in Fig. 37.

To summarize, we can write the energy numerator terms containing one contraction the $\{k_{n}\}$ and the $\{i_{n},j_{n}\}$ groups in the closed form

[TABLE]

The two variational derivatives with respect to $\ell_{u}(r)$ and $\ell_{u}(r)$ can be done with the same diagrammatic techniques as the ones used in 23 for the calculation of the single-particle energies. It leads to the co-called “cyclic chain” diagrams already mentioned in Section 3.2. We have in the further analysis obtained only the leading terms which have been spelled out in Eq. (5.31).

The procedure can be easily extended for the case of two, three, …contractions between the $\{k_{n}\}$ and the $\{i_{n},j_{n}\}$ groups. Instead of a diagonal one-body operator $\hat{\epsilon}^{(v)}$, one obtains two, three, …operators in Eq. (5.31). We will not pursue this further because we did not include these terms in our calculations, and they disappear in the weakly coupled limit whereas the terms examined here lead to the energy numerator terms in Eq. (5.6).

The same consideration can be used to calculate the correction to the superfluid density due to correlations, note that we have pointed out in connection with Eq. 5.15 that the quantity defined there is not the same as the density

[TABLE]

where $\hat{N}$ is the density operator. We can write

[TABLE]

where N_{0}=\bigl{\langle}{\mathrm{BCS}}\bigr{|}\hat{N}\bigl{|}{\mathrm{BCS}}\bigr{\rangle}. Since the particle number operator $\hat{N}$ commutes with the correlation operator $F$, the second term in Eq. (C.19) has exactly the same structure as Eq. (C.1). We can then go through the further steps of the calculation, with the simplification that the number operator is a constant one-body operator. The result is therefore immediately obtained by replacing the one-body operator $\hat{\epsilon}^{(v)}$ by the number operator. Since this is the same as taking just the constant term $\mu$, we can immediately conclude that

[TABLE]

with

[TABLE]

Appendix D Calculation of exchange diagrams

A working formula for these exchange diagrams is [154]

[TABLE]

$I(k;r)$ is conveniently calculated by an expansion in spherical harmonics:

[TABLE]

which gives

[TABLE]

with

[TABLE]

Here

[TABLE]

With that, we get

[TABLE]

and

[TABLE]

Acknowledgment

This work was supported, in part, by the College of Arts and Sciences, University at Buffalo SUNY, and the Austrian Science Fund project I602 (to EK). We would like to thank Grigori Astrakharchik, Jordi Boronat, John Clark and Peter Schuck for useful discussions on the subject of this paper.