Bogolubov-Hartree-Fock theory for strongly interacting fermions in the low density limit

We consider the Bogolubov-Hartree-Fock functional for a fermionic many-body system with two-body interactions. For suitable interaction potentials that have a strong enough attractive tail in order to allow for two-body bound states, but are otherwise sufficiently repulsive to guarantee stability of the system, we show that in the low-density limit the ground state of this model consists of a Bose-Einstein condensate of fermion pairs. The latter can be described by means of the Gross-Pitaevskii energy functional.


Introduction
We consider a gas of fermions confined in an external trap at zero temperature. The particles interact through a two-body potential V which admits a negative energy bound state. At zero temperature and low particle densities, this leads to the formation of diatomic molecules forming a Bose-Einstein condensate. It was realized in the '80s [14,17] that BCS theory can be adequately applied to such types of tightly bound fermions. It was pointed out in [20,4,18,19] that in the low density limit the macroscopic variations in the pair density should be well captured by the Gross-Pitaevskii (GP) equation. From a mathematical point of view, the emergence of the GP functional in the low density limit was recently proven in [12] for the static case, and the dynamical case was subsequently treated in [10]. The assumption, that the two-body interaction potential allows for a bound state plays a crucial role. In the case of weak coupling where the potential is not strong enough to form a bound state, the pairing mechanism may still play an important role for the macroscopic behavior of the system, but the separation of paired particles can be much larger, in this case, than the average particle spacing. In fact this is the case in the usual BCS description of superconducting materials. Close to the critical temperature the macroscopic variation of the pairs is captured by the Ginzburg-Landau equation in this case, as pointed out by Gorkov [7] soon after the introduction of BCS theory.
The first mathematical proof of the emergence of Ginzburg-Landau theory from BCS theory was recently given in [6], which itself relied on earlier work on the BCS functional [8,11,5].
In the current paper our starting point is the full BCS Hartree-Fock functional. That is, we include the direct and exchange terms in the interaction energy. One also finds this functional under the name Bogolubov-Hartree-Fock (BHF) functional in the literature. The inclusion of the densitydensity interaction adds additional difficulties concerning stability of the system. It forces us to restrict to systems with a two-body potential V that, on the one hand, has an attractive tail deep enough to allows for a bound state and, on the other hand, is sufficiently repulsive at short distances to guarantee stability. This is consistent with typically considered interaction potentials [14]. We shall investigate the ground state energy of the BHF functional in the low density limit. We introduce a small parameter h playing the role of the inverse particle number, i.e., N = h −1 . We consider external potentials that confine the particles on a length scale of order h −1 , while the range of the interaction among the particles is of order one. This implies that the density of particles is of the order h 2 . Hence the small parameter h represents the square root of the density as well the ratio between the microscopic and macroscopic length scale. We are going to show that in the low density limit the fermions group together in pairs, such that the leading order in the energy is given by the number of pairs, 1/(2h), times the binding energy of a pair, −E b . The next to leading order is given by the energy of a repulsive Bose gas, consisting of fermion pairs, in a trap, and can be described in terms of the Gross-Pitaevskii energy functional. More precisely, if E BHF (h) denotes the BHF energy of 1/h fermions we shall obtain for small h, where E GP (g) denotes the Gross-Pitaevskii energy with appropriate interaction parameter g, which can be computed in terms of the microscopic quantities. The prefactor h/2 should be interpreted as N bos /L 2 , where N bos = N/2 = 1/(2h) is the number of fermion pairs and L = 1/h is the macroscopic length scale. We will also give a detailed description of the corresponding ground state of the BHF functional. Its minimizer turns out to be given, to leading order in h, in the form of the two particle wavefunction where α 0 is the ground state of a bound fermion pair with energy −E b , and ψ solves the GP equation and describes the density fluctuations of the pairs.
Our work is an extension of [12] in two directions. First, we include exchange and direct terms in the energy functional. Second, we avoid working with infinite, periodic systems, which allows us to significantly simplify the proof and also to improve the error bounds, utilizing ideas in [10]. In particular, we do not need to use here the rather involved semiclassical estimates of [6].
Our work presents the first proof of the occurrence of pairing in the ground state of a nontranslation invariant Bogolubov-Hartree-Fock system. (For a translation invariant system this was previously shown in [3].) The ground state properties of the BHF functional, in the context of Newtonian interaction, were studied in [15], see also [1]. Still it could not be shown that the fermions in the ground state exhibit pairing. Its occurrence was only shown numerically in [16]. In the low density limit, which we are studying here, the ground state actually predominately consists of pairs, in a sense to be made precise below. In particular, it is essential for our results that the pairing term is included in the energy functional; the Hartree-Fock functional for particle-number conserving states would lead to markedly different results, and is inappropriate for the description of low density gases.

