This paper is devoted to the study of an integrodifferential equation introduced in [Lieb 1963] in connection with the study of the Bose gas, a many-body problem in quantum mechanics. The equation is
(1-1) with x ∈ ޒ d and * denoting convolution: u * u(x) := u(xy)u(y) dy. Here, V is a given function (called the potential) in L 1 ޒ( d )∩ L p ޒ( d ), with p > d/2 for d 2 and p > 1 for d = 1. We assume V to be nonnegative. (This corresponds to a repulsive interaction between the particles in the underlying quantum system.) The two parameters e and ρ are nonnegative numbers, and they are related by a constraint, namely
(1-2)
We are interested in solutions of (1-1) that satisfy the constraint (1-2), or, in other words, solutions of the system (1-1)-(1-2). We are particularly interested in the case d = 3, though other dimensions are also of interest. As explained in [Lieb 1963], the parameter ρ corresponds to the particle density N /V of the underlying Bose gas in the large volume and large particle number limit, and e = E/N stands for the energy per particle.
One would like to fix a value ρ for the density, and then one expects, on the basis of the arguments in [Lieb 1963], that there will be a unique value of e = e(ρ) such that there is a solution of (1-1)-(1-2) with u taking values in [0,1]. This value of e is then the energy per particle of the Bose gas in its ground state.
The problem of determining this ground state energy per particle, as a function of the density, has attracted the attention of a great many researchers since the pioneering work [Lenz 1929]. In that paper and subsequent work [Bogolubov 1947;Lee et al. 1957], an asymptotic expansion of e(ρ) for d = 3 and small ρ was obtained:
where a, called the scattering length, is a property of the pair-interaction potential V(x) and is defined in (4-8)-(4-12) below. Here, we set both the mass m of the particle and Planck's constant h to 1. This early work was not mathematically rigorous, and it was not until [Lieb and Yngvason 1998] that the validity of the first term 2πρa was proved, and not until [Fournais and Solovej 2019] that the validity of second term was also proved, utilizing upper bounds proved earlier in [Dyson 1957;Yau and Yin 2009]. This timeline gives some idea of the complexity of the problem of directly studying the Bose gas ground state as a many-body problem. The complexity makes it very attractive to try to show that the system (1-1)-(1-2) provides a useful and illuminating route to the computation of the properties of the ground state for a Bose gas. Interest is piqued further by the fact that numerical studies show that the function e(ρ) computed using the system (1-1)-(1-2) is surprisingly accurate for all densities, not only low densities, as we discuss later in this paper. Until now, however, there has been no mathematically rigorous study of this system, and even the most basic questions concerning existence and uniqueness of solutions had remained open.
In this paper, we settle some of these basic questions and raise others. It may at first appear surprising that (1-1) poses any serious mathematical challenges. After all, if one replaced the convolution nonlinearity u * u in (1-1) by a power nonlinearity, say u 2 , one would have a familiar sort of local elliptic equation:
(1-4)
However, the convolution nonlinearity in (1-1) makes it nonlocal, and very different from (1-4).
As explained in [Lieb 1963] the solutions of physical interest are integrable and must satisfy u(x) 1 for all x. Our first result is that for integrable solutions of the system (1-1)-(1-2), the upper bound u 1 implies the lower bound u 0: Theorem 1.1 (positivity). Suppose that V is nonnegative and integrable and that u is an integrable solution of (1-1)-(1-2) such that u(x) 1 for all x. Then u(x) 0 for all x, and all such solutions have fairly slow decay at infinity in that they satisfy |x|u(x) d x = ∞.
(1-5)
Thus, any physical solutions of (1-1)-(1-2) must necessarily satisfy the pair of inequalities 0 u(x) 1 for all x.
(1-6) This a priori result, which we prove before we take up existence and uniqueness, relies on results [Carlen et al. 2020] obtained in collaboration with Michael Loss on the convolution inequality f f * f in L 1 ޒ( d ). While u(x) 1 is a physical requirement, u(x) 0 is not; see Section 6 for details.
The converse of Theorem 1.1 also holds, as stated in the following theorem.
) such that u(x) 0 for all x, then u(x) 1 for all x.
Remark. We have thus proved that u 0 if and only if u 1. This, in principle, leaves the door open to solutions that are sometimes > 1 and sometimes < 0, though we do not believe such solutions exist.
Before stating our main theorems, we make a few observations. 1A. The system (1-1)-(1-2) is actually equivalent to (1-1) and the constraint
(1-7)
To prove this, consider the operator
which is given by
where Y 4e is the Yukawa potential [Lieb and Loss 2001, Section 6.23], which is nonnegative and satisfies
(1-10)
Equation (1-1) can be rewritten as
(1-11)
Since u and V are assumed to be integrable, and u(x) is assumed to satisfy (1-6), all terms in (1-11) are integrable, and integrating yields
(1-12)
Thus, for integrable solutions u of (1-1) satisfying (1-6), the constraint (1-2) is equivalent to (1-7).
1B. There is another useful way to write the system (1-1)-(1-2). The damped heat semigroup e -t (-+4e) is a strongly continuous contraction semigroup on L p ޒ( d ), and the domain of its generator is D(-+4e) = W 2, p ޒ( d ). By the Sobolev embedding theorem [Lieb and Loss 2001, Theorem 10.2], since p > d/2, all functions f ∈ D(-+ 4e) are continuous and vanish at infinity. Since V 0, we know e -tV is also a strongly continuous contraction semigroup on L p ޒ( d ), and since V ∈ L p ޒ( d ), the domain of its generator, D(V), contains all bounded functions, and in particular W 2, p ޒ( d ). Writing V as the sum of a piece with a small norm in L p ޒ( d ) and another piece that is bounded, it is easy to see that there are numbers a, b > 0 with a < 1 2 such that for all f ∈ W 2, p ޒ( d ),
Then by the Banach space version of the Kato-Rellich theorem [Reed and Simon 1975, p. 244] the operator
Define K e to be the inverse operator
By the Trotter product formula, the operator K e has a positive kernel that we denote by K e (x, y); in particular, K e preserves positivity. By the resolvent identity
we conclude that 0 K e (x, y) G e (x, y) (1-16) for all x, y. Thus, the operator K e extends to a bounded operator on L 1 ޒ( d ) and all terms in the equation
are well-defined whenever u is integrable. Moreover, since V ∈ L p ޒ( d ), and since u * u ∈ L p ޒ( d ) when u is integrable and satisfies (1-6), every integrable solution u of (1-17) that satisfies (1-6) actually belongs to W 2, p ޒ( d ) and satisfies (1-1). Several simple bounds follow almost immediately from this form of the equation. First of all, since the last term on the right of (1-17) is nonnegative, we have an a priori lower bound on u(x), namely u(x) u 1 (x) := K e V(x).
(1-18)
Integrating both sides of (1-18), and using (1-7) yields an upper bound on ρ depending only on e, namely,
. By (1-2) and (1-18),
.
(1-20)
In fact, the left side of (1-20) is equal to half of the right side. To see this observe that u 1 = K d V satisfies (-+ 4e + V)u 1 = V, and hence u 1 = G e (V(1u 1 )). Integrating both sides yields u 1 d x = (1/(4e)) V(1u 1 ) d x. By (1-18), we obtain the simpler (albeit less sharp) bounds
In particular, this shows that the system (1-1)-(1-2) does not have a solution for arbitrary values of ρ and e: when either is small, a solution of the type we seek can only exist if the other is correspondingly small, as specified by (1-21) and (1-22). In fact, as is stated in the following theorem, ρ and e are constrained to be related by a functional equation.
Theorem 1.3 (existence and uniqueness).
Then there is a constructively defined continuous function ρ(e) on (0, ∞) such that lim e→0 ρ(e) = 0 and lim e→∞ ρ(e) = ∞ and such that for any e 0 and ρ = ρ(e), the system (1-1) and (1-2) has a unique integrable solution u(x) satisfying u(x) 1. Moreover, if ρ = ρ(e), the system (1-1) and (1-2) has no integrable solution u(x) satisfying (1-6).
Remarks. • We do not assume here that the potential is radially symmetric. However, the uniqueness statement implies that u is radially symmetric whenever V is radially symmetric.
• The function ρ(e) is the density function, which specifies the density as a function of the energy. Thus, our system together with (1-6) constrains the parameters e and ρ to be related by a strict functional relation ρ = ρ(e). In most of the early literature on the Bose gas, ρ is taken as the independent parameter, as suggested by (1-3): One puts N particles in a box of volume N /ρ, and seeks to find the ground state energy per particle, e, as a function of ρ. Our theorem goes in the other direction, with ρ specified as a function of e. We prove that e → ρ(e) is continuous, and we conjecture that ρ(e) is a strictly monotone increasing function. In that case, the functional relation could be inverted, and we would have a well-defined function e(ρ).
• Since lim e→0 ρ(e) = 0 and lim e→∞ ρ(e) = ∞, the continuity of e → ρ(e) implies that for each ρ ∈ (0, ∞) there is at least one e such that ρ(e) = ρ.
Having proved that the solution to the simple equation is unique, our second main result is an asymptotic expression for e(ρ), both for low and for high density.
Theorem 1.4 (asymptotics of the energy for d = 3). Consider the case d = 3. Let V be nonnegative, integrable and square-integrable. Then, for each ρ > 0 there is at least one e > 0 such that ρ = ρ(e). For any such ρ and e we have the following bounds for low and high density (i.e., small and large ρ). For low density,
where a is the scattering length of the potential, which is defined in (4-11). For high density, in any dimension d 1,
(1-24)
Remark. For low densities in d = 3, the energy e predicted by the simple equation (1-1)-(1-2) is asymptotically equal to the ground state energy of the Bose gas [Lee et al. 1957;Yau and Yin 2009;Fournais and Solovej 2019]. For high densities, when the potential has a nonnegative Fourier transform, the asymptotic formula for the ground state energy of the Bose gas coincides with (1-24) [Lieb 1963, Appendix]. Thus, the simple equation yields the same asymptotes for both low and high densities as the Bose gas does (at least when the potential has a nonnegative Fourier transform, as in the example V(x) = e -|x| discussed in Section 6B).
Theorem 1.5 (decay of u at infinity). In all dimensions, provided V is spherically symmetric with
Under stronger assumptions on the potential, this is actually the case.
, and, for |x| > R,
for some A, B > 0, then there exists α > 0 such that
(1-27)
Remarks.
• This result is consistent with a prediction in [Lee et al. 1957] that the truncated 2-point correlation function in the ground state of the Bose gas decays like |x| -4 .
• To prove this theorem, we will use analytical properties of the Fourier transform V of V, which is why we assume that V decays exponentially at infinity. For potentials with slower decay, it seems that the decay of u should still be |x| -4 , except if V itself decays slower than |x| -4 , in which case u should decay like V.
• It is presumably not too difficult to extend this result to cases with potentials that are not spherically symmetric.
Remark. The simple equation (1-1) is actually an approximation of a richer equation for u [Lieb 1963], which should more accurately depict the Bose gas; see (7-2). Little is known about this richer equation.
The paper is organized as follows. We prove Theorem 1.1 in Section 2, Theorems 1.2 and 1.3 in Section 3, Theorem 1.4 in Section 4, and Theorem 1.5 in Section 5. In Section 6, we explain how the simple equation is related to the Bose gas, and present some numerical evidence that it is very good at predicting the ground state energy. In Section 7 we discuss a few open problems and extensions.
As explained in the Introduction, the solutions of (1-1)-(1-2) that are of physical interest are those that are integrable and satisfy u(x) 1 for all x. In this section we prove, making no assumptions on the potential V other than its positivity and integrability, that all such solutions are nonnegative and have slow decay so that |x|u(x) d x = ∞.
Our starting point is the form of (1-1) given in (1-11). For an integrable solution u, define
and (1-11) can be written as
Lemma 2.1. Let u(x) be an integrable solution of the system (1-1)-(1-2) such that u(x) 1 for all x. Let f be defined in terms of u, e and ρ by (2-1). If f (c) 0 for all x, then u(x) 0 for all x.
Proof. Since Y 4e * (V(1u(x))) 0, it follows that
(2-4)
Integrating, we find u -d x 1 2 u -d x, and this implies u -= 0.
Proof of Theorem 1.1. Multiply (2-3) through by 2eρ, and then convolve both sides with Y 4e . The result
for all x. It is proved in [Carlen et al. 2020] that all integrable solutions of (2-5) are nonnegative and have integral no greater than 1 2 , and that moreover, (2-2) and (2-3) together imply
(2-6) However,
(2-7)
Then since lim x→∞ (4e|x| -1 Y 4e * |x|) = 1, (1-5) follows.
As was shown in (1-11) and (1-17), there are at least two ways to write (1-1) as a fixed-point equation.
As it turns out, only the latter one
is adapted to solution by iteration, because of its monotonicity properties. Starting with u 0 (x) = 0, define
for n 1. It is easy to see that for arbitrary e, ρ 0, this produces a monotone increasing sequence of nonnegative integrable functions. Thus, u(x) := lim n→∞ u n (x) will exist, but it need not be integrable and it need not satisfy (1-2) or (1-6).
To bring (1-2) into the iteration scheme, we take e as the independent parameter, and define a sequence {ρ n } along with the sequence {u n (x)}, both depending on e, through
and
.
(3-4)
Comparing (3-3) to (3-1), note that the analog of now depends on n.
Both sequences {ρ n } and {u n } are well-defined and increasing, and for all n,
(3-5)
Proof. We proceed by induction. By definition, u 0 = 0 and
. Also by definition
. As noted in the discussion between (1-20) and (1-21), 2
Since t → t -1 is monotone decreasing on (0, ∞), this shows that ρ 1 > ρ 0 , and that (3-5) holds for n = 1. Now suppose that u n u n-1 0, ρ n ρ n-1 0, and
all of which we have just verified for n = 1. Then
and thus
Then since
This proves (3-5) for n + 1 and shows that (3)(4)(5)(6)(7)(8)(9)(10)(11) and then, as before,
Then for all n, the function u n (x) is continuous, vanishing at infinity, and 0 u n (x) 1.
Since K e maps L p ޒ( d ) into W 2, p ޒ( d ), u 1 is continuous and vanishes at infinity. Let A := {x : u 1 (x) > 1}.
Then A is open. If A is nonempty, then u 1 is subharmonic on A, and hence takes on its maximum on the boundary of A. Since u 1 would equal 1 on the boundary, this is impossible, and A is empty. This proves the assertion for n = 1. Now make the inductive hypothesis that 0 u n (x) 1 for all x. Then
By Young's inequality, u n * u n p u n p u 1 1 , and hence
It follows as before that u n+1 is continuous and vanishing at infinity, and in particular, bounded, and
where we have used ρ n u n 1 1, which is valid on account of (3-5). Define A := {x : u n+1 (x) > 1}. Then u n+1 is subharmonic on A, and maximal on the boundary of A, where u n (x) would equal 1. This contradiction shows that u n+1 ∞ 1.
(3-13)
Then both limits exist, u ∈ W 2, p ޒ( d ) and u satisfies (1-1), (1-2) and (1-6).
Proof. By Lemma 3.1, both limits exist, and by (3-5), ρ(e)
. Also by Lemma 3.1, Lemma 3.4. For all e ∈ (0, ∞), the solution u of the system (1-1)-(1-2) that we have constructed by iteration in Lemma 3.3 is the unique nonnegative integrable solution for ρ = ρ(e). Moreover, there does not exist such any such solution when ρ = ρ(e).
Proof. Consider any nonnegative solution integrable ũ, with
We first show that ũ u n by induction. We have
(3-16)
Since u 0 = 0, the positivity of ũ implies the positivity of ũ(x)u 1 (x). If ũ u n-1 , then, by (3-4), ρ ρ n-1 , from which ũ u n follows easily. This proves that both ρ ρ and ũ u. However, integrating both sides of the latter inequality yields
.
(3-17)
Since ρ ρ, equality must hold, and then since ũ u, it must be that so u = ũ. In particular, for each ρ ∈ (0, ∞), there is at least one e ∈ (0, ∞) such that ρ = ρ(e).
Proof. We now turn to the continuity of e → ρ(e). We close this section by remarking that if V is radially symmetric, then so is u 1 = K e V, and then by a simple induction, so is u n , and hence also the unique solution u provided by Theorem 1.3. This is consistent with the first remark following Theorem 1.3.
In this section, we prove Theorem 1.4. We will first prove the high-density asymptote (1-24), and then proceed to the low-density (1-23).
By Theorem 1.3, for each ρ > 0 there exists at least one e such that ρ(e) = ρ. If there is more than one, the theorems proved in this section apply to every such solution. Throughout this section, let u ρ denote the solution provided by Theorem 1.3 and any such choice of e.
4A. High-density ρ.
Note that this is not an optimal bound, as follows from (1-20).
Proof. By (1-2), it suffices to prove that
and take the decomposition
which, by (1-7), is bounded as
completing the proof.
4B. Low-density ρ. In this section, we only consider the dimension d = 3. As before, we suppose that V ∈ L 1 ޒ( 3 ) ∩ L p ޒ( 3 ), p > 3 2 , and V 0. We first recall the definition of the scattering length of the potential V and relate it to the solution of the system (1-1)-(1-2). The scattering equation is defined as
Note that (4-8) can be written as (-+ V)ϕ = V, and hence the solution is
where u 1 is the first term of the iteration introduced in the previous section. It follows from Lemma 3.2 that 0 ϕ(x) 1 for all x.
We now impose a mild localization hypothesis on V: For R > 0 define V R (x) = V(x) for |x| > R and otherwise V R (x) = 0. We require that, for some q > 1 and all sufficiently large R, Lemma 4.2. Let V ∈ L 1 ޒ( 3 ) ∩ L p ޒ( 3 ), p > 3 2 , and suppose that the localization condition (4-10) is satisfied. Let ϕ be the corresponding scattering solution given by (4-9). Then the scattering length a := lim |x|→∞ ϕ(x) exists and satisfies
(4-12)
Proof. By the resolvent identity, ϕ(x) = G * (V(1ϕ))(x) where G(x) = 1/(4π|x|). Since p > 3 2 , p ′ < 3, and it is easy to decompose G into the sum of two pieces
(4-13) For 0 < r < 1 and |y| < r |x|,
It follows that for all sufficiently large |x|,
Taking |x| → ∞ and then r → 0 proves (4-10).
Remark. The following lemma is valid if the scattering length a is defined by (4-12). For this reason, we do not impose the additional condition (4-10) in the statement of Theorem 1.4: Lemma 4.2 reconciles the stated definition with the formula (4-12). (4-15)
Proof. The scheme of the proof is as follows. We first approximate the solution u by w, which is defined as the decaying solution of
The energy of w ρ is defined to be (4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17) and, as we will show, it is close to e; more precisely,
In addition, (4-16) is quite similar to the scattering equation (4-8). In fact we will show that e w is close to the energy 2πρa of the scattering equation
Summing (4-18) and (4-19), we find
from which (4-15) follows. We are thus left with proving (4-18) and (4-19).
Proof of (4-18): By (1-2) and (4-17),
(4-21)
We will work in Fourier space (4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22) which satisfies, by (1-1),
Therefore,
Similarly, the Fourier transform of w ρ is
Note that, as |k| → ∞, we have û ∼ 2eS(k)/(ρk 2 ), so, while ûρ is not integrable, ûρ -ŵρ is. We invert the Fourier transform:
We change variables to k := k/(2 √ e):
is monotone increasing. In addition, by (4-24) and (1-1), and using the fact that u ρ (x) 1 (see Lemma 3.2) and V(x) 0,
By dominated convergence, and using the fact that S(0
Using (1-22), this proves (4-18). Incidentally, again by dominated convergence,
Proof of (4-19): Let ξ(r ) := w ρ (r )ϕ(r ).
(4-36)
By (4-16), (4-8) and (1-1),
Therefore, by (4-12),
By (4-34),
is bounded, we can use dominated convergence to show (4-19).
In this section, we prove Theorem 1.5. Our proof assumes that V decays exponentially, because we will use analyticity properties of the Fourier transform of the potential V. In particular, the theorem holds if V has compact support. We expect the result to hold for any potential that decays faster than |x| -4 . Algebraic decay for u seems natural: by (1-1), u * u must decay at infinity in the same way as u. This is the case if u decays algebraically, but would not be so if, say, it decayed exponentially.
Proof of Theorem 1.5. We begin by proving (1-25) in arbitrary dimension. Recall that the first part has already been proved in Theorem 1.1 without the additional assumption on the potential. For the second part, recall that by the first remark after Theorem 1.3, u is also radial, and hence V(1u) is nonnegative and radial. It then follows from the hypotheses on V that g
(5-1)
Then, as explained in Section 2, if f := 2eρY 4e * u, we have ff * f = g 0, and then by [Carlen et al. 2020, Theorem 4], the second part of (1-25) follows. Note that if
(5-2)
for some α > 0, then the only choice of m that is consistent with (1-25) is m = d + 1.
We now specialize to d = 3, and impose the additional assumption on the potential. Recall that the Fourier transform of u (4-22) satisfies (4-25),
where S was defined in (4-24),
(5-4)
We split û into
with
so that, taking the large |k| limit in (4-25),
5A. Decay of U 1 . We first show that
(5-8) decays exponentially in |x|. We have
with
(5-10) Therefore, by (1-26),
(5-11) so, setting b := min(B, 1), (5)(6)(7)(8)(9)(10)(11)(12) and since A 4π
e -b(|x-y|+|y|) |x -y| dy = Ae -b|x| 4b 2 (b|x| + 1), (5-13)
we have
(5-14)
In addition,
(5-38) so which implies (5-19) in this case.
5B2. We have thus proved that S is analytic in H τ , which implies that the singularities of U 2 in H τ all come from the branch points of
(5-41)
Therefore, since F is analytic in a strip around the real axis, there exists an open set containing the real axis in which F has one and only one root, at 0. Thus the only branch point of √ F on the real axis is 0. Thus, U 2 is analytic in H τ .
5C. Decay of U 2 . We deform the integral to the path {iηy : 0 < y < |x| -τ } ∪ {iη|x| -τ + y : y > 0}
(5-42) and find with
(5-44)
(5-45)
5C1. We first estimate I 1 . We expand S: with β > 0 (since S is analytic and symmetric, and |S(|k|)| 1). Therefore, y → U 2 (i y) is C 2 for y = 0, and
(5-47) (5-50)
).
(5-52)
5C2. We now bound I 2 . Recall that, for κ ∈ ,ޒ we have |S(κ)| 1. ,
(5-53) provided |x| τ > 2C. Therefore, for large κ,
(5-54) so
for some constant C ′ > 0.
5C3. Inserting (5-52) and (5-55) into (5-43) and (5-15), we find that
(5-56) which, using (5-14), concludes the proof of the theorem.
6. Comparison with the Bose gas 6A. Sketch of the derivation of the simple equation. The simple equation (1-1)-(1-2) was originally derived [Lieb 1963] to approximate the ground state energy E 0 of a repulsive Bose gas, which is a system of N quantum particles interacting via the repulsive potential V. The ground state energy of this system is the lowest eigenvalue of the Hamiltonian operator
acting on the space of L 2 functions on the torus ޔ V of volume V. The corresponding eigenfunction, which we will denote by ψ N , satisfies
As is well known, by a Perron-Frobenius argument, ψ N is unique, nonnegative, and hence symmetric under exchanges x i ↔ x j and under translations.
We can write E 0 by integrating both sides of (6-2),
). The computation of E 0 thus reduces to that of g (2) N . Note that the kinetic energy does not appear explicitly in (6-3).
To compute g (2) N , integrate both sides of (6-2) with respect to x 3 , . . . , x N . This yields an equation relating g (2) N , g (3) N and g (4) N . The main approximation made in [Lieb 1963] is to write g (3) N and g (4) N as products of g (2) N factors: roughly,
This is a sensible approximation in the case of low-density ρ = N /V ≪ 1. Indeed, in this regime, one might expect ψ N to be approximately a Bijl-Dingle-Jastrow function,
for some appropriately chosen real function φ. Thus, ψ N is approximated by the partition function of a classical statistical mechanical model of particles interacting via the pair-potential φ. In this setting, g
N is the p-point correlation function of the canonical Gibbs distribution of this model. When (6-5) holds asymptotically as the particles move away from each other (remember, the density is low), the statistical mechanics system is said to satisfy the clustering property. There is a long literature on proving the clustering property for a large class of potentials φ; see, among many others, [Ruelle 1969;Gallavotti 1999;Pulvirenti and Tsagkarogiannis 2012].
Assuming the clustering property for the potential φ, the assumption (6-5) does not seem far fetched. This product structure leads to an equation for g (2) N . At this stage, one takes the thermodynamic limit: N → ∞ and ρ = N /V fixed. There are some subtleties to taking this limit, which are explained in [Lieb 1963]. Defining u := 1g (2) ∞ , the equation for u is [Lieb 1963, (3.29)]. After a few extra reasonable approximations, this equation reduces to (1-1). The equation for the energy (1-2) is simply the N → ∞ limit of (6-3). Plot of e/(4πρ) as a function of ρ on a log scale. The potential is V(r ) = e -r , in which case the scattering length is a ≈ 1.25. The solid curve is the energy computed from the simple equation (1-1)-(1-2), and the discrete points are the values of the energy of the Bose gas computed by M. Holzmann using a Monte Carlo algorithm. The gray area corresponds to a 5% error on the value of the energy. At low densities, we recover the Lenz asymptote e/(4πρ) ∼ a/2 and at high densities, we recover e/(4πρ) ∼ 1. The difference between the Monte Carlo simulation and the solution of the simple equation is smaller than 5%.
In particular, u is related to the correlation function g (2) of the Bose gas. The condition (1-6) that u(x) 1 is necessary to ensure that g (2) (x) 0. However, u(x) 0 is not a physical requirement, as g (2) (x) could, in principle, be > 1 for some x.
6B. Numerical comparison. One of the motivations for studying the simple equation is that it provides a simple tool to approximate the ground state energy of the Bose gas. In [Lieb and Liniger 1964], it was found that in one dimension the simple equation gives a value for the energy that differs from the Bose gas ground state energy by at most 69% (a more complete form of the equation yields an even better result with a maximal error of 19%). In one dimension, the difference is larger at high density.
In three dimensions, by Theorem 1.4, the simple equation predicts the correct low-density asymptote as the Bose gas. This is a not so surprising, since the derivation of the simple equation from the ground state equation of the Bose gas sketched above seems somewhat sensible when the density is low. However, when the density is high, at least in the case in which the potential has a nonnegative Fourier transform, the simple equation also yields the same asymptote as the Bose gas. In fact, considering the case
(6-7)
(which has a positive Fourier transform), we compared the ground state energy of the simple equation with values from a Monte Carlo simulation of the Bose gas computed by M. Holzmann (work in preparation), to whom we are most grateful for sharing his unpublished work. The comparison is in Figure 1, in which we found that the maximal error made by the simple equation, over the entire range of densities, is 5%! This is a promising result, which we will investigate in more depth and with more rigor in a later publication. . The asymptotic values are a ≈ 1.25 for ρ → 0 and 2 for ρ → ∞. This second derivative seems to be clearly positive, so ρe appears to be convex.
7A. Monotonicity. An important open problem is to show that e → ρ(e) is an increasing function. If the solution of the simple equation is in any way related to the ground state wave function of the Bose gas, then this should hold: if the density increases, the energy should increase. In addition, it would enable us to prove the uniqueness of the solution of the simple equation with fixed ρ, and might even allow us to generalize our result to potentials with hard core components, as well as to relax the constraint that V decays exponentially in Theorem 1.5. By running a few numerical computations, it seems clear that ρ(e) should be increasing; see Figure 1. Using a modified iteration in which ρ is fixed, we have proved that eρ(e) is strictly monotone increasing in e, but the proof that ρ(e) is as well has eluded us thus far.
Convexity. Another open problem is to prove that ρe(ρ) is a convex function, or, equivalently, that 1/ρ(e) is convex. In a physical setting, one expects ρe(ρ) to be convex. Indeed if ρe =: e v were not convex, there would exist ρ 1 < ρ < ρ 2 such that (ρ 1 + ρ 2 )/2 = ρ and e v (ρ 1 ) + e v (ρ 2 ) < 2e v (ρ). Furthermore, e v is the energy per unit volume, and, considering a volume V that is split into two equal halves, we find that a configuration in which one half of the volume holds a density ρ 1 of particles and the other holds ρ 2 would have energy V 2 (e v (ρ 1 ) + e v (ρ 2 )) < V e v (ρ).
(7-1)
Therefore, it would pay to have more particles in one half than in the other, which is unstable. Numerically, it seems quite clear that ρe(ρ) is convex; see Note that e appears only as the integral of S; see (1-2). While little is known rigorously about this equation, we have been studying it numerically in collaboration with M. Holzmann, and have found it to be remarkably accurate. These results will be detailed in a future publication.