Gross-Pitaevskii Limit of a Homogeneous Bose Gas at Positive Temperature

We consider a dilute, homogeneous Bose gas at positive temperature. The system is investigated in the Gross-Pitaevskii limit, where the scattering length a is so small that the interaction energy is of the same order of magnitude as the spectral gap of the Laplacian, and for temperatures that are comparable to the critical temperature of the ideal gas. We show that the difference between the specific free energy of the interacting system and the one of the ideal gas is to leading order given by 4πa ( 2̺2 − ̺0 ) . Here ̺ denotes the density of the system and ̺0 is the expected condensate density of the ideal gas. Additionally, we show that the one-particle density matrix of any approximate minimizer of the Gibbs free energy functional is to leading order given by the one of the ideal gas. This in particular proves Bose-Einstein condensation with critical temperature given by the one of the ideal gas to leading order. One key ingredient of our proof is a novel use of the Gibbs variational principle that goes hand in hand with the c-number substitution.


Background and summary
The experimental realization of the first Bose-Einstein condensate (BEC) in an alkali gas in 1995 [1,6] triggered numerous mathematical investigations on the properties of dilute Bose gases. The starting point was a work by Lieb and Yngvason [26] who proved a lower bound for the ground state energy of a dilute Bose gas in the thermodynamic limit. Together with the upper bound given in [24], it rigorously establishes its leading order behavior. In the case of hard-core bosons, the correct upper bound had already been proven in 1957 by Dyson [8]. Also the next-to-leading order correction to the ground state energy predicted by Lee, Huang and Yang in 1957 [16] could recently be proven, see [36] for the upper bound and [11] for the lower bound.
Bose gases in experiments are usually prepared in a trapping potential and such a set-up is well-described by the Gross-Pitaevskii (GP) limit. As shown in [24,20,21,27], the ground state energy of a Bose gas in this limit is to leading order given by the minimum of the GP energy functional. Additionally, a convex combination of projections onto the minimizers of this functional approximates the one-particle density matrix of the gas to leading order. Also in the GP limit the next to leading order correction to the ground state energy predicted by Bogoliubov in 1947 could be justified [5]. The accuracy reached in this work allows for an approximate computation of the ground state wave function and for a characterization of the low lying excitation spectrum. The dynamics of a system in the GP limit, on the other hand, can be described by the time-dependent GP equation, which was established in [9,10,3,28]. For a more extensive list of references we refer to [23,29,4].
While ground states provide a good description of quantum gases at very low temperatures, positive temperature effects are crucial for a complete understanding of modern experiments. In this case one is interested in the free energy and the Gibbs state of the system rather than in its ground state energy and in the ground state wave function. For the dilute Bose gas in the thermodynamic limit, the leading order behavior of its free energy per unit volume has been established, see [37] for the upper bound and [32] for the lower bound. The techniques developed in [26,24] have also been extended to treat fermions, both for the ground state energy [22] and for the free energy at positive temperature [30]. We mention also the papers [17,18,19,12] and [13] where Gibbs states of Bose gases with mean-field interactions are studied.
In a recent work [7], the trapped Bose gas at positive temperature is studied in a combination of thermodynamic limit in the trap and GP limit. It was shown that the difference between the free energy of the interacting system and the one of the ideal gas is to leading order given by the minimum of the GP energy functional. Additionally, the one-particle density matrix of any approximate minimizer of the Gibbs free energy functional is to leading order given by the one of the ideal Bose gas, but with the condensate wave function replaced by the minimizer of the GP functional. This in particular proves the existence of a BEC phase transition in the system. The proof of these statements relies heavily on the fact that particles in the thermal cloud have a much larger energy per particle, and therefore live on a much larger length scale than particles in the condensate. As a consequence, the interaction can be seen to leading order only in the condensate. The case of the homogeneous gas in a box, where the condensate and the thermal cloud live on the same length scale, was left as an open problem.
In the present work we consider this case, that is, we consider a homogeneous Bose gas (a gas in a box) at positive temperature in the GP limit. In this system the condensate and the thermal cloud necessarily live on the same length scale and interactions between them are relevant. We prove similar statements as in the case of the trapped gas in [7], in particular, we show the existence of a BEC phase transition with critical temperature given by the one of the ideal gas to leading order.

Notation
For functions a and b of the particle number and other parameters of the system, we use the notation a b to say that there exists a constant C > 0 independent of the parameters such that a ≤ Cb. If a b and b a we write a ∼ b, and a ≃ b means that a and b are equal to leading order in the limit considered.

The model
We consider a system of N bosons confined to a three-dimensional flat torus Λ of side length L (we could set L = 1 but we prefer to keep a length scale to explicitly display units in formulas). The one-particle Hilbert space is thus H = L 2 (Λ, dx), with dx denoting Lebesgue measure, and the Hilbert space of the N-particle system is the N-fold symmetric tensor product H N = L 2 sym (Λ N , dx). That is, H N is the space of square integrable functions of N variables that are invariant under exchange of any pair of variables. On H N we define the Hamiltonian of the system by 1 (1.1) Here ∆ denotes the Laplacian on the torus and d(x, y) is the distance between two points x, y ∈ Λ. The interaction potential is of the form The scattering length is a combined measure for the range and the strength of a potential and its definition is recalled in [23,Appendix C]. We are interested in the choice a v ∼ 1, i.e. a N /L ∼ N −1 . By definition, v is allowed to take the value +∞ on a set of positive measure which corresponds to hard core interactions. We will assume that v vanishes outside the ball with radius R 0 , that is, it is of finite range.
In the concrete realization of Λ as the set [0, L] 3 ⊂ R 3 , ∆ is the usual Laplacian with periodic boundary conditions and the distance function d(x, y) is given by d(x, y) = min k∈Z 3 |x − y − kL|. We also note that v N (d(x, y)) = k∈Z 3 The canonical free energy related to the Hamiltonian H N at inverse temperature β is defined by The unique minimizer of the Gibbs free energy functional is given by the canonical Gibbs state (1.7) Apart from the particle number N and the scattering length a v of the unscaled potential v, our system depends on the density ̺ = N/|Λ| and on the inverse temperature β. We are interested in the free energy of the system as N tends to infinity and for temperatures that are comparable to or smaller than the critical temperature of the ideal Bose gas, or equivalently, such that β β c . Here β c = (4π) −1 (ζ(3/2)/̺) 2/3 denotes the inverse critical temperature of the ideal Bose gas in a box of side length L and ζ denotes the Riemann zeta function. Since a N ∼ L/N this limit can be interpreted as a combined thermodynamic and GP limit, see also Remark 11 in Section 1.5 below.
For a given state Γ N we denote by γ N its one-particle density matrix (1-pdm) which we define via its integral kernel γ N (x, y) = Tr a * y a x Γ N . (1.8) In the above equation a * x and a x denote the usual creation and annihilation operators (actually operator-valued distributions) of a particle at point x, fulfilling the canonical commutation relations a x , a * y = δ(x − y). Here δ denotes the Dirac delta distribution. Equivalently, the integral kernel of γ N can be defined via the integral kernel Γ N (x 1 , ..., x N ; y 1 , ..., y N ) of the state Γ N by integrating out all but one coordinates and multiplying the result with N: Γ N (x, q 1 , ..., q N−1 ; y, q 1 , ..., q N−1 ) d(q 1 , ..., q N−1 ). (1.9) By definition, a sequence of states Γ N shows BEC if the related 1-pdms γ N have at least one eigenvalue of order N, i.e., (1.10)

The ideal Bose gas on the torus
In this section we recall some basic facts and formulas concerning the ideal Bose gas on a flat torus Λ. Since no explicit formulas are available for the canonical ensemble we state all results for the grand canonical ensemble. This is justified because the discussion in the Appendix shows that the free energy and the expected number of particles in the condensate, when computed in the two ensembles, agree with a precision that is sufficient for our purposes.
The expected number of particles in the condensate and in the thermal cloud (all particles outside the condensate) are given by . The ideal Bose gas shows a BEC phase transition in the limit of large particle number. More precisely, for N → ∞ one has and ζ the Riemann zeta function. Here [x] + = max{x, 0} denotes the positive part. For inverse temperatures such that β > β c (1 + ǫ) with ǫ > 0, the chemical potential is to leading order given by µ 0 ≃ −(βN gc 0 ) −1 ∼ −(βN) −1 , while for β < β c (1 − ǫ) it scales as µ 0 ∼ −β −1 . Finally, the free energy of the system is given by (1.13)