Main Results
As in BCS theory, the state of a fermionic system is described by a self-adjoint operator Γ ∈ L L 2 (R 3 ) ⊕ L 2 (R 3 ) , satisfying 0 ≤ Γ ≤ 1. It is determined by two operators γ, α ∈ L L 2 (R 3 ) and has the form where 0 ≤ γ ≤ 1 is trace class and α is Hilbert-Schmidt and symmetric, i.e. α(x, y) = α(y, x), which implies that α * =ᾱ. We denote by γ, α the operators with kernels γ(x, y) and α(x, y), respectively. We note that we do not include spin variables here, but rather assume SU (2)-invariance of the states [13]. The full, spin-dependent Cooper-pair wave function is the product of α with an anti-symmetric spin singlet. Since α is symmetric, the latter is thus anti-symmetric, as appropriate for fermions. Given an external potential W and a two-particle interaction potential V , the corresponding Bogolubov-Hartree-Fock functional (BHF) is given by (2.1) We note that the terms in the first line represent the BCS functional, while the last line contains the additional exchange and direct terms in the interaction energy. A formal derivation of this functional from quantum mechanics can be obtained via restriction to quasi-free states, see [2], [8,Appendix] or [13]. Let us mention that our methods also allow to include a magnetic external vector potential, but for simplicity we shall not do so here. We study a system of h −1 fermions interacting by means of a two-body interaction V = V (x − y), confined in an external potential of the form W (hx). I.e., the external potential varies on a scale of order 1/h whereas V varies on a scale of order one. Since the trap W confines the particles within a volume of order 1/h 3 , the particle density is of the order h 2 . Hence the limit of small h corresponds to a dilute or low density limit.
Since we expect the interaction energy per particle pair to be of the order of the density, we shall also consider suitably weak external potentials, i.e., we replace W by h 2 W . It is convenient to use macroscopic variables instead of microscopic ones, i.e., we define x h = hx, y h = hy, α h (x, y) = The resulting BHF functional is then given by (now dropping the subscripts h) The corresponding ground state energy is denoted as For ψ ∈ H 1 (R 3 ) the GP functional is defined as The factors 1/2 and 2, respectively, in the first two terms result from the fact that (2.4) describes fermion pairs. The interaction parameter g > 0 will be determined by the BHF functional and represents the interaction strength among different pairs. We denote the ground state energy of the GP functional as We shall consider the minimization problem (2.3) and show that its value in the limit h → 0 is to leading order given by the binding energy of the fermion pairs, i.e. −E b 1 2h . This assumes, of course, that the two-body interaction potential V allows for a negative energy bound state, which is part of the following assumption.
and −2∆ + V has a normalized ground state α 0 with corresponding ground state energy −E b < 0.
Including direct and exchange term into the BCS functional gives rise to a new problem. A priori it is not clear whether the functional guarantees stability of the second kind. To ensure it we impose the following further assumption on V .
In other words, we consider potentials which, after cutting its positive part in half, can be bounded from below by functions with a non-negative Fourier transform. In particular, this means that the potentials have to have a strong enough repulsive core and a relatively small attractive tail, which still has to be large enough to allow for bound states. Remark 1. The following construction shows that it is easy to find potentials V with the desired properties of Assumptions 1 and 2: Choose a potential U which is strictly negative on an open set Ω ⊂ R 3 , such that U ≥ 0. The latter property can be ensured, e.g., by taking U to be the convolution of some function u with its reflection u(− · ). Now set V = 2U + − U − . Obviously this V fulfills Assumption 2. Finally, scale V according to V → λV until the negative part is deep enough for a bound state to appear.
With these assumptions we are ready to formulate our main theorem.
Under Assumptions 1 and 2, we have for small h, where g > 0 is given by Moreover, if Γ is an approximate minimizer of E BHF , in the sense that for some > 0, then the corresponding α can be decomposed as and ψ is an approximate minimizer of E GP in the sense that Remark 2. In contrast to the case of the usual BCS functional [12,10], where the coupling constant g only consists only of the BCS term it receives here two additional contributions from the direct and exchange energies, respectively. It is easy to see that our Assumption 2 implies that g dir + g ex ≥ 0, hence g > 0.
Remark 3. The proof of Thm. 1 partly relies on ideas in [10], where the corresponding time-dependent problem was studied for the BCS functional, i.e., in the absence of direct and exchange term. A similar result can also be shown to hold in the case of the time-dependent BHF equation, which in a different context was studied in [9]. By following the strategy of [10] and handling the exchange and direct terms in a similar way as done here, one can derive the time-dependent GP equation with interaction parameter g. Notation: In the following we often write a b to denote a ≤ Cb for some generic constant C > 0.