The main theorem
Our main result is the following statement: is a nonnegative and measurable function with compact support. Denote by ̺ 0 (β, N, L) the expected condensate density in the canonical Gibbs state of the ideal Bose. In the combined limit N → ∞, β̺ 2/3 ∼ 1 and a N given by (1.3) with a v > 0 fixed, we have for some α > 0. Moreover, for any sequence of states Γ N ∈ S N with 1-pdms γ N and Here γ N,0 denotes the 1-pdm of the canonical Gibbs state of the ideal gas in Λ and · 1 is the trace norm.
The fact that the difference between the specific free energy of the interacting system and the one of the ideal gas is given by 4πa(2̺ 2 − ̺ 2 0 (β, N, L)) also holds for the dilute Bose gas in the thermodynamic limit, see [37,32]. The formulas look the same because, as already mentioned above, our limit can be interpreted as a combined thermodynamic and dilute limit. We emphasize that (1.16) holds for any approximate minimizer of the Gibbs free energy functional in the sense of (1.15), and not only for the interacting Gibbs state (1.7) of the system. This in particular proves BEC for this class of states (see also Remark 9 below).
(For η → 0 the constants in the error terms in (1.14) and (1.16) blow up.) Our rate in Eq. (1.14) is the same as the one for the lower bound for the dilute Bose gas in the thermodynamic limit [32]. The known result for the ground state is implied by our result in the limit β → ∞. The error term is worse, however.
2. The result in (1.16) does not assume translation invariance of the states Γ N . If one assumes that Γ N is translation invariant the rate of convergence can be improved. In particular, one finds that the error term is bounded in terms of δ(N) 1/2 instead of on δ(N) 1/8 in this case.
3. Our result is uniform in the unscaled scattering length a v as long as a v ∈ (0, d] with 0 < d < ∞.
4. For β ∼ β c and a N ∼ L/N, we have F 0 ∼ |Λ|̺ 5/3 = L −2 N 5/3 for the free energy of the ideal gas, whereas the interaction energy is given by |Λ|a N ̺ 2 ∼ L −2 N. Up to this scale we control the free energy of the interacting gas.
5. The interaction energy is for β ≤ β c given by 8π|Λ|a N ̺ 2 to leading order, which has to be compared to 4π|Λ|a N ̺ 2 , its value at zero temperature. The additional factor of two is an exchange effect due to the symmetrization of the wave function, which only plays a role if the particles occupy two different one-particle orbitals. Above the critical temperature this is essentially always the case but particles inside a condensate do not experience this effect. This leads to the dependence of the interaction energy in Eq. (1.14) on ̺ 0 (β, N, L).
6. The free energy F 0 (β, N, L), the condensate density ̺ 0 (β, N, L) and the 1-pdm γ N,0 are the ones of the ideal gas in the canonical ensemble for which no explicit formulas are available. Our results are still valid if these three quantities are replaced by their corresponding grand canonical versions, as can be seen from the discussion in the Appendix.
7. Our bounds depend, apart from the scattering length a v , only on the range R 0 of the interaction potential. This dependence could be displayed explicitly. By cutting v in a suitable way one can extend the result to infinite range potentials which are integrable outside some ball with finite radius, that is, to all nonnegative potentials with a finite scattering length.
8. Our proof allows for the incorporation of internal degrees of freedom such as spin. For simplicity we only treat the case of spinless particles here.
9. Eq. (1.16) implies BEC into the constant function |Λ| −1/2 on the torus with condensate fraction given by the one of the ideal Bose gas to leading order. The statement follows from the fact that the trace norm · 1 bounds the operator norm · , and hence (1. 16) implies The critical temperature does not depend on the interaction in the dilute GP limit considered here. Deviations from this limiting value can be observed in experiments [34], however.
10. In the initial experiments BECs could only be prepared in harmonic traps. More recent set-ups also allow for the preparation of such systems in a box type potential with approximate hard wall boundary conditions, see [14]. The inclusion of these boundary conditions into our setting will be discussed in the next section.
11. Let us compare our setting to the one of the dilute Bose gas in the thermodynamic limit, which was considered in [32,37]. In the latter case one first takes the thermodynamic limit N, L → ∞ with ̺ = N/L 3 and β fixed, and afterwards considers the dilute limit where a̺ 1/3 ≪ 1 and β̺ 2/3 1. In our case we take the limit N → ∞ with a N ∼ L/N and β̺ 2/3 1. That is, we take a combined thermodynamic and dilute limit, where a N ̺ 1/3 ∼ N −2/3 . In this limit the spectral gap of the Laplacian is of the same order of magnitude as the interaction energy per particle, ∼ L −2 . This allows for a proof of BEC based on the coercivity of the relevant (free) energy functional, see (1.16), [20] and [7]. In contrast, proving BEC for the dilute Bose gas in the thermodynamic limit has been a major open problem in mathematical physics for almost a century. The spectral gap of the Laplacian closes in the thermodynamic limit and the system is expected to have Goldstone modes (sound waves in the case of the dilute Bose gas), that is, excitations with arbitrarily small energy. Accordingly, the relevant coercivity of the (free) energy is lost. We refer to [2] for an overview of an ambitious longterm project aimed at proving BEC in the thermodynamic limit with renormalization group techniques. Although it has so far not been possible to prove BEC for the dilute Bose gas in the thermodynamic limit, it is possible with our methods to obtain information on the 1-pdm when it is projected to suitably chosen high momentum modes. This should be compared to the case of the dilute Fermi gas in the thermodynamic limit, where comparable bounds have been proven in [30].

Extension to the case of Dirichlet boundary conditions
The methods developed for the proof of Theorem 1.1 also allow for the proof of a similar statement when the periodic boundary conditions are replaced by Dirichlet boundary conditions. In this case the system does not condense into the constant function but into the minimizer of the GP energy functional. To state this result, we first need to introduce some notation.
, we introduce the GP energy functional Here H 1 0 ([0, L] 3 ) denotes the usual Sobolev space of functions with zero boundary conditions. We denote by its ground state energy and by φ GP N,a its minimizer, which is unique up to a phase. One readily checks the scaling relations E GP (N, a, L) = L −2 NE GP (1, Na/L, 1) and φ GP N,a = N 1/2 φ GP 1,Na . Additionally, let H D N be the Hamiltonian (1.1), where the Laplacian on Λ is replaced by ∆ D , the Dirichlet Laplacian on [0, L] 3 . Similarly let F D (β, N, L) be the canonical free energy for H D N and let F D 0 (β, N, L) be the same quantity in the case v = 0. By ̺ D (x) we denote the density of the canonical Gibbs state of the ideal gas, and ϕ 0 is the ground state of −∆ D . Finally, we introduce with N D 0 the expected number of particles in the condensate of the Gibbs state of the ideal gas and denote by ̺ D th (x) = ̺ D (x) − N D 0 |ϕ 0 (x)| 2 the density of its thermal cloud. For simplicity we suppress the dependence of the densities on β, N and L.
The analogue of Theorem 1.1 in the case of Dirichlet boundary conditions is the following statement: is a nonnegative and measurable function with finite scattering length a v . In the combined limit N → ∞, β̺ 2/3 ∼ 1 and a N given by (1.3) with a v > 0 fixed, we have Moreover, for any sequence of states Γ N ∈ S N with 1-pdms γ N and 22) where γ N,0 denotes the 1-pdm of the canonical Gibbs state of the ideal gas and lim is our combined limit. Finally, lim 1 where · denotes the operator norm. In particular, we have BEC with the same condensate fraction and the same critical temperature as in the case of the ideal Bose gas to leading order.

Remarks:
1. In the case of periodic boundary conditions, the condensate wave function is given by a constant and therefore minimizes the GP energy functional on the torus, that is, (1.18) for functions φ ∈ H 1 per (Λ) (the usual Sobolev space of functions with periodic boundary conditions). For Dirichlet boundary conditions this picture changes because they force the minimizer of the GP functional to vary on the length scale of the box, that is, on L. This results in a macroscopic change of the energy of the condensate compared to the case with periodic boundary conditions. Although the Dirichlet boundary conditions do also change the free energy of the ideal gas compared to the case of periodic boundary conditions (not to leading order but on the scale we are interested in), they do not affect the density of the thermal cloud to leading order. This is because the energy per particle inside the thermal cloud is for β ∼ β c given by ̺ 2/3 , where ̺ = N/L 3 . Its density therefore varies on a length scale of order ̺ −1/3 which is much smaller than the length scale of the box: ̺ −1/3 /L ∼ N −1/3 . Hence, the density of the thermal cloud is essentially a constant until close to the boundary. Since the expected number of particles in the condensate does not depend on the boundary conditions to leading order this, in particular, implies that the second term in the bracket on the right-hand side of (1.20) satisfies The term on the right-hand side depends on the expected condensate density of the Gibbs state of the ideal gas in the case of periodic boundary conditions ̺ 0 and on ̺ th = ̺ − ̺ 0 . This should be compared to (1.14), where the same terms appear.
2. In the remaining part of the paper we will prove Theorem 1.1 but we will not prove Theorem 1.2.
The methods developed to prove Theorem 1.1 can, however, be adjusted to also obtain a proof for Theorem 1.2. Let us mention the main points to consider. Concerning the lower bound, the main point is that the technique from [32], that we use for the proof of the lower bound, naturally translates to the case of Dirichlet boundary conditions. This is because the c-number substitution is done with a sufficient number of momentum modes such that the GP minimizer, which varies on the length scale of the box, can be efficiently approximated with them. Additionally, as explained in the previous remark, the density of the thermal cloud of the ideal gas is constant to leading order. This allows one to use essentially the same technique to compute the free energy related to the modes that are not affected by the c-number substitution as in the case of periodic boundary conditions. To extend the proof of the upper bound, one has to cut the Fock space into high and low momentum modes, as it has been done in the proof of the lower bound. In the Fock space related to the low momentum modes one chooses the trial state to be a product wave function with N 0 (β, N, L) particles sitting in an approximate version of the GP minimizer. As above, N 0 (β, N, L) denotes the expected number of particles in the condensate of the ideal Bose gas. The overall trial state is then given by the symmetric tensor product of this function and a non-interacting canonical Gibbs state acting on the Fock space related to the high momentum modes (at the correct temperature and with N − N 0 (β, N, L) particles). In order to obtain the leading order behavior of the interaction energy, which depends on the scattering length, one has to, as in the case of periodic boundary conditions, add a correlation structure. The proof of the asymptotics of the 1-pdm remains up to minor adjustments unchanged. Here the main point is that the Griffith argument has to be done with an approximate version of the GP minimizer, which depends only on the low momentum modes of the c-number substitution, instead of with the constant function. Since the concrete implementation of the above strategy would considerably increase the length of the proof compared to the case of periodic boundary conditions, without adding substantial new difficulties, we only give the proof of Theorem 1.1 here.
3. We expect that the error bounds one obtains by following the strategy indicated by Remark 2 to prove Theorem 1.2 are not worse than those appearing in Theorem (1.1). In particular, we expect the same uniformity of the remainder in the inverse temperature as long as β β c .
4. Apart from periodic and Dirichlet boundary conditions we could also treat Neumann boundary conditions. Since the condensate function is a constant in this case one obtains the same statement as for periodic boundary conditions.

The proof strategy
Before we come to the proof of Theorem 1.1, we first briefly present the main steps to guide the reader.
Sec. 2 contains a proof of the upper bound for the free energy of the interacting gas. It is based on the Gibbs variational principle and the construction of a trial state whose free energy can be bounded from above by the desired expression. As a trial state we use the canonical Gibbs state of the ideal Bose gas. In order to obtain the scattering length in the interaction energy, we have to add a correlation structure which decreases the probability of finding two particles close together. Our ansatz yields a much simpler and shorter proof of the upper bound than the related proof in case of the dilute Bose gas in the thermodynamic limit [37]. This is only possible, however, because the scattering length scales as a N ∼ L/N, and hence the system is extremely dilute.
For the proof of the lower bound for F(β, N, L) in Sec. 3, we adjust the techniques developed for the related proof for the dilute Bose gas in the thermodynamic limit [32]. One key ingredient of our approach is a novel use of the Gibbs variational principle that goes hand in hand with the c-number substitution, which is a central ingredient of the proof in [32]. In comparison to [32], this allows us to work with a general state Γ instead of with a version of the grand canonical Gibbs state. In particular, we can keep the information that Γ has exactly N particles. This adjustment is essential for the proof of the asymptotics of the 1-pdm of approximate minimizers of the Gibbs free energy functional that we give in Sec. 4. Also in view of Sec. 4, we have to prove the lower bound for a slightly generalized Hamiltonian in which the energy of the lowest eigenfunction of the Laplacian is shifted by λ ≥ 0. We remark that, if one is only interested in the lower bound for the free energy, the technique from [32] can be applied essentially without modifications. More precisely, one would only need to consider all those terms in [32] that do not grow proportionally to the volume in the thermodynamic limit, and check that they are also of subleading order in the GP limit considered here (which is true).
The proof of the asymptotics of the 1-pdm γ N of an approximate minimizer of the Gibbs free energy functional is based on the novel use of the Gibbs variational principle mentioned above and has two main ingredients. The first ingredient is an estimate showing that γ N is, when projected to high momentum modes, given by the 1-pdm of the ideal gas to leading order. This part of the proof is motivated by a similar proof in [7] and is based on certain lower bounds for the bosonic relative entropy (the difference between two free energies) quantifying its coercivity. One main novelty in this part of our proof is a new lower bound for the bosonic relative entropy that allows us to simplify this part substantially w.r.t. the related part in [7]. In particular, it allows one to obtain better rates for the trace norm convergence of the relevant 1-pdm for given bounds on the relative entropy. In order to show the same statement for γ N projected to the low momentum modes, which is the second main ingredient of our proof, we apply a Griffith argument. Such arguments are based on the fact that differentiation of the free energy w.r.t. a parameter in the Hamiltonian yields the quantity one is interested in. In our case the parameter is the shift of the lowest eigenvalue of the Laplacian and the quantity of interest is the expected number of particles in the constant function, that is, in the condensate.

The variational ansatz
As trial state for the upper bound we choose the canonical Gibbs state of the ideal Bose gas on the torus and add a correlation structure. This is motivated by the following three observations: Firstly, the condensate wave function of the ideal gas on Λ is given by |Λ| −1/2 . If we turn on a repulsive interaction this is not going to change. Secondly, the free energy of the ideal gas is for β ∼ β c much larger than the interaction energy given the second term on the right-hand side of (1.14). This tells us that an approach based on first order perturbation theory should lead to the correct interaction energy. Finally, since v may contain a hard core repulsion and because the scaled pair interaction v N becomes very singular for large N, we need to assure that the probability of finding two particle close together is reduced compared to the ideal gas. This is achieved with the correlation structure in the spirit of [15]. In particular, it allows us to obtain the correct leading order of the interaction energy which is proportional to the scattering length. The idea to use a correlation structure in order to obtain the dependence of the energy of a dilute Bose gas on the scattering length has for the first time been used in [8] in the homogeneous case and in [24] in the inhomogeneous case.
The canonical Gibbs state of the ideal Bose gas on the torus Λ is given by As correlation structure we choose the Jastrow-like function [15,8,10] where b > 0 is a parameter to be determined and f 0 (|x|) is the unique solution of the zero-energy scattering equation see also [23,Appendix C]. We expand the canonical Gibbs state as where the functions Ψ α are chosen as symmetrized products of real eigenfunctions of the Laplacian on Λ. The final trial state with correlation structure is defined by and its free energy equals The remainder of this section is devoted to finding an appropriate upper bound for F ( Γ G N,0 ). We start with the computation of the energy.