Stability
Before giving a sketch of the proof of Theorem 1 we show how Assumption 2 gives rise to stability of the second kind. In fact we simply show that the assumption guarantees that the sum of the direct and exchange terms is non-negative. To this aim we first consider the exchange term and estimate using |γ(x, y)| 2 ≤ γ(x, x)γ(y, y). Hence we have for the sum of direct and exchange term where we used the assumption V − 1 2 V + ≥ U . Since U ≥ 0 the last term is non-negative. Hence the question of stability is reduced to the corresponding problem for the BCS functional, and is easily seen to hold under our assumptions on V .

Sketch of the proof of Theorem 1
The proof of (2.6) consists of deriving appropriate upper and lower bounds on E BHF (h).

Upper bound
For the upper bound one has to construct a suitable trial state. We shall proceed similarly to [10] and define the trial state Γ ψ via the pair wavefunction Since we expect that the system in its ground state energy consists predominantly of pairs we define the one particle density γ ψ such that to leading order it equals α ψ α ψ . More precisely, we choose the trial state The function ψ here is only approximately normalized, i.e., ψ 2 = 1 + O(h 2 ), to ensure that Tr γ ψ = 1/h. We will see below that for small enough h the definition (4.3) guarantees that 0 ≤ Γ ψ ≤ 1.
In the limit of small h the GP energy functional emerges from the BHF functional E BHF (Γ ψ ) as follows. If we consider the kinetic energy term plus the pairing term and subtract the total binding energy, 2 Tr γ ψ , the contribution to coming from the α ψ α ψ term in (4.3) can be written as Since α ψ (x, y) is symmetric we can replace ∆ x by 1 2 (∆ x +∆ y ). In terms of center of mass X = (x+y)/2 and relative coordinates r = x − y the kinetic energy has the form ∆ x + ∆ y = 1 2 ∆ X + 2∆ r , such that in these coordinates the last term has the form where we used the fact that α 0 is the normalized ground state of −∆ term is due to the contribution of α ψ α ψ in the direct and exchange interaction terms. The estimation of these terms is straightforward but tedious and occupies the main part of the proof.
Furthermore, it will be easy to show that Consequently we shall obtain Finally, we remark that the constraint Tr γ ψ = 1/h implies for ψ that ψ 2 The precise derivation of this bound will be given in Section 6.
• Since the infimum of E BHF is attained by a projection [2], it would be natural to chose the trial state Γ ψ as a projection. The operator γ ψ would then satisfy γ ψ = γ 2 ψ + α ψ α ψ . The expansion of γ ψ in terms of α ψ α ψ would be more complicated, however, and we find the choice (4.3) more convenient.
• Our trial state satisfies 0 ≤ Γ ψ ≤ 1 for small enough h. To see this note that the condition is If γ ψ is of the special form (4.3), which is a function of α ψ α ψ , the off-diagonals of vanish and thus the statement is equivalent to Plugging in the expression (4.3) for γ ψ (4.8) is equivalent to In Corollary 1 below we shall show that the operator norm of α ψ satisfies α ψ ∞ h 1/2 , which guarantees that (4.9) is satisfied for h small enough. In fact, h 1/2 in (4.3) could be replaced by any factor large compared to h, but a different choice would not improve our error bounds.