The energy
We use the definition of our trial state (2.4) and write Bearing in mind that all eigenfunctions Ψ α of H N,0 are chosen to be real-valued, we integrate by parts once to rewrite the kinetic energy of the i-th particle as where dX is short for d(x 1 , . . . , x N ). For the energy of a single function Ψ α this implies where E α denotes its energy w.r.t. H 0 N , and the whole energy can be written as (2.9) The following Lemma provides a lower bound for the norm of FΨ α , and thereby an upper bound on the second term on the right-hand side of Eq. (2.9), as long as (4π/3)|Λ|̺ 2 a N b 2 < 1.
Proof. Spelled out in more detail, the norm of FΨ α reads where ̺ (2) Ψ α (x, y) denotes the two-particle density of Ψ α . Next, we use the fact that the Ψ α are symmetrized products of one-particle orbitals to conclude that holds. Here ̺ Ψ α is the one-particle density of Ψ α . Since the density of each Ψ α is a constant, we have ̺ Ψ α = ̺. This allows us to bound the integral on the right-hand side of Eq. (2.12) in the following way: (2.14) To obtain the bound for the integral of η b , we used its explicit form and the lower bound f 0 (|x|) ≥ [1 − a N /|x|] + , see [23,Appendix C]. In combination with (2.12), this proves the claim.
Next we analyze the numerator of the second term on the right-hand side of Eq. (2.9). We compute The square of this expression is given by These terms need to be inserted into the numerator of the second term on the right-hand side of Eq. (2.9) and we start with the first term on the right-hand side of the above equation. Introducing the function The off-diagonal terms in Eq. (2.16) can be bounded from above in a similar way by Combining these two bounds with (2.9) and (2.10), we obtain as an upper bound for the energy, where To derive this bound we had to assume that (4π/3)|Λ|̺ 2 a N b 2 < 1.
Let us denote by a * p and a p the usual creation and annihilation operators of a plane wave state ϕ p (x) = |Λ| −1/2 e ipx in Λ for p ∈ 2π L Z 3 . Also let n p = a * p a p be the related occupation number operator. To bound the first term on the right-hand side of (2.20), we use that Γ G N,0 has a fixed number of particles, and hence p∈ 2π L Z 3 n p can always be replaced by N when acting on Γ G N,0 . This implies When we use that all eigenfunctions ϕ p of −∆ are in absolute value equal to |Λ| −1/2 (they are plane waves), we come to the second line in the following inequality: By S 2 we denote the group of permutations of two elements. To arrive at the third line, we simply estimated n p (n p − 1) ≤ n 2 p and to come to the last line we used (2.21) and The second term on the right-hand side of Eq. (2.20) can be treated with a rough bound that we derive now. An application of the Cauchy-Schwarz inequality tells us that where S 3 denotes the group of permutations of three elements. We insert this bound into the second term on the right-hand side of Eq. (2.20) and obtain For the derivation of this result, we assumed (4π/3)|Λ|̺ 2 a N b 2 < 1 and a N < bη with 0 < η < 1. In the next step we will estimate the entropy of the state Γ G N,0 in terms of the entropy of Γ G N,0 and compute the final upper bound.

The entropy and the final upper bound
To relate the entropy of the state Γ G N,0 to the one of Γ G N,0 , we use [30, Lemma 2] which we spell out here for the sake of completeness.
Having the bound for the entropy at hand, we compute the free energy. With Eqs. (2.5), (2.26), (2.27) and (2.30) we find To obtain the result we assumed 4π 3 |Λ|̺ 2 a N b 2 < 1 and a N < bη with 0 < η < 1. Optimization yields b = (a N /(Na N ̺ + β −1 )) 1/3 and the bound We note that this bound is uniform in the parameter space β̺ 2/3 1. When we use a N LN −1 and β̺ 2/3 1, This completes the proof of the upper bound.

Proof of the lower bound
The proof of the lower bound proceeds along similar lines as the proof of the lower bound for the free energy in the thermodynamic limit [32]. One crucial ingredient in this work is a c-number substitution for momentum modes smaller than some cutoff which allows one to include a condensate. From a technical point of view, the proof in [32] is written in terms of the interacting Gibbs state of the system and uses the Berezin-Lieb inequality. The main difference between our setting and the one in [32] is that we also want to make a statement about the 1-pdm of approximate minimizers of the Gibbs free energy functional. To that end, we develop an alternative approach that is based on the Gibbs variational principle and goes hand in hand with the c-number substitution, and therefore also with the approach in [32]. To prove the statement about the 1-pdm of approximate minimizers of the Gibbs free energy functional, it will be necessary to prove the lower bound for the free energy related to the more general Hamiltonian where Φ(x) = |Λ| −1/2 and the index i indicates that the projection acts on the i-th particle. By adding this term we shift the energy of the lowest eigenvalue of −∆ by λ. In the following, we will assume that λ ∈ [0, (2π/L) 2 η] for some 0 < η < 1.
Before presenting the details of the lower bound, we give for the convenience of the reader a short summary of the main ideas of the proof in [32] in the context of the present setting.

Strategy of the proof of the lower bound
A key ingredient in the proof of the lower bound for the free energy of the interacting system in [32] is the observation that the interaction energy in (1.14) is, as long as β ∼ β c , much smaller than the free energy of the ideal gas F 0 (β, N, L) (compare with Remark 4 in Section 1.5). A naive version of first order perturbation theory is, however, not applicable because the interaction energy of the Gibbs state of the ideal gas is too large (it is even infinite if hard spheres are considered), see also the discussion in the beginning of Section 2. In the case of the GP limit with Dirichlet boundary conditions there is also a second obstacle, namely that the condensate wave function of the interacting system is not given by the ground state of the Laplacian in the box (as in case of the ideal gas), but rather by the minimizer of the GP energy functional (1.18).
The first problem is solved in [32] with the help of a Dyson Lemma. It allows to replace the singular and short ranged interaction potential v N by a softer potential with longer range. The price to pay for this replacement is a certain amount of kinetic energy. In the positive temperature setting it is important that only modes with momenta much larger than β −1/2 are used in the Dyson Lemma, since the other modes are needed to obtain the free energy of the ideal gas in (1.14).
After the replacement of the interaction potential a rigorous version of first order perturbation theory is applied. It is based on a correlation inequality [31] that is applicable to fermionic systems and to bosonic system above the critical temperature for BEC of the ideal gas. It allows to replace a general state in the expectation of the interaction energy by the Gibbs state of the related ideal gas and to estimate the error. An essential ingredient for this method is that the reference state in the perturbative analysis shows an approximate tensor product structure w.r.t. localization in different regions in position space. For a quasi-free state this is true if the off-diagonal of its 1-pdm decays sufficiently fast in position space. In order to overcome this shortcoming, coherent states are used in [32] to replace creation and annihilation operators of certain low momentum modes by complex numbers. In particular, this allows for the description of a condensate. Coherent states show an exact tensor product structure w.r.t. spatial localization in different regions in space, and therefore fit seamlessly into the above framework. In the case of Dirichlet boundary conditions, this approach also allows us to take into account that the condensate wave function is given by the GP minimizer, which solves the second problem from above.
The statement in Theorem 1.1 is uniform in the temperature as long as β β c . If the temperature is sufficiently low the free energy of the ideal gas in (1.14) is much smaller than the interaction energy and the approach from above cannot be expected to work. To extend the proof to this regime, we apply a different technique that uses in an essential way the zero temperature result in [26].

Reduction to integrable potential
In the proof of the lower bound we will make use of Fock spaces and, in particular, it will be required that the interaction potential has finite Fourier coefficients. As in [32, Section 2.1] we are therefore going to replace the potential v by an integrable potential. This is achieved with the following lemma whose proof can be found in [32, Sec. Since 0 ≤ v(r) ≤ v(r) for all r, we can replace v by v in the Hamiltonian H λ N for a lower bound. The Hamiltonian we obtain by this procedure will be denoted by H λ N . We also define a N to be the scattering length of the scaled potential v N .

Fock space
In the proof of the lower bound it is convenient to give up the restriction on the number of particles and to work in Fock space instead of in H N . In this section we introduce the necessary notation for this analysis. By µ(λ) we denote the chemical potential of the ideal Bose gas related to the one-particle Hamiltonian −∆ + λ|Φ Φ|, leading to an expected number of N particles, and we define µ 0 = µ(0). Let F be the Fock space over L 2 (Λ). We define the Hamiltonian H λ on F by where the kinetic and the potential energy operators are given by where ϕ has been introduced in the previous section. In the following we will denote the grand canonical kinetic energy operator for λ = 0 by T and similarly for the full Hamiltonian.

Coherent states and the Gibbs variational principle
In this section we introduce a formalism that allows us to apply a c-number substitution while still keeping information on a given state Γ whose free energy we want to investigate. We start by introducing notation for the c-number substitution. Let us pick some p c > 0 and decompose the Fock space as F F < ⊗ F > , where F < and F > denote the Fock spaces of the momentum modes with |p| < p c and |p| ≥ p c , respectively. The trace over F < will be denoted by Tr < and similarly for F > . To keep the notation simple and because we do not expect it to cause confusion, we will denote the traces over F and H N by the same symbol Tr. By M we denote the number of momenta p ∈ 2π L Z 3 with |p| < p c . For a vector z ∈ C M we introduce the coherent state |z ∈ F < by where |vac denotes the vacuum in F < . Coherent states of this kind form an overcomplete basis with C M |z z| dz = 1 F < . Here dz = M |p|<p c dz p with z p = x p + iy p and dz p = dx p dy p π . For every state Γ on the Fock space F , we define the operator Γ z acting on F > by Since Γ is a state, ζ Γ (z) dz is a probability measure on C M . The entropy of the classical distribution ζ Γ is defined by On the level of the Hamiltonian, we will need the lower symbol of H λ which is defined by H λ s ≔ z, H λ z . It is an operator-valued function from C M into the unbounded operators on F > . Since a p |z = z p |z , the lower symbol can be obtain from H λ by simply replacing a p by z p and a * p by z p for all |p| < p c . By H λ,s (z) we denote the upper symbol of the Hamiltonian H λ which is defined by the identity (3.9) To compute it, one has to replace |z p | 2 by |z p | 2 − 1 in the lower symbol and similarly with other polynomials in z p , see [25].
The following Lemma shows that the entropy of a state Γ can be bounded from above in terms of the expectation of the entropies of the states acting on F > w.r.t. the probability measure ζ Γ (z) dz, plus one additional term that quantifies the entropy of the classical distribution ζ Γ (z).
Lemma 3.2. Let Γ be a state on F . The entropy of Γ is bounded in the following way: Proof. We write the first term on the right-hand side of (3.11) as To prove the result, we need to show that holds. To that end, we expand We apply Jensen's inequality to show that Ψ z . Hence, the left-hand side of (3.14) is bounded from below by The measure λ α Ψ z α 2 dz is a probability measure with respect to summation over α ∈ N and integration over C M in z. Another application of Jensen's inequality therefore tells us that To come to the last line, we used ∞ α=1 |Ψ z α Ψ z α | = 1 F > . This proves the claim.
With the definitions from above and Lemma 3.2, we can derive a lower bound for the Gibbs free energy functional. Let Γ be a state on H N ⊂ F . Eq. (3.9) allows us to write the expectation of the energy as In combination with Lemma 3.2, this implies Although the upper symbol naturally appears in the above inequality, it is more convenient to work with the lower symbol instead. Let N s = |z| 2 + |p|≥p c a * p a p denote the lower symbol of the particle number operator. The difference between the upper and the lower symbol ∆H λ (z) = H λ s (z) − H λ,s (z) can be written as as well as We will later choose the parameters p c and ϕ such that Z (1) ≪ |Λ|a N ̺ 2 . Eq. (3.23) is the formula we were looking for. It should be compared to (2.3.9) and (2.3.10) in [32], in which a version of the grand canonical Gibbs state of the interacting system appears. In contrast to that, (3.23) allows to use the c-number substitution while still working with a given state Γ. The Gibbs variational principle applied to Tr H λ Υ z − 1 β S (Υ z ) will later allow us to obtain information on an approximate minimizer of the Gibbs free energy functional (1.5), see Sec. 4.
Remark 3.1. In [32] the additional term is added to the second quantized Hamiltonian before relaxing the restriction on the particle number. Like this one obtains a strong control on the expected number of particles in the system. We do not need this term in our approach because the information that the state Γ has exactly N particles is still encoded in the Fock space formalism through the state Γ z and the measure ζ Γ (z) dz.
In the remaining part of this section we will go through the proof in [32], mention changes due to our approach and collect the necessary results. The following sections will be named like the ones in [32].

Relative entropy and a-priori bounds
In this section we derive an a-priori bound for states Γ whose free energy is small in an appropriate sense. This bound is the only information we are going to need about the state to prove the lower bound.
For two general states Γ and Γ ′ on Fock space we denote by the relative entropy of Γ with respect to Γ ′ . It is a nonnegative functional that equals zero if and only if Γ = Γ ′ . Let Γ 0 be the Gibbs state corresponding to T s (z) at inverse temperature β on F > , which is independent of z.
We emphasize that T s (z) is the lower symbol of the grand canonical kinetic energy operator with λ = 0. Since the interaction potential v N and λ are nonnegative, we have Let us integrate both sides of the above equation with ζ Γ (z) dz over C M . The first term on the right-hand side equals The chemical potential µ 0 is negative because the lowest eigenvalue of −∆ equals zero. From the Gibbs variational principle we know that To arrive at the last line we used the inequality x ≥ 1 − e −x . Eqs. (3.26)-(3.28) therefore imply Note that we have chosen µ 0 such that the expected number of particles in the grand canonical system equals N. Concerning the lower bound, it is sufficient to consider states Γ with free energy bounded from above by Here F 0 (β, N, L, λ) denotes the canonical free energy for the Hamiltonian H λ N with v = 0. The actual lower bound we are going to prove will be smaller than the right-hand side of (3.30), that is, the statement will hold independently of this assumption. We use Lemma A.1 in the Appendix to obtain an upper bound for the canonical free energy in Eq. (3.30) in terms of the grand canonical free energy, that is, (1.13) with µ 0 replaced by µ(λ) in the first term and p 2 − µ 0 replaced by p 2 + δ p,0 λ − µ(λ) in the second term. Together with (3.23) and (3.30), this implies To obtain an upper bound on the difference between the two grand canonical potentials in the first line, we write The second estimate follows from µ 0 ≤ µ(λ) which is implied by the monotonicity of the map λ → µ(λ). In combination with (3.31) and This is the a-priori bound we were looking for. To compute the interaction energy, we will use (3.34) to replace Γ z by Γ 0 in a controlled way. In other words, the lower bound represents a rigorous version of first order perturbation theory.
Remark 3.2. The interacting free energy corresponding to H λ N depends on λ only through the free energy of the ideal gas to leading order. This is because the interaction energy depends, apart from |Λ|, a N and ̺ (which are independent of λ), only on the expected density of the condensate ̺ 0 (β, N, L, λ). It can be checked, however, that ̺ 0 (β, N, L, λ) does not depend on λ to leading order if λ ∈ [0, (2π/L) 2 η] with 0 < η < 1. This justifies the use of Γ 0 in the computation of the interaction energy.
We also derive a second a-priori bound. It is a simple estimate for the variance of the probability measure ζ Γ (z) dz which counts the number of particles in the Fock space with momenta smaller than or equal to p c and reads

Replacing vacuum
In order to prove the lower bound, we have to estimate the kinetic energy and the interaction energy of states of the form Υ z = U(z)|vac vac|U(z) * ⊗ Γ z with U(z) defined in (3.5), and where Γ z obeys the a-priori bound (3.34). We find it necessary for this analysis to replace the vacuum in the formula for Υ z by a more general quasi-free state, which we do in a controlled way in this section. This will become important below when the interaction energy of Υ z is computed. For this purpose the latter will be replaced by a quasi-free state, whose one-particle density matrix should show rapid off-diagonal decay in order for the localization technique of the relative entropy to be applicable, see [32, Sections 2.8, 2.13]. Hence the momentum distribution needs to be sufficiently smooth and cannot vanish identically for the low momentum modes (as it does in the vacuum state).
We denote by Π a particle-number conserving quasi-free state on F < . It is fully determined by its 1-pdm π = |p|<p c π p |p p|. Here |p denotes a plane wave state in L 2 (Λ) with momentum p. We also define P = |p|<p c π p = tr π. Here and in the following we denote by tr[·] the trace over the one-particle Hilbert space H. Finally, let us introduce the state Υ z π = U(z)ΠU(z) * ⊗ Γ z (3.37) on F . In order to replace Υ z by Υ z π in a controlled way, we have to estimate the effect of this replacement on the kinetic and the potential energy. Our analysis follows the one in [32, Section 2.5] with the only difference that we control the particle number with the measure ζ Γ (z) dz and not with the help of the operator K, see Remark 3.1. More concretely, we use the identity C M Tr > [N s Γ z ]ζ Γ (z) dz = N + M. When we go through the analysis in [32, Section 2.5], we obtain We also have to replace Υ z by Υ z π in the kinetic energy which can be done with the identity In combination with (3.23), we obtain as a lower bound for the free energy of Γ.

Dyson Lemma and Filling the Holes
The sections 2.6 (Dyson Lemma) and 2.7 (Filling the Holes) in [32] remain basically unchanged. To introduce several quantities that are needed later and to mention the necessary changes due to the term λa * 0 a 0 in the Hamiltonian, we collect the main result here. The Dyson Lemma [32, Lemma 2] is used to replace the singular and short ranged potential v N by a softer potential with a longer range at the expense of a certain amount of kinetic energy. To be precise, only the high momentum modes are used for the Dyson Lemma. This is necessary because the low momentum modes are used to obtain the free energy of the ideal Bose gas. The Dyson Lemma naturally leads to an effective interaction potential with a hole around zero. Because it will be necessary for the computation of the interaction energy, this potential is replaced by a slightly different one without a hole.
By R we denote the length scale of the effective potential from the Dyson Lemma satisfying 10R 0 L/N < R < L/2. When this potential is replaced by a potential without a hole in the middle, one obtains a potential with a slightly reduced scattering length with two parameters 0 < κ < 1 and ǫ > 0 that are related to the Dyson Lemma, see [32,Sec. 2.6]. The definition of the function j is given by We apply the Dyson Lemma and the analysis to replace the relevant potential by one without a hole in the same way as in [32, Sec. 2.7] and compute the term λa * 0 a 0 separately. The final result of this analysis reads In the above equation T c s (z) denotes the lower symbol of the operator We will later choose κ and R such that a N (R 0 L/N) 2 /R 3 ≪ κ. This in particular implies κ ′ > 0. The function ν : R 3 → R + is chosen such that ν(p) = 0 for |p| ≤ 1, ν(p) = 1 for |p| ≥ 2, and 0 ≤ ν(p) ≤ 1 in-between. It is used to implement the fact that only the high momentum modes are used in the Dyson Lemma. The parameter s obeys s ≥ R and will later be chosen such that s ≫ R. We will also choose κ ≪ 1. In combination with λ ≤ (2π/L) 2 η with 0 < η < 1, this implies ǫ(p) > 0 for all p. The effective interaction potential W will not be specified here because we will use the same estimate for Tr WΥ z π as in [32]. Its definition can be found in [32,Sec. 2.7]. Note that, compared to [32, (2.7.15)], we have the additional term β −1 S (Γ z , Γ 0,λ c ) in our lower bound (3.43). Here Γ 0,λ c denotes the grand canonical Gibbs state for the kinetic energy operator T c s (z) which is independent of z and depends on λ only through the chemical potential µ(λ). The additional term is not important for the lower bound (it is positive and could be dropped), but it will be important for the proof of the asymptotics of the 1-pdm of approximate minimizers of the Gibbs free energy functional in Sec. 4. When we insert the above result into (3.40) and argue as in (3.27) and (3.28), we find From the first two terms on the right-hand side of (3.46), we will obtain the free energy of the ideal gas.

Localization of Relative Entropy
In this section we introduce notation that will be important for the following. The main result from the related section in [32] will not be stated since we will not explicitly need it. It is used only in parts of the proof in [32] that we do not have to adjust.

Interaction Energy Part 1 -3
The expectation of the effective potential W in the state Υ z π is estimated as in [32,]. The result of the analysis in these sections is the following lower bound: The scattering length a ′ N has been defined in (3.41), m(r) is an explicitly given smooth function that vanishes faster than any power for |r| → ∞, compare with [32,Sec. 2.10], and c > 0. Additionally, where j is defined in (3.41). To obtain the result, we started with Eqs. [32, (2.11.19-21)] and the same choice of the parameters ǫ and D as in [32]. We also applied Jensen's inequality to the term proportional to the relative entropy and used (3.35) We assumed that p c ̺ 1/3 , ̺ ω ̺ and R s. The bound is valid for any choice of the parameter 0 < b ≤ L/2 that has been introduced in the previous section. We will later choose b such that bp c ≫ 1 and βb −2 ≪ 1.

A bound on the number of particles
The lower bound on the interaction energy contains a term of the form Recall that we will later choose R ≪ s, that is, the term Tr N Υ z π − Ω z b is multiplied by a positive constant. In this section we will first rewrite the integral over the trace on the right-hand side of (3.53). This way it will be apparent that the term in (3.53) can be combined with another error term that we will find in Sec. 3.11. This term will be of the same form but it will be multiplied by a negative constant that is much smaller than the one in the equation above. Accordingly, we only have to derive a lower bound for the integral on the right-hand side of (3.53) to finally estimate the sum of these two terms.
Let us start by rewriting the integral on the right-hand side of (3.53). We note that Tr[NΩ z b ] = |z| 2 + Tr[NΩ b ] = |z| 2 + Tr[NΩ π ] as well as Tr[NΥ z π ] = |z| 2 + Tr[NΥ π ] and we denote by N > = |p|≥p c a * p a p the particle number operator on the Fock space F > . Since Υ π = Π ⊗ Γ z and Ω π = Π ⊗ Γ 0 we can write We also know that In combination, Eqs. (3.54) and (3.55) imply This is the first result we were looking for.
Next, we will derive a lower bound for the right-hand side of (3.54). It implies a lower bound on the righthand side of (3.56) that will later allow us to estimate the relevant error term in Sec. 3.11. A bound of this kind was proved in [32,Sec. 2.12] in the case of the dilute Bose gas in the thermodynamic limit. In this limit momentum space sums can be replaced by integrals because the relevant errors do not grow proportionally to the volume and are therefore irrelevant. In the GP limit that we consider these error terms have to be quantified, however. To that end, we have to adjust the estimates in Eqs. (2.12.9)-(2.12.12) in [32], which will be done with the help of the following Lemma. holds.
Proof. Assume first that κ − √ 3 2π L ≤ 0. In this case we drop the characteristic function on the left-hand side of (3.57) for an upper bound. Next, we write the sum over p as the sum over those p that are an element of one of the coordinate planes, i.e. with one coordinate p i equal to zero, plus a sum over all remaining p. To estimate the sum over all remaining p, we interpret the sum as a lower Riemann sum and find that it is bounded from above by L 2π 3 R 3 f (|p|) dp. (3.58) The sum over those p that are an element of the coordinate planes can be estimated similarly. Here we write the whole sum as a sum over those p that are an element of one of the coordinate axes plus the sum over all remaining p. The sum over the remaining p is estimated again by interpreting it as a lower Riemann sum. For one such coordinate plane, we find Because there are three coordinate planes, we have three such terms. It remains to estimate the sum over those p that are an element of one of the coordinate axes of R 3 . Again by interpreting the sums over the three coordinate axis as lower Riemann sums, we find that they are bounded from above by In order to write the two-dimensional integral from (3.59) in terms of a three-dimensional integral, we use A similar computation can be done for the term in Eq. (3.60). Putting these estimates together proves (3.57) in this case. The bound in the case κ − √ 3 2π L > 0 can be obtained similarly. Here we only have to realize that κ − √ 3 2π L is the radius of the largest ball such that the integral over its complement is an upper bound to the relevant three-dimensional lower Riemann sum. This proves the claim.
To adjust the analysis in [32] after (2.12.8), we have to find an upper bound for the sum An application of Lemma 3.3 tells us that it is bounded from above by A short computation shows that the expression in the above equation cannot be larger than This bound replaces (2.12.12) in [32]. Using the above and (2.12.8) in [32], we conclude that holds.
If we combine (3.65) with (3.54) and (3.56) we obtain the bound An application of Jensen's inequality and the a-priori bound (3.34) therefore imply This is the bound we were looking for. It will later be used to bound the relevant error term in (3.51).

Relative Entropy, Effect of Cutoff
In this section we estimate the relative entropy S (Υ z π , Ω z b ) = S (Π ⊗ Γ z , Ω b ), which appears in the lower bound (3.51) for the interaction energy, in terms of S (Π ⊗ Γ z , Ω π ) = S (Γ z , Γ 0 ). Since we have an a-priori bound for the integral w.r.t. ζ Γ (z) dz over the latter expression at hand this will allow us to finalize the lower bound for the interaction energy. Compared to [32], we have to adjust how the momentum space sum related to [32, (2.13.21)] is estimated.
We are faced with estimating the term Here with some parameters D, t and q that are specified in [32] and that are chosen such that β(p 2 − µ 0 ) − Dβtq 2 > 0 for all p ∈ 2π L Z 3 . When we insert the estimate [32, (2.13.24)] for ω t into (3.68), we see that it is bounded from above by a constant times (3.71) We will later choose p c and b such that bp c ≫ 1, and hence τ > 0. Next, we bound the summands in (3.70) from above by a monotone function with the same behavior at zero and at infinity. Afterwards we use Lemma 3.3 to see that (3.70) is bounded from above by a constant times This is the estimate for the term in (3.68) we intended to show. The remaining part of the analysis in [32] can be done similarly. With the a-priori bound (3.34) and the estimate a N ≤ a N , we finally arrive at To obtain the result, we used that βb −2 is small enough and that bp c is large enough.

Final lower bound
We have obtained all necessary estimates to complete the lower bound for the free energy of Γ. To that end, we collect the estimates from the previous sections, that is, (3.46), (3.51), (3.56) and (3.73) and find To obtain this result, we used the definition (3.41) of a ′ N and γ b ≤ ̺ ω ≤ ̺. The first part of the inequality for ̺ ω follows from the definition of γ b (3.52) and the second part from the choice of π p in [32] after (2.13.15). Using the definition of π p again, we estimated as in [32,Sec. 2.14]. To replace a N by a N in the term in the second line in (3.74) we applied Lemma 3.1 with the choice ǫ = a N N/(Lϕ). We will later choose ϕL/N ≫ a N . The error terms Z (1) and Z (2) are defined in (3.21) and (3.38), respectively.
To obtain a bound for the interaction term in (3.74), we write The last term in the first line of (3.77) can be dropped for a lower bound. Next we combine the first term in the third line of (3.74) with 4πa N |Λ| times the term in the second line of (3.78) integrated with ζ Γ (z) dz over C M . Together, they read We will later choose R ≪ ̺ −1/3 , that is, the term in the second bracket in (3.79) is negative and we need a lower bound for Such a bound is provided by (3.67). In combination, the results of this paragraph imply (3.64). To obtain the result, we also used M |Λ|p 3 c .
In the following, we assume p c 0. The case where p c = 0 will follow easily from the analysis of this case. Using the definition for ̺ gc 0 (β, N, L) in Sec. 1.4, we see that ̺ − ̺ 0 ≥ ̺ gc 0 which implies that we obtain a lower bound for the terms in the second line in (3.81) when we replace ̺ − ̺ 0 by ̺ gc 0 . In order to derive a lower bound for ̺ ω , we estimate To obtain the last bound, we used Lemma 3.3. The integral in the second line is not larger than a constant times p c /β and the sum is bounded by a constant times (βL) −1 . Hence, When we follow the argumentation in [32, Sec. 2.14] and invoke Lemma 3.3, we see that It remains to replace the grand canonical condensate density ̺ gc 0 (β, N, L) by its canonical version ̺ 0 (β, N, L). This can be achieved with the help of Lemma A.3 in the Appendix which implies We recall that A has been defined in (3.64). The result has been obtained under the assumption that R̺ 1/3 is small enough.

The non-interacting free energy
In this Section we derive a lower bound for the first two terms on the right-hand side of (3.74). The dispersion relation ǫ(p) has been defined in (3.44). The following Lemma will be necessary to derive a lower bound for the second term.
Lemma 3.4. The chemical potential µ(λ) satisfies Proof. The lower bound follows from the fact that the map λ → µ(λ) is monotone increasing. The lowest eigenvalue of the operator p 2 + δ p,0 λ is given by min (2π/L) 2 , λ . To prove the upper bound, we realize that 1 exp β min 2π The above statement follows from the fact that we have N particles in the system, and hence the expected number of particles in the condensate cannot exceed N. Eq. (3.90) is equivalent to the second inequality in (3.89) and proves the claim.
A long as |p| ≤ 1/s we have ǫ(p) = (1 − κ + κ ′ )p 2 + δ p,0 λ − µ(λ) and ǫ(p) ≥ κ ′ p 2 − µ(λ) holds for |p| > 1/s. We will later choose the parameters such that s 2 /(L 2 κ ′ ) is much smaller than one. Together with λ ≤ ( 2π L ) 2 η for 0 < η < 1, this in particular implies that ǫ(p) > 0 for |p| ≤ 1/s and κ ′ p 2 − µ(λ) > 0 for |p| > 1/s. In accordance with this decomposition of the momenta, we split the sum in the first line on the right-hand side of (3.74) into two parts. The first part is given by To arrive at the right-hand side, we used the concavity of the map x → ln(1 − e −x ). An application of Lemma 3.4 together with the assumption that |κ − κ ′ | is small enough tells us that the absolute value of the term in the second line is bounded from above by a constant times The summands times β are bounded from above by a constant times e −βp 2 /8 . An application of Lemma 3.3 therefore tells us that the sum in (3.92) cannot be larger than a constant times |Λ|/β 5/2 .
The part of the sum in the first line of (3.74) coming from the momenta with |p| > 1/s is bounded from below by As already mentioned in the discussion after Lemma 3.4, we will later choose κ ′ and s such that is a negative and monotone increasing functions, we can use Lemma 3.3 to show that the right-hand side of (3.93) is bounded from below by a constant times We will later choose κ ′ and s such that βκ ′ /s 2 ≫ 1. We also note that the term in (3.94) is an exponentially decaying function of this parameter. Putting the results of this section and the definition of κ ′ in (3.45) together, in particular, (3.91)-(3.94), we find To obtain the result, we also used s ≤ L and Lemma A.1 in the Appendix to replace the grand canonical free energy F gc 0 (β, µ(λ), L, λ) by the canonical free energy F 0 (β, N, L, λ). The bound has been derived under the assumptions that s 2 /(L 2 κ ′ ), s 2 /(βκ ′ ) and |κ − κ ′ | are small enough.

Choice of Parameters
Optimization under the assumptions a N = a v L/N with fixed a v > 0, λ ≤ ( 2π L ) 2 η with fixed 0 < η < 1 and β̺ 2/3 1 leads to the same choice of parameters as in [32,Sec. 2.16] and implies the lower bound and some function s → C δ (s) that is uniformly bounded for s ∈ [d, ∞) with d > 0. For δ → 0, the function C δ blows up. The error term is of lower order as long as N −2/3 (β̺ 2/3 ) 5/2 ≪ 1. The bound is uniform for 0 < a v 1. The desired uniformity in the temperature will be achieved in the next section.
Since it will be needed in Sec. 4, we also state here the choices of p c and R resulting from the optimization. They are given by

Uniformity in the temperature
We follow the analysis in [32,Sec. 2.17] until equation (2.17.4) and arrive at To obtain this bound, it has been assumed that κ = (a 3 N ̺) 1/17 and R = a N (a 3 N ̺) −5/17 . From this point on we have to adjust the analysis in [32]. This is necessary because we cannot replace sums by integrals, we have to add the term δ p,0 λ to the one-particle Hamiltonian, and we want to obtain the canonical free energy and the canonical condensate density (in the thermodynamic limit the canonical and the grand canonical free energies and condensate densities are the same).
In order to obtain the final estimate, we also need to replace the interaction energy in Eq. (3.100) by the formula we have in Theorem 1.1. As above we denote by ̺ 0 (β, N, L) the expected condensate density of the ideal gas in the canonical ensemble in the case λ = 0 and we define by ̺ th = ̺ − ̺ 0 the expected density of the thermal cloud. We then have for β β c To come to the second line, we used ̺ 0 ≤ ̺ as well as Lemma A.2 in the appendix to bound ̺ th ̺ gc th . To see that ̺ gc th (3.106) The first term on the right-hand side is bounded by a constant times L 2 /β. To bound the second term, we invoke Lemma 3.3 to see that it is bounded from above by a constant times |Λ|/β 3/2 .
In combination, (3.104) and (3.105) imply: In the derivation of this bound, we assumed that s 2 /(L 2 κ) and s 2 /(βκ) are small enough. Before we optimize under the assumption λ ≤ ( 2π L ) 2 η with fixed 0 < η < 1 and a N = a v L/N with a v > 0 fixed, we insert κ and R from above, see the discussion after (3.100). With some δ > 0 we choose and ǫ 2 = 1 This fulfills the condition on s and κ stated after (3.107) and it assures the smallness of the term in the third line in (3.107). It implies We have to combine the two bounds (3.96) and (3.109) in the same way as in [32,Sec. 217] to obtain the optimal rate, that is, we use (

Proof of the asymptotics of the one-particle density matrix
In this section we prove the claimed asymptotics for the 1-pdm of approximate minimizers of the Gibbs free energy functional. A crucial input for the analysis in this section are the lower bounds (3.96) and (3.109). The proof is split into four parts: In the first part we consider the 1-pdm projected onto the subspace of the oneparticle Hilbert space with momenta at least p c and we show that it equals the one of the non-interacting Gibbs state to leading order. In the second step we consider the 1-pdm projected to the orthogonal complement of that subspace and show that also there it is to leading order given by the one of the non-interacting Gibbs state.
In the third step we estimate the off-diagonal contributions and in the fourth part we prove the uniformity in the temperature. We highlight that off-diagonal contributions to the 1-pdm have to be estimated because we do not assume that the states under consideration are translation invariant. With this assumption we would obtain a better rate. An important example of a translation invariant state is the interacting Gibbs state (1.7).

The one-particle density matrix of the thermal cloud
Let Γ be an approximate minimizers of the Gibbs free energy functional in the sense that with 0 ≤ δ(N) = o(1) in the considered limit. Together with the lower bound (3.96) with λ = 0, this implies The state Γ 0 c ≡ Γ 0,λ=0 c was defined in Sec. 3.6. The index c refers to the fact that the relevant dispersion relation is not p 2 − µ 0 but the one we obtained after applying the Dyson Lemma, see (3.44). The goal of this section is to obtain quantitative information on the 1-pdm γ of Γ from this bound. Let us define P = 1 (−∆ < p c ) and Q = 1 − P. When projected to the high momentum modes, γ reads where γ z is the 1-pdm of the state Γ z . Hence, if we denote by γ 0 c the 1-pdm of Γ 0 c , we have In the following, we will derive a bound on the right-hand side.
The starting point of our analysis is (4.2). Since Γ 0 c is a quasi-free state, the left-hand side (4.2) can be bounded from below in terms of the bosonic relative entropy. For two nonnegative operators a, b with finite trace, it is defined by where σ(x) = x ln(x) − (1 + x) ln(1 + x) and {λ i , ψ i } and {η j , ϕ j } denote the eigenvalues and eigenfunctions of a and b, respectively. We also denote by S(a) = − tr [σ (a)] (4.7) the bosonic entropy of a. We then have S(γ z ) ≥ S (Γ z ), see [35, 2.5.14.5], as well as S (Γ 0 c ) = S(γ 0 c ), and therefore conclude that In combination with (4.2), this implies In order to obtain quantitative information from Eq. (4.9), we need the following Lemma which quantifies the coercivity of the bosonic relative entropy. It is an improved version of [7, Lemma 4.1].
Lemma 4.1. There exists a constant C > 0 such that for any two nonnegative trace-class operators a and b we have (4.10) Proof of Lemma 4.1. Let f (x, y) = σ(x) − σ(y) − σ ′ (y)(x − y). In the proof of Lemma 4.1 in [7] it has been shown that there is a C > 0 such that .
x + √ y 2 which allows us to bound the right-hand side from below by To obtain the final estimate, we assumed that y ∈ [0, y max ]. In combination with (4.6), this proves Next, we write the difference between the two density matrices as and estimate their trace norm difference by Here · 2 denotes the Hilbert-Schmidt norm. Together with (4.13) and this proves the claim.   Given the size of the remainder in (3.96), an improved rate of convergence is not of particular relevance for the analysis here, however. Bounds for the trace norm difference of two one-particle density matrices in terms of the relative entropy of the related states have recently also been proven in [19, Theorem 6.1].
Before we apply Lemma 4.1, we will show that tr γ z ≤ N holds. To that end, we write 16) where N > = |p|≥p c a * p a p , as before. Let us denote by P N the projection onto the N-particle sector of the Fock space F . It is sufficient to show that holds. With [N > , |z z|] = 0 = [N > , P N ] we check that [P N N > P N , P N |z z|P N ] = 0. But this implies the claim. With this information at hand, we apply Lemma 4.1 to the left-hand side of (4.9) and additionally use tr γ 0 c N as well as 1 + γ 0 c L 2 /β. We find that holds. It remains to replace γ 0 c by the canonical 1-pdm of the ideal Bose gas.
To that end, we first replace the dispersion relation ǫ(p) (3.44) in the definition of γ 0 c by p 2 − µ 0 . This can be done with an analysis that is very similar to the one carried out between (3.91) and (3.94) in Sec. 3.12 and yields (4. 19) In order to replace γ gc 0 by its canonical analogue γ 0 , we invoke Lemma A.3 in the Appendix to show that L 2 /β one checks that the term in the second line of (4.20) is bounded from above by a constant times N 5/6 ln(N) 1/2 uniformly in β̺ 2/3 1. In combination (4.18)- (4.20) imply The bound yields the desired result as long as N −2/3 (β̺ 2/3 ) 5/2 ≪ 1.

The one-particle density matrix of the condensate
In order to investigate PγP and, in particular, to show the existence of a BEC, we apply a Griffith argument. From Eq. (3.110), we know that where the Hamiltonian H λ N was defined in (3.1) and 0 ≤ λ ≤ ( 2π L ) 2 η with some fixed 0 < η < 1. Together with (4.1) this implies for some 0 ≤ λ ≤ λ. As above we denoted by Φ the constant function with value |Λ| −1/2 on the torus. The second derivative in the above equation is nothing but (−β) times the variance Var λ (n 0 ) = n 2 0 λ − n 0 2 λ of the occupation of the p = 0 orbital. Here and in the following, we denote by · λ the expectation in the canonical ensemble with the energy of the p = 0 orbital shifted by λ. We also recall that n p = a * p a p . In order to bound the above variance, we need the following Lemma:  Together with (4.25) and (4.26), this proves the claim.
We use Lemma 4.2 to bound the second term on the right-hand side of (4.23). The choice λ = ( 2π L ) 2 /2 implies together with a N N/L the bound (4.28) We highlight that the right-hand side of (4.28) is uniform in β̺ 2/3 .

The off-diagonal of the one-particle density matrix and the final estimate
In this section we are going to control the off-diagonal parts of γ which will allow us to give the final estimate. Our analysis follows the lines of a similar analysis in [7,Sec. 4.3]. We write γ − γ 0 1 ≤ P (γ − γ 0 ) P 1 + 2 P (γ − γ 0 ) Q 1 + Q (γ − γ 0 ) Q 1 (4.39) and estimate the right-hand side term by term. A bound for the first term on the right-hand side was given in (4.38). From (4.21) we know that the last term is bounded by N( c ℓ (β, N, L) + δ(N)) 1/2 .
The first term in the bracket on the right-hand side can be bounded with (4.38), the second with (4.31). We find  This proves the claimed asymptotics for the 1-pdm as long as N −2/3 (β̺ 2/3 ) 5/2 ≪ 1. In the next section we discuss the uniformity in β̺ 2/3 .

Uniformity in the temperature
In order to show the desired uniformity in the temperature, we have to consider the case where β̺ 2/3 is so large that c ℓ (β, N, L) is no longer small, that is, β̺ 2/3 N 4/15 . In this case we have ̺ 0 (β, N, L) ≃ ̺, and hence the contribution of the thermal cloud to the 1-pdm of the ideal gas is of lower order. In combination with (4.28), this will imply a similar statement for γ. In particular, it will allow us to conclude that γ − γ 0 1 ≪ N uniformly in β̺ 2/3 .

A. Some properties of the ideal Bose gas
In this appendix we collect three Lemmas concerning properties of the ideal Bose gas, which have been proven in [7, Appendix A] or follow from a statement there, and which we need in the main text. In particular, they concern the comparison of relevant quantities when computed in the canonical and in the grand canonical ensemble. Although these statements hold more generally, we state them here only for the ideal Bose gas on the torus Λ.
As in the main text we denote by F 0 (β, N, L, λ) the canonical free energy of the ideal gas related to the Hamiltonian (3.1) with v = 0 and by F gc 0 (β, µ, L, λ) its grand canonical analogue. We recall that λ ≥ 0 denotes the shift of the lowest eigenvalue of the Laplacian on Λ. Similarly, · N and · gc,µ denote the expectations and N 0 and N gc 0 the expected number of particles in the condensate in the two ensembles (for simplicity we have suppressed the λ-dependence here). The expected number of particles in the grand canonical ensemble is denoted by N(µ) and γ 0 /γ gc 0 is the 1-pdm of the canonical/grand canonical ideal gas (which depend on λ). The following three statements hold.
Lemma A.1. Assume µ is such that N = N ∈ N. Then