Lower bound
From the upper bound we learn that for an approximate ground state Γ we can assume We will show in Lemma 3 that the corresponding α necessarily has to be of the form for an appropriate ψ ∈ H 1 (R 3 ), with ξ being small compared to α ψ , i.e., The function ψ is obtained by projecting α in the direction of α 0 with respect to the relative coordinates, We shall show that With this ψ at hand one can then define Γ ψ as in (4.2)-(4.3). With the help of the decomposition (4.10) one then argues that the difference between E BHF (Γ) and E BHF (Γ ψ ) is bounded from below by a term of higher order than the contribution from the GP functional. More precisely, This estimate is uniform in ψ, since it is possible to obtain a priori bounds on the H 1 -norm of ψ that are independent of h. Using now our calculation (4.6) from the upper bound immediately implies Together with (4.7) this combines to (2.6).

Useful properties of the pair-wavefunction
In the following we shall derive some useful properties of the pair wavefunction (4.3), which will be used throughout the proof. Recall that α 0 was defined in Assumption 1 to be the normalized ground state of −2∆ + V . It is a rapidly decaying function, and both |α 0 | and |∇α 0 | have smooth Fourier transforms which are in L p (R 3 ) for any p ≥ 2.
(i) For n ∈ {2, 4, 6}, α ψ n n h n−3 ψ n n |α 0 | n n , (5.1a) (iv) Let σ be a Hilbert-Schmidt operator. Then Let us mention that we use the symbol · p for the L p -norm of functions as well as for the operator norm in the corresponding Schatten class. E.g., the left side of (5.1a) concerns Schatten norms, while on the right side the norms are in L n (R 3 ).
Proof of Lemma 1, Part I. We postpone the proof of (5.1a), (5.1b) and (5.2) to the appendix. In order to see (5.3) we use the Hölder and Sobolev inequalities as Similarly, Proof. The estimates (5.5a)-(5.5c) are a consequence of (5.1a) and (5.1b). In the case of n = 6, we use the Sobolev's inequality and in the case of n = 4, we use Hölder combined with Sobolev to conclude 2 . Inequality (5.5d) follows immediately from α ψ ∞ ≤ α ψ 6 together with (5.5b).
Remark 5. Since γ ψ is to leading order equal to α ψ α ψ , we obtain as a corollary that the operator norm of γ ψ is at most O(h), meaning that the largest eigenvalue of the one-particle density matrix is of order h. However, the two-body density matrix corresponding to the state Γ ψ is to leading order of the form |α ψ α ψ |, and hence has one large eigenvalue of order h −1 . This is a manifestation of the Bose-Einstein condensation of the fermion pairs.
The desired upper bound (4.7) is then an immediate consequence of the following estimates: where the constants g bcs , g ex and g dir are defined in Remark 2. The remainder of this section will be devoted to the proof of these estimates.
(6.2) From the fundamental theorem of calculus we obtain Using the Cauchy-Schwarz inequality for the integration over the X variable, the last integral is bounded by Since α 0 is the ground state of the Schrödinger operator −2∆ + V and hence rapidly decaying, | · |α 0 2 is finite. This shows (6.1b).
6.3 Direct and exchange term (Proof of (6.1c) and (6.1d)) We first argue that the leading order contribution of the direct and exchange terms originates from replacing γ ψ by α ψ α ψ . To see this, we simply estimate the differences (6.3b) Both expressions can be bounded using the following lemma, whose proof is elementary.
For (6.4b) we follow a similar strategy and first split Applying to σ, δ, and σ + δ the fact that for positive trace class operators a its kernel satisfies together with the Cauchy-Schwarz inequality, we obtain the stated inequality.
In order to recover the ψ 4 4 contribution we inspect the remaining parts of the direct and the exchange terms separately. We begin with the exchange term and write explicitly Introducing new variables and rescaling r/h → r, s/h → s, t/h → t, the last expression becomes The latter equals

This can be bounded by
using the Hölder, Sobolev and Cauchy-Schwarz inequalities. This shows (6.1c). We continue with the direct term. Its remaining part is given by where we changed to the variables and rescaled r, s, t. By proceeding as above, we see that this expression equals This shows (6.1d), and thus concludes the proof of the upper bound.

Lower bound
Our proof of the lower bound on E BHF (h) in Theorem 1 consists of two parts. As a first step we obtain a priori bounds on approximate ground states.

Then these functions satisfy the bounds
Note that our definition implies that ξ(X, · ) is orthogonal to α 0 ( · /h) for almost all X. The norms in (7.2e)-(7.2g) are in L 2 (R 6 ).
Proof. We have seen in Section 3 that the sum of the direct and exchange terms is non-negative. Consequently, We bring the term h 2 W ∞ Tr(γ) h to the left side. Adding and subtracting the expression The last two terms on the right side can be expressed via center-of-mass and relative coordinates as where we used that α 0 is the normalized zero energy eigenvector of the operator −∆ + V /2 + E b /2, as well as the fact that ξ(X, ·) is orthogonal to α 0 (·/h) for almost every X ∈ R 3 . Hence (7.3) implies Since all terms on the right side are non-negative, and γ − γ 2 ≥ αᾱ, we immediately obtain the estimates (7.2a), (7.2b), (7.2d) and (7.2f). According to Assumption 1 the operator −∆ + V /2 has a spectral gap between the ground state energy −E b /2 and the next eigenvalue. This implies that we can find a κ > 0 and an ε > 0 such that In combination with (7.5) this yields the estimates (7.2e) and (7.2g).
With the aid of the function ψ we can define a corresponding state Γ ψ as in (4.1)-(4.3). By multiplying ψ with a factor λ = 1 + O(h 2 ) we can assume that Tr γ λψ = 1/h. The second step now consists of proving that for a lower bound we can replace E BHF (Γ) by E BHF (Γ λψ ) up to higher order terms. Together with the calculations from the upper bound this implies the lower bound stated in Theorem 1.
It remains to prove the bound (7.8), which is an immediate consequence of the following estimates: (7.9d) The remainder of this section will be dedicated to proving these estimates. 7.1 Kinetic and potential energy (Proof of (7.9a)) Let us decompose γ according to where (γ − αα − γ 2 ) and (γ − αα) 2 are positive self-adjoint operators and thus Adding and subtracting the term we obtain The identity (7.4) immediately implies that and hence the sum of the second and third lines on the right side of (7.10) is bounded from below by . Hence the proof of (7.9a) reduces to establishing the estimates which we are going to show in the following. Inequality (7.14) is an immediate consequence of (5.2). It also implies that it is enough to show (7.12) for λ = 1. To do this, we rewrite αααα − α ψ α ψ α ψ α ψ = α ψ ααξ + ξααα ψ + ξααξ + α ψ (αα − α ψ α ψ )α ψ . (7.15) With H := −h 2 ∆ + E b 2 we obtain with the aid of Hölder's inequality for traces Note that for any operator T , we have where in the last line the operators ∇ X T and ∇ r T are defined via the kernels (∇ X T )(x, y) and (∇ r T )(x, y), respectively. By applying this to the terms in (7.16) we obtain where we used ∇ X α ψ 6 = α ∇ψ 6 ≤ α ∇ψ 2 = ∇ψ 2 h −1/2 , together with (5.5b), (5.5c), (7.2e), (7.2f) and (7.2g). The term αα − α ψ α ψ 3/2 in (7.16) can be bounded by where we used (5.5b) and (7.2e). By combining (7.16) with (7.17)-(7.19) we obtain (7.12).
To show (7.13), we can bound The first factor on the right side is bounded by O(h) according to (7.2a). Moreover, which is bounded by O(h 1/2 ) using (7.12) together with (5.2). This proves (7.13).
7.3 Direct and exchange term (Proof of (7.9c) and (7.9d)) Our strategy is as follows. As a first step we reduce the direct term and exchange term to corresponding expressions involving α only, and show that As a second step, we show that up to an error O(h 3/2 ) we are able to replace α by α ψ in the corresponding expressions, i.e., These two steps together, in combination with λ = 1 + O(h 2 ), lead to (7.9c) and (7.9d), respectively. The estimates (7.21) and (7.22) can be obtained by applying Lemma 2 with σ = αα and δ = γ − αα. As a result we obtain that the left sides of both (7.21) and (7.22) are bounded by (7.25b) By (7.2b), the term (7.25a) is bounded by For (7.25b), we are going to use the decomposition α = α ψ + ξ in the form αᾱ = α ψ α ψ + ξα ψ + α ψ ξ + ξξ and we thus have to bound four terms. First, observe that Second, using (5.3), For the remaining terms we use (5.4) with σ = ξ and the Cauchy-Schwarz inequality to obtain We now turn to the estimates (7.23) and (7.24). The difference of the exchange terms in (7.23) is bounded by The 2-norm of αα − α ψ α ψ can be bounded by the 3/2-norm, which in turn is bounded by O(h) according to (7.19). Moreover, α ψ α ψ 2 = α ψ 2 4 h 1/2 by (5.1a), proving (7.23). For the direct term we insert the decomposition α = α ψ + ξ into the difference in (7.24), yielding 15 terms. However, due to symmetry, it suffices to estimate the following 5 terms We begin with (7.27a). Obviously For (7.27b) we obtain with the help of (5.3) For the last three terms we invoke equation (5.4) from Lemma 1 with σ = ξ and the Cauchy-Schwarz inequality. For (7.27c) this gives The desired bound O(h 3/2 ) then follows from the fact that the last factor α ψ α ψ (·, ·) 2 on the right side is of order O(h −1 ). To see this, we write Changing to the variables r = x − y, s = x − z and x and using Cauchy-Schwarz in x, we indeed obtain α ψ α ψ (·, ·) 2 2 = h −8 This concludes the proof of (7.9c) and (7.9d).
A Appendix: Proof of Lemma 1 Proof of Lemma 1, Part II. We first prove (5.1a) and (5.1b). For n ∈ 2N, we can write We switch to the following coordinates X = 1 n n k=1 x k r k = x k+1 − x k , k = 1, . . . , n − 1.

(A.2)
It is easy to see that the corresponding Jacobi determinant is equal to 1. Moreover, we can recover the original coordinates via i.e.
Using integration by parts, we can bound A 1 as where the last inequality follows from α ψ ∞ ≤ α ψ 6 , which is h 1/2 ψ 6 as shown above.
The desired result then follows from the Sobolev inequality ψ 6 ∇ψ 2 . This concludes the proof of Lemma 1.