Bose-Einstein Condensation in a Dilute, Trapped Gas at Positive Temperature

We consider an interacting, dilute Bose gas trapped in a harmonic potential at a positive temperature. The system is analyzed in a combination of a thermodynamic and a Gross-Pitaevskii (GP) limit where the trap frequency $\omega$, the temperature $T$ and the particle number $N$ are related by $N \sim (T / \omega)^{3} \to\infty$ while the scattering length is so small that the interaction energy per particle around the center of the trap is of the same order of magnitude as the spectral gap in the trap. We prove that the difference between the canonical free energy of the interacting gas and the one of the noninteracting system can be obtained by minimizing the GP energy functional. We also prove Bose-Einstein condensation in the following sense: The one-particle density matrix of any approximate minimizer of the canonical free energy functional is to leading order given by that of the noninteracting gas but with the free condensate wavefunction replaced by the GP minimizer.


Background and summary
Proving Bose-Einstein condensation (BEC) rigorously for systems of interacting particles has for a long time been a major challenge in Mathematical Physics. The experimental realization of BEC in trapped alkali gases in 1995 [1,8] triggered numerous mathematical investigations of the properties of dilute Bose gases. Building on work of Dyson on hard-core bosons from 1957 [11], the first proof of an asymptotically accurate lower bound for the ground state energy of a dilute Bose gas in the thermodynamic limit was achieved in [30]. Together with the upper bound given in [29], it established firmly the leading order behavior of the ground state energy. Perhaps more important than the result itself were the techniques of the proof which formed the basis of much subsequent work.
For dilute, trapped gases as prepared in typical experiments the Gross-Pitaevskii (GP) limit for the ground state is relevant. This limit is characterized by the requirement that the interaction energy per particle is kept of the same order of magnitude as the spectral gap in the trap. Mathematically, this can be achieved by either scaling the trap or the interaction potential suitably as the particle number tends to infinity. In [29,25,26] it was proved that the ground state energy of a Bose gas is in this limit equal to the minimum of the GP energy functional. Additionally, the projection onto the minimizer of this functional is the limit of the one-particle density matrix of the gas, proving complete Bose-Einstein condensation in the ground state. The dynamics of a system in the GP limit, on the other hand, can be described by the time-dependent GP equation, see [12,13,2,34]. For a more extensive list of references to the mathematical analysis of dilute Bose gases we refer to [28,37,3].
While ground states provide a good description of quantum gases at very low temperatures in first approximation, the understanding of finite temperature effects in cold gases is also essential. For instance, the spectacular emergence of a peak in the momentum distribution out of a maxwellian thermal cloud in the experiments [1,8] as the temperature falls below a critical value, cannot be explained from the ground state properties alone. For describing such phenomena one has to consider the Gibbs state and the free energy of the system rather than the ground state and the corresponding energy. For a dilute, homogeneous Bose gas, the leading order behavior of the free energy has been established, see [40] for the lower bound and [46] for the upper bound. The techniques developed in [30,29] have also been extended to treat fermions, both for the ground state [27] and for the free energy at finite temperature [39]. We mention also the papers [20,21] and [14] where Gibbs states of Bose gases with mean-field interactions are studied. A general proof of Bose-Einstein condensation in dilute gases remains elusive, however.
In this paper we consider the Gibbs state of an interacting Bose gas in a harmonic trap at positive temperatures in a combination of a thermodynamic and a GP limit. We show that in this limit the free energy becomes equal to that of the ideal gas plus a correction given by the GP energy of the condensate. Moreover, we show that the one-particle density matrix of the system is well approximated by that of the ideal gas with the noninteracting condensate wavefunction replaced by the minimizer of the GP energy functional. This proves, in particular, Bose-Einstein condensation at positive temperatures with the same transition temperature and the same condensate fraction as for the ideal gas to leading order.

Notation
For functions a and b depending on the particle number or other parameters, 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 trapped in a three-dimensional harmonic potential with trap frequency ω. The one-particle Hilbert space is H 1 = L 2 (R 3 ) and that for the whole system is the N -fold symmetric tensor product H N = ⊗ N s L 2 (R 3 ). On H N we define the Hamiltonian of the system by 1 In this formula ∆ i denotes the Laplacian acting on the i-th particle and we have subtracted 3 2 ω for convenience so that the ground state energy of the harmonic oscillator is zero. The interaction potential is given by with a nonnegative, radial and measurable function v, independent of N . A simple scaling argument shows that if a v is the (dimensionless) scattering length of v then the scattering length a N of v N is The scattering length is a combined measure of the range and the strength of a potential and its definition is recalled in Section 1.9. To be able to include hard-core interactions, we allow v to take the value +∞.
If v happens to be infinite on a set of nonzero measure the domain of the Hamiltonian has to be restricted to functions that vanish on this set. We require that v is integrable outside some finite ball in order for the scattering length a v to be finite.
To motivate the scaling (1.2) we recall that the energy per particle of a dilute gas of density ρ with interaction potential v N is to first approximation proportional to ρa N [30,29]. If ρ ∼ N ℓ −3 osc , where ℓ osc ∼ ω −1/2 denotes the length scale of the trap, we see that i.e., the interaction energy per particle is of the order of the spectral gap ω as N → ∞. Note also that the dimensionless "gas parameter" ρa 3 N tends to zero as N −2 so the scaling (1.2) amounts to considering a special case of a dilute limit.

The Gross-Pitaevskii energy functional and the GP limit
The GP energy functional E GP with trapping potential as in Eq. (1.1) and scattering length a is defined by Its ground state energy is E GP (N, a, ω) = inf The unique minimizer of E GP will be denoted by φ GP N,a . To keep the notation simple, we suppress its dependence on ω. The energy and the minimizer satisfy the scaling relations For a detailed discussion of the mathematical properties of the GP functional and its relation to the ground state of the Hamiltonian (1.1) we refer to [25] (see also [26,28,29,32]) where the following is proved: In the GP limit, where a = a N as in Eq. (1.3) and N → ∞, the ground state energy per particle of the many-body Hamiltonian (1.1) converges to ωE GP (1, N ω 1/2 a N , 1) = ωE GP (1, a v , 1). Moreover, the normalized one-particle density matrix of the ground state wavefunction converges in trace norm to the projector onto φ GP 1,avω −1/2 .
Since the GP minimizer differs considerably from the gaussian ground state of the harmonic oscillator if a v is large, the interaction can leave a clear mark on the density profile of the gas despite the high dilution imposed by the GP limit, as seen in experiments [17,7]. In fact, for large a v ("Thomas-Fermi limit") the profile has approximately an inverse parabolic shape of extension ∼ a 2/5 v ℓ osc .
Our goal is to generalize these results to equilibrium states at positive temperatures when the GP limit is combined with the natural thermodynamic limit in the trap. The definition of the latter and the heuristics behind our main results can be deduced from a comparison of the length scales involved in the problem, as discussed next.

Length scales, thermodynamic limit, heuristics
For a noninteracting gas at inverse temperature β = T −1 the following length scales are relevant: • The extension ℓ osc ∼ ω −1/2 of the ground state of the harmonic oscillator.
• The extension of the thermal cloud in the trap, R th ∼ ω −1 β −1/2 , obtained by equating the potential energy ω 2 R 2 th and the thermal kinetic energy β −1 .
• The mean particle distance d th ∼ N −1/3 R th in the thermal cloud.
The thermodynamic limit is defined by keeping the ratio d th /ℓ th fixed as N → ∞, i.e. by the condition N (βω) 3 fixed. (1.8) The thermodynamic limit requires in particular (βω) ∼ N −1/3 → 0 and thus ω ≪ T . If d th ℓ th , i.e., N (βω) 3 1, thermal de Broglie wave packets overlap and condensation can be expected. This heuristics is confirmed by the standard analysis of the ideal Bose gas in the harmonic trap [35,33,5]: Bose-Einstein condensation takes place if the temperature T is smaller than the critical temperature T c given by Here ζ is the zeta-function and ζ(3) = 1.202 . . . To be more precise, denote by N 0 (β, N, ω) the expected number of particles occupying the ground state of the harmonic oscillator in the canonical ensemble. If N → ∞ with N (βω) 3 fixed the condensate fraction is given by with [t] + = max{t, 0}. The condition for the right-hand side of Eq. (1.10) to be larger than zero, i.e., T < T c , is equivalent to lim N (βω) 3 > 1.202. (1.11) In this case the ground state of the harmonic oscillator is macroscopically occupied in the ideal Bose gas. For T > T c on the other hand, i.e., if there is no condensation. These formulas are most conveniently derived in the grand canonical ensemble.
For an assessment of the effects of interactions the following observation is crucial: The length scales R th and ℓ osc become separated if N → ∞, The average density of the condensate, 14) and that of the thermal cloud, 15) with N th := N − N 0 , are therefore widely different in the condensation regime 0 < T < T c where N 0 and N th are comparable: The same holds for the ratios of the interaction energy per particle, a N ρ th and a N ρ 0 respectively. We remark that in the presence of a condensate described by a GP minimizer with large a v it would be more precise to replace ℓ osc in Eq. (1.14)-(1.15) by a 2/5 v ℓ osc . Also, since the ball |x| a 2/5 v ℓ osc is essentially excluded for the thermal cloud, R th in Eq. (1.15) should be replaced approximately by R th + a 2/5 v ℓ osc . As long as a v stays O(1), however, the asymptotic behavior in Eq. (1.16) remains valid.
The separation of scales expressed by Eq. (1.13) leads to the following expectations for the combined GP and thermodynamic limit: • The thermal cloud of the ideal gas remains essentially intact.
• BEC takes place for T < T c and the condensate can be described by the GP minimizer, residing close to the center of the trap.
Transforming this heuristic picture into a mathematical proof is the subject of this paper. In order to state our main results precisely we need a few more definitions.

Gibbs state, free energy and the concept of BEC
The canonical Gibbs state for the Hamiltonian (1.1) is with Z(β, N, ω) = Tr H N e −βH N the canonical partition function. The free energy of the system at inverse temperature β = T −1 is given by The trace in Eqs. (1.17)-(1.18) is taken over H N , that is, over the subspace of permutation symmetric functions in L 2 (R 3N ). In the following we will drop the index H N and just write Tr for this trace. By F 0 (β, N, ω) we denote the free energy of the ideal Bose gas in the harmonic trap, that is, the one for v = 0. and using this definition, the free energy can be written as (1.20) The unique minimizer of F N is given by the canonical Gibbs state defined in Eq. (1.17).
The reduced one-particle density matrix of a state Γ N ∈ S N is defined via the integral kernel Here a * x and a x denote creation and the annihilation operators (actually operator-valued distributions) of a particle at point x, fulfilling the canonical commutation relation [a x , a * y ] = δ(x−y). Alternatively, γ N can be defined as N times the partial trace of Γ N over N − 1 particle variables.
Finally, Bose-Einstein condensation for a sequence of states Γ N means, by definition, that . Let F (β, N, ω) be the corresponding free energy, F 0 (β, N, ω) the free energy of the ideal gas, and N 0 (β, N, ω) the expected number of particles occupying the ground state of the harmonic oscillator in the canonical Gibbs state of the ideal Bose gas. In the combined thermodynamic and GP limit, that is, for N → ∞, (βω) −3 ∼ N and a N as in Eq. (1.3) with a v fixed, we have Moreover, for any sequence of states Γ N ∈ S N with Here γ N,0 denotes the one-particle density matrix of the noninteracting canonical Gibbs state, ϕ 0 is the normalized ground state wavefunction of the harmonic oscillator Hamiltonian h = −∆ + 1 4 ω 2 x 2 − 3 2 ω, and · 1 stands for the trace norm. Finally, where · is the operator norm. In particular, Bose-Einstein condensation takes place with the same transition temperature T c and the same condensate fraction as for the ideal gas to leading order, with the GP minimizer macroscopically occupied while the occupation of every state orthogonal to it is o(N ).

Remarks
1. The bounds leading to Theorem 1.1 are uniform in (βω) −3 N −1 as long as this quantity remains in a compact interval [c, d], 0 ≤ c < d < ∞ of the real line. That is, we need not require lim T /T c to exist, only to stay bounded. In particular, Theorem 1.1 continuously extrapolates to the known result at T = 0.

2.
We have the following uniformity of our bounds in the scattering length: Assume v(x) = a −2 vṽ (x/a v ) for some potentialṽ with scattering length equal to one. By scaling, the scattering length of v is given by a v . Then the bounds in Theorem 1.1 are uniform for a v ∈ (0, d] with 0 < d < ∞. 3. The free energy F 0 (β, N, ω) of the ideal Bose gas in the harmonic trap is of order N β −1 ∼ ω(βω) −4 , see Section 1.10 below. The GP energy, on the other hand, is of order N ω ∼ ω(βω) −3 which is the scale up to which we have to control the free energy F (β, N, ω) of the interacting gas.
4. Theorem 1.1 is stated and proven for the explicit choice of 1 4 ω 2 x 2 as a trapping potential. This choice is mainly for notational simplicity, but is also motivated by the fact that the harmonic trap is the physically most relevant one. Rotational symmetry is not important, however. The treatment of an anisotropic trap requires only slight notational modifications and the interpretation of ω as the geometric mean of the principal frequencies of the parabolic potential.

5.
The impressive cover picture of the first Bose-Einstein condensates in the July 1995 issue of Science, where the paper [1] appeared, shows the momentum distribution rather than the spatial distribution of the trapped gas. The momentum distribution of the thermal cloud is approximately an isotropic maxwellian of width ∼ β −1/2 . The condensate momentum distribution, which is the modulus squared of the Fourier transform of the GP minimizer, is anisotropic because the trap was anisotropic. The width of the peak in momentum space is ∼ ω 1/2 and thus narrower than the thermal cloud by a factor (βω) 1/2 . Our Theorem confirms this picture rigorously for the first time.
6. The techniques used in the proof of Theorem 1.1 carry over with moderate adjustments to the case of trapping potentials behaving as |x| α with α < ∞ for large |x|. The key point is that such potentials still lead to an asymptotic power law behavior of the eigenvalues of the related Schrödinger operator and cause a separation of length scales between the condensate and the thermal cloud. It should be noted, however, that if α > 2 the exponent 1/6 in (1.13) is replaced by a smaller exponent and a clear separation of scales thus requires even larger values of N . 3 If α → ∞ the whole system is confined in a box and the condensate and the thermal cloud are no longer spatially separated. Treating such systems will require a different approach from the one of the present paper. This is an important open problem because traps with very large α have recently become available in experiments [15,31].
7. The fact that the transition temperature for BEC and the condensate fraction for the interacting gas stay the same as for the ideal gas relies essentially on the diluteness of the system, expressed through the scaling (1.2), and the N → ∞ limit. Under less restrictive conditions finite size corrections can be expected and have been seen in experiments [43]. Extending our results to capture these effects requires a proof of BEC beyond the GP limit, a difficult unsolved problem.
8. Since we work in the canonical ensemble, explicit expressions for the free energy F 0 (β, N, ω), the one-particle density matrix γ N,0 and the condensate fraction N 0 /N in the ideal Bose gas are not available. However, Theorem 1.1 remains valid if these expressions are replaced by their corresponding grand canonical versions, which we recall in Section 1.10. This is due to the fact that the difference between the two ensembles is negligible in the limit of consideration here (see Section 1.10 and the discussion in the Appendix for details).
9. Theorem 1.1 is also valid if we replace F (β, N, ω) by its grand canonical analogue F gc (β, N, ω) where a chemical potential µ is chosen such that the expected number of particles equals N . We know from the Gibbs variational principle that F gc (β, N, ω) ≤ F (β, N, ω) holds. Hence the upper bound for the canonical free energy directly implies the upper bound for its grand canonical version. On the other hand, the proofs of the lower bound for the free energy and of the asymptotics of the one-particle density matrix are carried out in a way that is directly applicable to the grand canonical ensemble.
10. The strategy of the proof of Theorem 1.1 also applies in two space dimensions with the obvious modifications, compare with [28, Chapter 6.2].

Supplementary 1: The scattering length
Let us quickly summarize the basic facts about the scattering length. A more detailed discussion can be found in [28,Appendix C]. Our assumptions on v guarantee that the zero energy scattering equation has a unique solution. It satisfies for some constant a > 0 which is called the scattering length of v. The scattering length has the natural interpretation of a combined measure for the range and the strength of the interaction potential v. If v happens to be a hard core potential with range r for example, one finds a = r, whereas for weak, integrable potentials a ≈ (8π) −1 R 3 v(x)dx. For dilute quantum gases where collisions can be described in a low energy approximation, the leading order contribution to the scattering amplitude comes from s-wave scattering. In this approximation particles are scattered in every direction with the same probability and the scattering cross-section is given by 4πa 2 .
In the case of nonnegative potentials, the scattering length can be characterized via the following variational principle. Denote by X the set of functions in H 1 loc (R 3 ) with φ(x) → 1 for |x| → ∞. Then the scattering length is given by By a trial state argument one can improve this inequality and show that it is strict if v is not identically zero. For typical quantities of interest, as for example the free energy, there do not exist simple closed form expressions in the canonical ensemble. Nevertheless, in the thermodynamic limit as defined in Section 1.5, these quantities are very close to those computed in the grand canonical ensemble, which allows for explicit computations. We start by introducing several grand canonical quantities, and subsequently discuss their relation to their canonical versions.
The ideal Bose gas in the harmonic trap is described by the one-particle Hamiltonian The expected number of particles in the condensate and the thermal cloud are given by respectively. Here g(n) = (n + 1)(n + 2)/2 is the degeneracy of the energy level ωn of h. For βω ≪ 1 the sum in the above equation can be interpreted as a Riemann sum and one finds N gc th ∼ (βω) −3 . The expected number of particles in the gas is N = N gc 0 + N gc th . It will be adjusted by the chemical potential µ 0 such that N = N holds. To be more precise, we choose µ 0 = µ 0 (β, N, ω) such that (1.10). On the other hand, for T > T c (1 + ǫ) one has −µ 0 ≃ ηT , where η is the unique solution of the equation The grand canonical free energy is given by In the same way as for the expected number of particles in the thermal cloud, the sum in the above equation can be interpreted as a Riemann sum and one finds F gc 0 (β, µ 0 , ω) ∼ −ω(βω) −4 , compare with [35,Eqs. (10.19)- (10.22)]. Note, however, that the error one makes by approximating the sum in Eq. (1.34) by an integral is of order ω(βω) −3 ∼ ωN which is the order of the GP energy. Hence, we cannot do this replacement and have to work with the sum. The same is true in the case of N gc th .
The canonical free energy F 0 is given by (1.18) with the last term much smaller than the main contribution in the thermodynamic limit. Moreover, Corollary A.2 tells us that |N 0 − N gc 0 | (βω) −3/2 (ln N ) 1/2 + (βω) −1 ln N . The same bound holds for the expected numbers of particles in the thermal cloud.

The proof strategy
The rigorous mathematical implementation of the heuristic picture behind Theorem 1.1 is technically rather involved. For the convenience of the reader we describe here briefly the main steps before turning to the proof in the remaining sections. Section 2 contains an upper bound on the free energy F (β, N, ω) that has the correct asymptotic form (1.23). It utilizes the Gibbs variational principle (1.20), so the task is to construct an appropriate trial state. The latter consists on the one hand of a pure state, describing the condensate particles, which are confined to a ball of some size R with ω −1/2 ≪ R ≪ ω −1 β −1/2 , i.e., large compared to the oscillator length but small compared to the length scale of the thermal cloud. On the other hand, the remaining particles, constituting the thermal cloud, will be described by a suitably modified Gibbs state confined to the complement of the ball. For the condensate, we can use the known zero-temperature results to obtain the GP energy. Even though the particle interactions affect the free energy to the order we are interested in only through the condensate, we cannot simply use the noninteracting Gibbs state as a trial state for the thermal cloud, since we want to allow for nonintegrable interaction potentials (e.g., having a hard-core) which can have infinite energy in that state. We thus have to appropriately modify the trial state, avoiding configurations where the particles are too close. This creates some technical complications which have to be dealt with carefully. For various bounds, it turns out to be necessary to compare certain expressions for the ideal Bose in the canonical and grand canonical ensembles, respectively. The relevant estimates are collected in Appendix A.
In Section 3 we shall give a lower bound on the free energy, which together with the upper bound proves Eq. (1.23). We shall use the technique of Fock-space localization to spatially divide the system into two, one confined to a ball of radius R (chosen as above to satisfy ω −1/2 ≪ R ≪ ω −1 β −1/2 ) and one confined to the complement. Inside the ball, the effect of the positive temperature is of lower order, and we can again utilize the known zero-temperature results to obtain a bound on the energy of these particles, as well as on the one-particle density matrix, which displays Bose-Einstein condensation into the GP minimizer. For the system in the complement of the ball, we can drop the interaction terms (using their positivity) to obtain the free energy of the ideal gas as a lower bound.
To obtain information on the one-particle density matrix of the interacting Gibbs state (or approximate Gibbs state) we develop in Section 4 a novel lower bound on the free energy functional for an ideal Bose gas quantifying its coercivity. More precisely, in Lemma 4.1 we show that any approximate minimizer of the Gibbs free energy functional is, in a suitable sense, close to the actual minimizer. In combination with the result on Bose-Einstein condensation for the system inside the ball of radius R, this allows us to prove Eqs. (1.25) and (1.26).
Finally, Appendix A collects certain properties of the ideal Bose gas that we need in our proofs.

The variational ansatz
In this section we construct a trial state Γ N whose free energy has the correct asymptotics (1.23). In the case of the ideal Bose gas in the harmonic trap, the characteristic length scale of the condensate is ω −1/2 while for the thermal cloud it is ω −1/2 (βω) −1/2 which is much larger in the limit we consider.
The main idea of the proof is based on the expectation that in the GP limit this picture does not change if an interaction is turned on. What also does not change to leading order is the expected number of particles in the condensate and in the thermal cloud. Since the thermal cloud is therefore much more dilute than the condensate, the free energy of the particles outside the condensate is not affected by the interaction to the same order of magnitude as the condensate. The following analysis makes this intuition precise.
The first step in the construction of our trial state Γ N is to decompose space into three disjoint parts, a ball B(R) with radius R > 0, an annulus A(R, R + ℓ) with radii R and R + ℓ and the complement of a ball with radius R + ℓ. All those sets are assumed to be centered around zero. For later convenience we will refer to them as Region 1, 2 and 3, respectively. The length R will be chosen such that ω −1/2 ≪ R ≪ ω −1/2 (βω) −1/2 , i.e., between the length scale of the condensate and the one of the thermal cloud. Choosing R like this, we will be able to spatially separate the system into two parts, a condensate living in B(R) and a thermal cloud living in B(R + ℓ) c without affecting each of them too much. In Region 2 there will be no particles in the trial state. The length ℓ ∼ ω −1/2 N −1 is chosen such that there is no hard core interaction between particles in Region 1 and particles in Region 3.
The one-particle Hilbert space and the Fock space naturally decompose as , respectively. The heuristics in Section 1.5 tells us that we can neglect the contribution from the thermal cloud in Region 1 since it is too dilute to contribute with a macroscopic number of particles. The condensate will be described by the ground state of the Hamiltonian acting on L 2 sym (B(R) N 0 ), the space of permutation symmetric square integrable functions depending on N 0 variables. Here, ∆ D i,≤R denotes the Laplacian on B(R) with Dirichlet boundary conditions, acting on the i-th particle, and N 0 = N 0 (β, N, ω) is the expected number of particles in the condensate of an ideal Bose gas in the canonical ensemble. Strictly speaking N 0 is not necessarily an integer and we should rather choose ⌈N 0 ⌉, the smallest integer larger than or equal to N 0 , in the definition of H D ≤R instead. In the end, this will lead to 1/N corrections which do not cause any additional difficulties, however. In order not to complicate the presentation unnecessarily, we therefore assume that N 0 is an integer. By Ψ D ≤R we denote the unique ground state wavefunction of H D ≤R with energy E D ≤R . The above construction will allow us to use existing results for the ground state of a Bose gas in a trap.
In Region 3 on the other hand, the condensate does not contribute to the free energy to leading order because its extension ω −1/2 is much smaller than R. To describe the thermal cloud, we define the one-particle Hamiltonian where ∆ D i,≥R+ℓ denotes the Laplacian on B(R + ℓ) c with Dirichlet boundary conditions. We also define the noninteracting N th -particle operator with energy cut-off Λ 3) acting on the i-th particle. By 1(h D ≥R+ℓ,i ≤ Λ) we denote the spectral projection onto the subspace of H 3 where h D ≥R+ℓ,i is at most Λ. The cut-off Λ in Eq. (2.4) is introduced for technical reasons, which will be explained in the text preceding Lemma 2.5 in Subsection 2.3 below. Let the many-particle projection P N th be defined by Its range consists of linear combinations of symmetrized products of eigenfunctions of h D ≥R+ℓ , where each of these one-particle functions has energy at most Λ. By we denote the canonical Gibbs state associated with the Hamiltonian H D,Λ ≥R+ℓ .
Since the interaction potential may include a hard core repulsion between the particles we have to add a correlation structure to the state Γ D,Λ ≥R+ℓ . For this purpose we define the Jastrow-type function [18] F where b is a parameter to be determined and f 0 (|x|) is the unique solution of the zero-energy scattering The parameter b will be chosen to be larger than the scattering length a N but of the same order of magnitude.
We expand Γ D,Λ ≥R+ℓ = ∞ α=1 λ α |Ψ α Ψ α | where we choose the functions Ψ α as symmetrized products of eigenfunctions of the one-particle Hamiltonian h D ≥R+ℓ with energy at most Λ. The energy of Ψ α is denoted by E α , that is, Here and in the following, Ψ denotes the L 2 -norm of Ψ.
After these preparations we can now finally define our trial state where Ω denotes the Fock space vacuum in F(H 2 ). Let be the noninteracting N th -particle Hamiltonian in the region B(R + ℓ) c without the cut-off. Because the first two factors in Eq. (2.9) do not contribute to the entropy of Γ N , its free energy with respect to the original Hamiltonian (1.1) is given by (2.11) Here V ij denotes the interaction between Regions i and j. The remaining part of this section will be devoted to finding an appropriate upper bound to the right-hand side of Eq. (2.11), which is an upper bound for F (β, N, ω) by the Gibbs variational principle (1.20). In order to simplify the notation, we will from now on replace R + ℓ everywhere by R. Since ℓ ≪ R the number ℓ does not enter the proofs explicitly, except when we bound the interaction energy between the condensate and the thermal cloud.

Preparatory lemmas
The thermal cloud in Region 3, that is, in B(R) c , is described by an ideal Bose gas with the oneparticle We would like to relate its canonical free energy to that of the gas living in all of R 3 and without a cut-off. To that end, we will first compare it to the grand canonical free energy with the help of Corollary A.1. Here explicit formulas are available which allow us to quantify the change in energy caused by the Dirichlet boundary conditions at ∂B(R) and by the energy cut-off Λ. As a preparation we state and prove in this subsection four Lemmas.
The first one, Lemma 2.1, is a general statement allowing to compare traces of functions of Schrödinger operators with different boundary conditions. Lemma 2.2 estimates the differences between traces of functions of Schrödinger operators with Neumann boundary conditions acting on L 2 (B(R) c ) and those acting on L 2 (R 3 ) without boundary conditions. Together these two lemmas are used in the sequel to quantify the difference between the grand canonical free energy of the system living in Region 3 and that of the system living in R 3 , and also the difference of the expected particle numbers. Lemmas 2.1 and 2.2 also enter the proof of Lemma 2.3 which concerns the effects of the boundary condition at ∂B(R) and the cut-off Λ on the chemical potential.
To show that the interaction energy in the thermal cloud and between the thermal cloud and the condensate is of lower order we need an estimate on the L ∞ -norm of the density of the canonical ideal gas in Region 3. Using Proposition A.2, we can estimate the canonical density in terms of the grand canonical density. To close the argument, we need a bound on the L ∞ -norm of the latter showing that the system is dilute in a suitable sense. Such a bound is given in Lemma 2.4 whose proof uses also Lemma 2.3.
with Neumann/Dirichlet boundary conditions at ∂B(R) and choose µ with µ ≤ Cω for some C > 0. We then have (2.13) Proof. By the assumptions on f we can write where We further assume that j 1 (x) equals one for x ∈ B(R) and zero for x ∈ B(2R) c , as well as An application of the IMS localization formula (see e.g. [6]) and the inequality holds. Here χ i denotes the characteristic function of the support of j i and h N,D ≥R,≤2R is the harmonic oscillator Hamiltonian in the annulus A(R, 2R) with Neumann boundary conditions on ∂B(R) and Dirichlet boundary conditions on ∂B(2R). To arrive at Eq. (2.15), we used that the cut-off functions j i introduce additional Dirichlet boundary conditions for the operators under the trace. Together with Eq. (2.14) and holds.
Let us continue with the first term on the right-hand side of the above equation. Using the convexity of f , we see that In order to give an upper bound for the second term on the right-hand side of Eq. (2.17), we note that for x → 0, we see that the first term on the right-hand side of Eq. (2.18) is at most of the order (βω 2 R 2 ) −3 . Our assumptions on R imply (βω 2 R 2 ) −3 ≪ (βω) −3 . To treat the second term, we recall that the eigenvalue ωn of h is g(n)-fold degenerate, where g(n) = (n + 2)(n + 1)/2. Hence, for an appropriately chosen integer m 0 > 0, we can write the second term on the right-hand side of Eq. (2.18) as To obtain the last line, we used that the sum in the line above can be interpreted as a Riemann sum approximating the integral in the last line, as well as the bound e α (h) − µ − 3R −2 ≥ ω. We now collect the results of Eqs. (2.17)- (2.19) and obtain as an upper bound on the first term on the right-hand side of Eq. (2.16).
It remains to give a bound on the second term on the right-hand side of Eq. (2.16). Using the Weyl asymptotics [45, Satz XI] for the eigenvalues of the Laplacian in A(R, 2R), one sees that holds for some appropriately chosen constant C > 0. This allows us to estimate the second term on the right-hand side of Eq. (2.16) by In combination, Eqs. Lemma 2.2. Let f : R + → R be a nonnegative and monotone decreasing measurable function. We also assume that Proof. Denote by h N ≤R the harmonic oscillator Hamiltonian inside B(R) with Neumann boundary conditions. Using the Weyl asymptotics [45, Satz XI] of the eigenvalues of the Neumann Laplacian −∆ N ≤R inside B(R), we see that there exists a C > 0 such that for α ≥ 0. Choose α 0 to be the smallest positive integer for which the right-hand side of Eq. (2.25) is larger than or equal to ω + µ. Our assumption µ ≤ cω and the monotonicity of f , the trace on the right-hand side of Eq. (2.24) can be bounded from above by The asymptotic behavior (2.23) of f at 0 and µ ≤ cω with c < 1 imply that the first term on the right-hand side of the above equation can be bounded from above by Next, we investigate the second term on the right-hand side of Eq. (2.26). The α 0 -th eigenvalue of h N R is bounded from above by To obtain the second inequality, we used that the eigenvalue ωn of h is (n+1)(n+2)/2-fold degenerate.
On the other hand, e α (h N ≥R ) ≥ ω 2 R 2 4 − 3ω 2 which is much larger than the right-hand side of Eq. (2.27) by the assumption ω −1/2 ≪ R. Hence, the second term on the right-hand side of Eq. (2.26) can be written as It remains to estimate the second term on the right-hand side of Eq. (2.28). We invoke Eq. (2.25) and the monotonicity of f to find To arrive at the second line, we used µ ≤ cω and the fact that the sum in the first line on the right-hand side can be interpreted as a Riemann sum approximating the integral in the second line. This proves the claim.

Lemma 2.3. Define µ via the equation
holds. Here N gc th denotes the expected number of particles in the grand canonical thermal cloud, defined in (1.31).
Proof. Since e α (h) ≤ e α (h D ≥R ), the chemical potential can only increase if we cut out B(R) and impose Dirichlet boundary conditions at ∂B(R). Since the cut-off Λ also causes the chemical potential to increase, we have µ ≥ µ 0 .
In order to find an upper bound on µ − µ 0 , we start by estimating the influence of the cut-off Λ and write To be able to proceed, we need a rough upper bound on the chemical potential. Certainly it cannot be larger than the lowest eigenvalue of h D ≥R . A simple trial state argument gives an upper bound ω 2 R 2 on this lowest eigenvalue and thereby on µ.
holds. The trace on the right-hand side of the above equation can be bounded by a constant times (βω) −3 , and hence Tr 1 Let us have a closer look at the first term on the right-hand side of Eq. (2.36). By convexity of Using that µ 0 < 0, we know from Lemma 2.2 with the choice f (x) = (e x − 1) −1 that holds. Note that the trace on the right-hand side of the above equation equals N gc th ∼ (βω) −3 , the expected number of particles in the grand canonical thermal cloud. It dominates the error term in (2.38) the one-particle density in the grand canonical ensemble in Region 3 with cut-off Λ and let the chemical potential µ be chosen as in (2.30) such that Assume that e −βΛ/4 βω and βω ln N 1 holds. Then, for Proof. We start by noting that we obtain an upper bound on ̺ D,Λ ≥R (x) if we drop the cut-off Λ on the right-hand side of Eq. (2.39). To be able to continue, we need an upper bound on the chemical potential µ which we are going to derive now. Lemma 2.3 tells us that µ − µ 0 is bounded from above by the right-hand side of Eq. (2.31). Since µ 0 < 0 the same expression also bounds µ, that is, The first term on the right-hand side of Eq. (2.42) is much smaller than ω because R ≫ ω −1/2 . The second term is smaller than λω 2 R 2 because R ≤ λω −1 β −1/2 . For the third term we have β −1 e −βΛ/4 ≪ ω 2 R 2 by assumption and with Lemma A.2 the fourth term can be estimated by . This is much smaller than λω 2 R 2 because of ω 2 R 2 ≫ ω and our assumption βω ln N 1. We conclude that µ λω 2 R 2 holds.
By choosing λ small enough, the upper bound on µ above allows to conclude that V (x) ≥ 0 holds for |x| ≥ R. Using the Feynman-Kac formula, see e.g. [41,4], we have for all x, y ∈ B(R) c . Here dW x,y denotes the Wiener measure on paths with startpoint x and endpoint y and Ω is the set of those paths that do not leave B(R) c . By 1 Ω (q) we denote its characteristic function. Since V (x) ≥ 0, the exponential function in the second line on the right-hand side of Eq. (2.43) is bounded from above by 1. We can also drop the characteristic function of Ω in the Wiener integral to obtain an upper bound. Together with Eq. (2.43) this implies for all x, y ∈ B(R) c . Here e βα∆ (x, y) denotes the heat kernel of the Laplacian on R 3 which is bounded from above by (4πβα) −3/2 , see e.g. [24]. Hence, which proves the claim.

The thermal cloud
With the above preparations at hand we start our discussion of the free energy of the trial state in Eq. (2.11) by considering the part representing the energy of the thermal cloud. In terms of the spectral decomposition of the density matrixΓ D,Λ ≥R this energy can be written as Bearing in mind that all eigenfunctions Ψ α of H D ≥R can be chosen to be real-valued, we integrate by parts once to rewrite the kinetic energy for the i-th coordinate as where dX is short for d(x 1 , . . . , x N th ). For the energy of a single Ψ α this implies and the whole energy can be written as 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.49) as long as (ω + Λ) 3/2 a N b 3 N 2 th is small enough. We need the cut-off Λ in the definition of the stateΓ D,Λ ≥R in the proof of Lemma 2.5. It could be avoided if we knew that the L 4 (R 3 )-norm of the eigenfunctions of the operator h D ≥R are bounded independently of the energy, which we expect to be true. (Compare with the result in [19] for h on the whole space R 3 .) The bound would most likely grow with R, however, which would need to be quantified. We do not have such a bound at our disposal, and therefore need the cut-off. It will be chosen such that ω ≪ T ≪ Λ holds. Lemma 2.5. There exists a constant C > 0 independent of α such that Proof. Spelled out in more detail, the norm of F Ψ α reads where ̺ Ψα (x, y) denotes the two-particle density of Ψ α . We use the fact that the Ψ α are symmetrized products of one-particle orbitals to conclude that ̺ Ψα (x, y) ≤ ̺ Ψα (x)̺ Ψα (y) holds, where ̺ Ψα is the one-particle density of Ψ α . This allows us to bound the integral on the right-hand side of Eq. (2.52) in the following way: To obtain the bound for the integral of η b , we used its explicit form and the lower bound f 0 (|x| Together with Eqs. (2.52)-(2.53) this proves the claim.
Next we analyze the numerator of the second term on the right-hand side of Eq. (2.49). 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.49) and we start with the first term on the right-hand side of the above equation. Introducing the function Ψα (x, y) = ̺ Similarly, the off-diagonal terms in Eq. (2.57) are bounded from above by Combining this with Eqs. (2.49)-(2.50), we obtain for an appropriately chosen constant C > 0 and C(ω + Λ) 3/2 a N b 2 N 2 th < 1 as an upper bound for the energy of the thermal cloud, where With the help of Proposition A.2 and Lemma 2.4, one readily estimates the first term on the righthand side of Eq. (2.61) by The terms below the curly brackets are obtained from the explicit form of f b , see [28,Appendix C].
To obtain the final bound we used b ≥ ca N for some c > 1. Note that we assumed e −βΛ/4 βω and βω ln N 1 in order to be able to apply Lemma 2.4.
The second term on the right-hand side of Eq. (2.61) can be treated with a rough bound that we derive now. Let a α and a * α denote the usual creation and annihilation operators on Fock space corresponding to an eigenfunction ϕ D α of h D ≥R . An application of the Cauchy-Schwarz inequality tells us that where S 3 denotes the group of permutations of three elements. Using abc ≤ 1 (2.64) We insert this bound into the second term on the right-hand side Eq. (2.61) and obtain Combining this with Eqs. (2.60)-(2.62) we finally obtain for some appropriately chosen C > 0 as an upper bound on the energy of the thermal cloud.
Next, we investigate the entropy of the thermal cloud. To that end, we use [39, Lemma 2] which we spell out here for the sake of completeness.
What remains to be done is to get rid of the Dirichlet boundary condition and the cut-off in the canonical free energy F D,Λ ≥R (β, N, ω) = Tr H D ≥R Γ D,Λ ≥R − T S(Γ D,Λ ≥R ) of the ideal gas. For that purpose, we will first relate it to its grand canonical analogue where explicit formulas are available. Using Corollary A.1, we can bound the canonical free energy from above by where the chemical potential µ is chosen as in (2.30) such that the particle number of the grand canonical Gibbs state associated to 1(h D ≥R ≤ Λ)h D ≥R equals N th . In order to replace the chemical potential µ by µ 0 we use the convexity inequalities To bound the last term from above, we can drop the projection 1(h D ≥R ≤ Λ). Moreover, since µ λω 2 R 2 (as argued in the proof of Lemma 2.4) and h D ≥R for a suitable constant c > 0 for small λ and large ωR 2 . Since the eigenvalues of h D ≥R are larger than the ones of h, this implies that Next, we estimate the influence of the cut-off. To do so, we use the same argumentation as the one in the proof of Lemma 2.3 that leads to Eq. (2.35). We find that In the final step we relate the right-hand side of Eq. (2.75) to the canonical free energy F 0 (β, N, ω) of the ideal gas. First of all, we note that we can drop the spectral projection 1 (h ≥ ω) in the first term on the right-hand side of Eq. (2.75) at the cost of an error of the size β −1 ln(1 − e βµ 0 ) β −1 ln N gc 0 ≤ β −1 ln N . Secondly, we add the missing term µ 0 N 0 which is O(T ). Finally, Corollary A.1 tells us that |F 0 (β, N, ω) − F gc 0 (β, N, ω)| ≤ T O(ln N ). Together with Eq. (2.75), we conclude that Recall that µ − µ 0 was estimated in Lemma 2.3. To obtain the result, we assumed ω −1/2 ≪ R ≤ λω −1 β −1/2 with λ > 0 small enough, e −βΛ/4 βω and βω ln N 1 as well as that Λ 3/2 a N b 2 N 2 th is small enough.

The condensate energy and the interaction between the condensate and the thermal cloud
We recall that E D ≤R denotes the ground state energy of the Hamiltonian H D ≤R in (2.2). In the following we write it as E D ≤R (N 0 ) to explicitly highlight its dependence on the number of particles in the condensate. The strategy here is to use existing results for the ground state energy to relate E D ≤R (N 0 ) to the GP energy E GP (N 0 , a N , ω) defined in Eq. (1.6).
If we go through the proof of the upper bound in [29], we obtain Here E GP,D ≤R (N 0 , a N , ω) denotes the GP energy when we minimize only over functions in H 1 0 (B(R)), that is, over functions that vanish outside B(R). It is therefore sufficient to find an upper bound for E GP,D ≤R (N 0 , a N , ω) in terms of the GP energy E GP (N 0 , a N , ω) without additional boundary conditions. Let j 1 , j 2 ∈ C ∞ (R 3 ) be a partition of unity in the sense that j 1 (x) 2 + j 2 (x) 2 = 1 for all x ∈ R 3 . We assume that j 1 equals one for |x| ≤ R/2 and zero for |x| ≥ R and that |∇j 1 (x)| 2 + |∇j 2 (x)| 2 ≤ 12R −2 . The IMS localization formula (see e.g. [6]) tells us that h = j 1 hj 1 + j 2 hj 2 − |∇j 1 | 2 − |∇j 2 | 2 where h is the harmonic oscillator Hamiltonian (1.30). For the GP energy this implies By minimizing over j 1 φ and j 2 φ separately, keeping the constraint j 1 φ 2 + j 2 φ 2 = N 0 , we obtain a lower bound on the right-hand side of Eq. (2.79). Since a N ∼ N −1 and ω 2 R 2 ≫ ω, one easily sees that the minimum is attained if we put all L 2 -mass into the function j 1 φ. For the energy, we therefore obtain This finally proves which is the desired bound for the first term on the right-hand side of (2.11).
The last term in (2.11) to consider is the interaction between the condensate and the thermal cloud, given by where ̺ΓD,Λ ≥R+ℓ is the one-particle density of the stateΓ D,Λ ≥R+ℓ defined in (2.8). Note that we have inserted the missing ℓ in B(R + ℓ) again (compare with the discussion in Subsection 2.1). When we use that F ≤ 1 and apply Lemma 2.5, we find for some C > 0 and C(ω + Λ) 3/2 a N b 2 N th < 1.

An application of Proposition A.2 and Lemma 2.4 tells us that
As already mentioned at the beginning of Section 2.1, we choose ℓ such that ℓN ω 1/2 is larger than the radius of the hard-core part of the interaction potential. We thus conclude that Tr (V 13 Γ N ) ω −1/2 β −3/2 . In combination with Eq. (2.81), this yields the upper bound for the condensate energy plus the interaction energy between the condensate and the thermal cloud. To obtain Eq. (2.85) we assumed in addition to ω −1/2 ≪ R ≤ λω −1 β −1/2 with λ > 0 small enough that e −βΛ/4 βω and βω ln N 1 holds, as required for the application of Lemma 2.4. We also assumed that Λ 3/2 a N b 2 N 2 th is small enough.

The final estimate for the upper bound
We combine the results of Eqs. (2.77) and (2.85). Together with an application of Lemma 2.3 and Corollary A.2, which allows us to obtain a bound on |µ − µ 0 |, we obtain We assumed ω −1/2 ≪ R ≤ λω −1 β −1/2 with λ > 0 small enough, e −βΛ/4 βω, βω ln N 1, b ≥ ca N for some c > 1 as well as that Λ 3/2 a N b 2 N 2 th is small enough. To conclude the proof of the upper bound we will distinguish two cases, one where βω (ln N ) −1 and the other one where (βω) (ln N ) −1 . We start with the first one.
To obtain a better bound for relatively large βω and, in particular, to cover the case βω (ln N ) −1 , we proceed as follows. Let E(N ) denote the ground state energy of H N in Eq. (1.1). We use the upper bound on E(N ) in [29] and estimate ). Moreover, F 0 (β, N, ω) −ω(βω) −4 , see Section 1.10 and Appendix A. We therefore have This completes the proof of the upper bound.

Proof of the lower bound
As for the upper bound, the main idea of the proof of the lower bound is to make use of the two different length scales on which the condensate and the thermal cloud live. Anticipating the result that their respective particle numbers are equal to those of the ideal gas to leading order, this implies that the thermal cloud is much more dilute than the condensate and therefore does not see the interaction to the same order as the condensate does. For a more detailed discussion of these issues, see Section 1.5. The main technique to implement this mathematically is geometric localization in Fock space which has been introduced in [9] in the context of bosonic quantum field theory. It allows, for the purpose of a lower bound, to replace the free energy of the whole system by a sum of two free energies, one of a system localized in a ball with radius 2R and another one of a system living in the complement of a ball with radius R. At this point the two systems are still correlated, however. The overlap of the two regions comes from the fact that we have to use smooth cut-off functions.
The radius R is, as in the proof of the upper bound, chosen such that ω −1/2 ≪ R ≪ ω −1 β −1/2 holds. In the following step, we minimize these two free energies separately which again results in a lower bound. For this it is necessary to drop the restriction on the particle number and work in Fock space. The minimization procedure in the ball with radius 2R yields the GP energy (as for the upper bound we use the known result at zero temperature here) plus lower order corrections coming from the entropy. In the complement of the ball with radius R we drop the interaction (it is positive by assumption) and obtain the free energy of an ideal gas. Throughout the proof of the lower bound we shall assume that (βω) −1 N 1/3 (i.e., T T c ) as well as that Rω 1/2 is large enough.
The case where N 0 = o(N ) and/or N 0 a N = o(1) is quite simple for the lower bound so we discuss it first. On the one hand, the GP energy is of order o(ωN ) in this case. On the other hand, our interaction potential is nonnegative which allows us to drop it to obtain a lower bound on the free energy, i.e., F (β, N, ω) ≥ F 0 (β, N, ω), which is sufficient in this case. The case where N 0 ∼ N is considerably more difficult and will be treated in the remaining part of this Section.
Our first task is to replace the free energy of a given state Γ N by the sum of the free energies of two localized versions of Γ N and to show that this yields a lower bound. Let j 1 , j 2 ∈ C ∞ (R 3 ) be a partition of unity in the sense that j 1 (x) 2 + j 2 (x) 2 = 1 for all x ∈ R 3 and choose j 1 such that it equals one for x ∈ B(R) and zero for x ∈ B(2R) c . We can also assume that |∇j 1 (x)| 2 + |∇j 2 (x)| 2 ≤ 3R −2 holds. By γ (k) Γ we denote the k-particle reduced density matrix of a state Γ on the bosonic Fock space F(H). It is defined via the integral kernel γ (k) Γ (y 1 , . . . , y k ; x 1 , . . . , x k ) = Tr a * x 1 · · · a * x k a y k · · · a y 1 Γ . (3.1) The following well-known Lemma, see [9,16], concerns localizations of a given state. Because its proof is instructive we sketch it here.
Lemma 3.1. Let Γ be a state on F(H) with k-particle density matrices γ (k) Γ , k ≥ 1 and j 1 , j 2 as above. Then there exist unique states Γ j 1 and Γ j 2 on F(H) with k-particle density matrices j ⊗k 1 γ (k) Γ j ⊗k 1 and j ⊗k 2 γ (k) Γ j ⊗k 2 , respectively. Moreover, the entropies of these states are related by Proof. Define the map J : H → H ⊕ H, ψ → (j 1 ψ, j 2 ψ), (3.3) where H = L 2 (R 3 ). That is, J splits the wavefunction into two parts, one supported in B(2R) and one supported in B(R) c . Let A be an operator on H. We denote by Υ(A) the second quantization of A acting on F(H). On basis vectors its action is defined by where Ω denotes the Fock space vacuum. It extends to all Fock space vectors by linearity. Let U be the unitary map that identifies F(H ⊕ H) and F(H) ⊗ F(H) given by Here {ψ i } ∞ i=0 denotes a basis for the first copy of H and {ϕ i } ∞ 0=1 is a basis for the second copy. By c * /d * we denote the creation operator acting on the first/second factor of F(H) ⊗ F(H) and Ω 1 /Ω 2 is the vacuum in the first/second factor. Having those two operators at hand, we can define the J-extension Γ J of Γ by It follows that the states Γ j 1 = Tr 2 (Γ J ) and Γ j 2 = Tr 1 (Γ J ) (3.7) have the desired property. Here Tr 1,2 denotes the partial trace over the first/second factor of F(H) ⊗ F(H). The uniqueness part follows from the fact that states are uniquely determined by their reduced density matrices, which we are not going to discuss here. The proof of Eq. (3.2) follows from Eq. (3.7) and the subadditivity of the entropy.
Before we continue, let us introduce some notation. For an operator A on H let the operator dΥ(A) on F(H) be given by dΥ(A)| F N = N i=1 A i where A i stands for A acting on the i-th particle in the the Fock space sector F N with N particles. We also define on F(L 2 (B(2R)), where h D ≤2R denotes the operator (1.30) restricted to B(2R) with Dirichlet boundary conditions. ByN we denote the particle number operator on Fock space.
Let Γ N ∈ S N be an N -particle state. Using the IMS localization formula (see e.g. [6]) and the fact that v(x) ≥ 0 for all x ∈ R 3 , we find (3.9) Lemma 3.1 tells us that we can bound the entropy of Γ N by the ones of Γ j 1 and Γ j 2 , the localized states related to Γ N . It also allows us to write the energies in Eq. (3.9) in terms of Γ j 1 and Γ j 2 . This implies ≥R the operator h in (1.30) restricted to B(R) c with Dirichlet boundary conditions. The states Γ j 1 and Γ j 2 are related via the fact that they both are constructed by localizing the state Γ N . In the next step we will minimize over each of them separately which results in a lower bound. We will also drop the restriction on the particle number. To do so, we have to introduce a chemical potential (and subtract it again).
The intuition from Section 1.5 tells us that only the condensate is affected by the interaction to leading order. Since the interaction energy per particle inside the condensate can be expected to be of order ω, while −µ 0 ∼ β −4 ω −3 is much smaller when N 0 ∼ N (see Section 1.10), the chemical potential µ for the interacting system will be of order ω, too. In fact, in the GP limit the GP chemical potential µ GP = µ GP (N 0 , a N , ω), defined by is the correct choice for µ. We note that because of scaling (see Eqs. (1.3) and (1.7)) µ GP = ωµ GP (1, a v , 1). In particular, µ GP ∼ ω for fixed a v > 0. The chemical potential µ GP necessarily also appears in the thermal cloud, but as we will see below, this does not affect its free energy at the level of accuracy we are interested in.
Using the explicit form of the one-particle density matrices of Γ j 1 and Γ j 2 (see Lemma 3.1), we check that Tr(γ Γ j 1 + γ Γ j 2 ) = N holds. Together with Eq. (3.10), this implies From here on we estimate the two contributions on the right-hand side of Eq. (3.12) separately, and start with the first line.

The condensate
For the minimization problem inside B(2R) we take a small amount of the kinetic energy to control the entropy, which results in a lower order contribution. The minimization problem for the remaining part of the energy then follows from the results in [29] for the ground state energy. Let 0 < ǫ < 1.
From the positivity of the interaction potential v, we conclude that We will later choose ǫ such that ǫ ≪ 1 holds.
We start by considering the first term on the right-hand side of Eq. (3.13 (3.14) To bound the infimum on the right-hand side of the above equation, we distinguish two different regimes for the particle number M . This is necessary because we cannot relate E(M ) to the GP energy (1.6) if M is too small.
We first consider the case where M ≪ N 0 . Choose δ > 0 such that δ ≪ 1 and still δN 0 ≫ 1 and assume M ≤ δN 0 . When we drop the positive energy E(M ), we obtain the lower bound for the energy inside B(2R). Next, we estimate the contribution coming from the entropy in the region B(2R). We use Tr[N Γ j 1 ] ≤ N to see that the term in the second line of Eq. (3.13) is bounded from below by (3.17) Note that we have added and subtracted a chemical potential of size − 3 2 ω. To bound the first term on the right-hand side, we note that the h D ≤2R is bounded from below by −∆ D ≤2R − 3 2 ω. Using the Weyl asymptotics [45, Satz XI] of the eigenvalues of −∆ D ≤2R , we see that there exists a constant C > 0 such that the αth eigenvalue of this operator satisfies e α (−∆ D ≤R ) ≥ C(α + 1) 2/3 /R 2 . This allows us to estimate The last sum can be interpreted a Riemann sum approximating the corresponding integral, and one readily checks that it is bounded below by −O(R 3 β −3/2 ǫ −3/2 ). That is, Together with Eqs. (3.13), (3.16) and (3.17) this yields as a lower bound for the free energy inside B(2R).

The thermal cloud
Next we consider the terms in the second line on the right-hand side of Eq. (3.12). Explicit minimization shows that To relate the right-hand side of Eq. (3.21) to F 0 (β, N, ω), we first need to replace µ GP by µ 0 , the chemical potential of the ideal gas leading to an expected number of N particles. After the chemical potential has been replaced, we have to get rid of the Dirichlet boundary conditions in the formula for the grand canonical free energy, and replace it by its canonical version.
To replace µ GP by µ 0 , we use convexity to bound By arguing as in Eq. (2.72), we see that the last trace is bounded by O((βω) −3 ). It remains to get rid of the Dirichlet boundary conditions in the various terms in Eq. (3.22). To that end, we use that the α-th eigenvalue of h D ≥R is bounded from below by e α (h D ≥R ) ≥ e α (h). For the second term on the left-hand side of (3.22), this implies Note that the above bound is rough in the sense that the right-hand side includes the grand canonical potential of the condensate. The latter is negligible in the limit considered, however. The first term on the right-hand side of Eq. (3.22) is proportional to the particle number and hence we have to be more careful. We use e α (h D ≥R ) ≥ e α (h) for α ≥ 1 and e 0 (h D ≥R ) ≥ ω 2 R 2 4 − 3ω 2 , which avoids adding the expected number of particles in the condensate, and find Tr 1 Note that the first term on the right-hand side of Eq. (3.24) is N gc th , the expected number of particles in the grand canonical thermal cloud. Combining Eqs. (3.21)-(3.24) and using that |µ 0 − µ GP | ω we find It remains to replace the grand canonical free energy by the canonical one, and N gc th by N th . Corollary A.1 tells us that the difference of the canonical and the grand canonical free energy is at most of order −T ln N . Moreover, |N th − N gc th | (βω) −3/2 (ln N ) 1/2 + (βω) −1 ln N by Corollary A.2. We therefore have as a lower bound for the free energy in region B(R) c . Note that we have added the additional negative term µ 0 N gc 0 on the right hand side.

The final estimate for the lower bound
We combine the results from Eqs. (3.12), (3.20) and (3.26) and the fact that To obtain the result we assumed that Rω 1/2 is large enough, and used that (βω) −1 N 1/3 to dominate some of the error terms by others. The optimal choice of the parameters R and ǫ turns out to be R ∼ ω −1/2 N 1/8 (βω) 5/16 and ǫ ∼ N −1/4 (βω) −5/8 . The three terms on the second line of the right hand side of Eq. (3.27) are thus bounded by ωN 23/24 . In particular, (3.28) This completes the proof of the lower bound.

Proof of the asymptotics of the one-particle density matrix
In the following discussion, we assume that N 0 ∼ N and N a N ≥ ǫω −1/2 , i.e. a v ≥ ǫ for some ǫ > 0 holds. The last assumption is necessary keeping in mind the second remark in Section 1.8. The case where one of these conditions is not fulfilled will be taken care of at the end.
Assume we are given a sequence of states Γ N with reduced one-particle density matrices γ N such that as N tends to infinity. Choose two functions j 1 and j 2 as in the proof of the lower bound, satisfying j 1 (x) 2 + j 2 (x) 2 = 1 for all x ∈ R 3 , j 1 (x) = 1 for x ∈ B(R) and j 1 (x) = 0 for x ∈ B(2R) c , and also |∇j 1 (x)| 2 + |∇j 2 (x)| 2 ≤ 3R −2 . Using Eqs. (3.12), (3.13), (3.20) and (3.26) and the choice of parameters from Section 3.3, we see that holds. The operator H D ≤2R was defined in Eq. (3.8), µ GP is given by (3.11) and the state Γ N,j 1 is related to Γ N in the way described in the proof of Lemma 3.1. Moreover, Eqs. (3.12) and (3.20) together with the fact that µ 0 N 0 < 0 tell us that where also Γ N,j 2 is related to Γ N in the way described in the proof of Lemma 3.1. Note that we have the grand canonical free energy on the right-hand side of Eq. (4.3) instead of the canonical free energy, which is allowed by Corollary A.1 in the Appendix.
Eqs. (4.2) and (4.3) will be used to deduce the desired bounds on γ N . As a first step we will derive asymptotic expressions for the one-particle density matrices of Γ N,j 1 and Γ N,j 2 , that is, for j 1 γ N j 1 and j 2 γ N j 2 , respectively. Afterwards, we consider the "off-diagonal" contribution coming from j 1 γ N j 2 and j 2 γ N j 1 and show that their trace norm is of order o(N ). As one would expect, j 1 γ N j 1 turns out to be close to To derive a bound for j 1 γ N j 1 , we make use of existing results [25,26,38] on the convergence of the one-particle density matrix of approximate minimizers of the ground state energy functional to the projection onto the GP minimizer. The main difficulty to overcome is that the particle number of the state Γ N,j 1 may fluctuate, that is, it is a state on the full Fock space.
As in Section 3.1 we choose 0 < δ ≪ 1 such that still δN 0 ≫ 1. Let P M be the projection onto the Fock space sector with M particles. Keeping in mind the normalization Tr[Γ N,j 1 ] = 1, Eq. (4.2) can be written as Let us again distinguish two cases: For 0 ≤ M ≤ δN 0 we drop H D ≤2R to obtain a lower bound and use µ GP δN 0 ≪ ωN 0 to show that there is a constant C 1 > 0 such that the expression in the curly brackets in the above equation is bounded from below by C 1 ωN 0 . In the case where δN 0 < M ≤ N , we invoke [38,Eq. (101) and Lemma 4] to see that there exists a constant C 2 > 0 such that the summand in Eq. (4.4) is bounded from below by Here γ P M Γ N,j 1 P M denotes the one-particle density matrix of the state P M Γ N, It is normalized to have Tr[γ P M Γ j 1 P M ] = M . By P GP M a N we denote the projection onto the GP minimizer φ GP 1,M a N . This estimate is in fact uniform in M a N for M in the range we consider.
We shall use the strict convexity of M → E GP (M, a N , ω) in order to obtain a lower bound on the first three terms in the above parentheses that is strictly positive for M = N 0 . Using that µ GP = d dN E GP (N 0 , a N , ω), as well as the convexity of ρ → |∇ √ ρ| 2 , we deduce that For a lower bound, we pick s > 0 and t ∈ R and estimate After optimizing over s and t, we thus obtain the lower bound 1 14π In particular, with (4.6) we conclude that Putting all this together, we obtain Next, we write j 1 γ N j 1 = N M =0 Tr [P M Γ N,j 1 P M ] γ P M Γ N,j 1 P M and estimate the trace norm difference of j 1 γ N j 1 and N 0 P GP N 0 a N in a first step by Eq. (4.12) tells us that the first two terms on the right-hand side of the above equation are of order o(N ). For the term in the last line, we insert N 0 P GP M a N − N 0 P GP M a N in the obvious place to see that it is bounded from above by To bound the right-hand side of Eq. (4.14), we choose 0 < κ < 1 and split the sum into two parts, one where |M − N 0 | ≤ κN 0 and another one where |M − N 0 | > κN 0 . We claim that there exists a function f : R + → R + with f (x) → 0 for x → 0 such that the first part of the sum is bounded by (κ + f (κ))N 0 . This follows from the continuity of the map M → φ GP 1,M a N in L 2 (R 3 ), which can easily be deduced from the uniqueness of the minimizer of the GP functional. To estimate the contribution to the sum of the terms with |M − N 0 | > κN 0 , we write Together with Eq. (4.12), this implies for κ 1 that Choosing κ ≪ 1 and δ appropriately, we see that the right-hand side of Eq. (4.16) as well as the part of the sum in Eq. (4.14) where |M − N 0 | ≤ κN 0 is of the order o(N ). Together with Eq. (4.13), this proves The main ingredient to derive a bound for j 2 γ N j 2 is a novel coercivity estimate for the bosonic relative entropy that we prove in Lemma 4.1 below. Using this estimate, we shall show that j 2 γ N j 2 is close to γ gc 0 − N gc 0 |ϕ 0 ϕ 0 | in trace norm, where denotes the grand canonical analogue of γ N,0 . This part of the proof is motivated by a related analysis for the one-particle density matrix of a dilute Fermi gas in [39].
For positive trace-class operators γ define We have [44, 2.5.14.5] Since we conclude that holds. Let us define ν(x) = max{x, µ} with cµ GP < µ ω for some c > 1. For what follows, it will be convenient to replace the Hamiltonian h by ν(h) on the right-hand side of Eq. (4.22). We write h = ∞ α=0 e α (h)|ϕ α ϕ α | and choose α 0 to be the largest integer such that e α 0 (h) < µ. Using α 0 = O(1) and µ ω we can estimate where C > 0 is some appropriately chosen constant. To obtain the last inequality on the right-hand side of Eq. (4.24), we used the decay of the eigenfunctions of h and the fact that the support of j 2 is given by B(R) c . Together with Eq. (4.3), Eq. (4.22) and ω 1/2 R ≫ 1, this implies To be able to compare the expressions in the first and in the second line on the right-hand side of Eq. (4.25), we will replace h by ν(h) and afterwards µ 0 by µ GP in the first term in the second line. In fact, since ν(h) ≥ h, Moreover, using that βµ 0 = O(N −1 0 ) and thus |µ 0 − µ GP | ω, we can proceed similarly to Eq. (3.22) to obtain holds.
The right-hand side of Eq. (4.28) can be written in terms of the relative entropy, which is defined as follows. For two nonnegative operators γ, γ 0 with finite trace, the bosonic relative entropy of γ with respect to γ 0 is given by with σ defined in Eq. (4.19). For matrices it is well-defined as long as γ 0 > 0. In case γ 0 has a nontrivial kernel and γ − γ 0 = 0 on ker(γ 0 ), one defines S(γ, γ 0 ) = ∞. If γ − γ 0 = 0 on ker(γ 0 ) the trace is by definition taken on the complement of that subspace. In the case of trace-class operators with γ 0 strictly positive, one can equivalently define where S(x, y) = σ(x) − σ(y) − σ ′ (y)(x − y) ≥ 0, and {λ i , ψ i } respectively {ν j , ϕ j } are the eigenvalues and eigenfunctions of γ and γ 0 , respectively. The definition (4.30) will be most convenient for our purpose. We note that S can be defined more generally even for non-compact operators by approximating the operators by matrices and taking limits, see [22,10].
Denote by the one-particle density matrix of the Gibbs state related to the grand canonical potential 1 β Tr[ln(1 − e −β(ν(h)−µ GP ) )]. A simple computation shows that Eq. (4.28) can equivalently be written as In order to get quantitative information out of Eq. (4.32), we need the following Lemma: There exists a constant C > 0 such that for any two nonnegative trace-class operators γ, γ 0 we have . Proof. We start with the proof of the second inequality. With S(x, y) = σ(x) − σ(y) − σ ′ (y)(x − y) we wish to show that S(x, y) ≥ C (x − y) 2 (x + y)(1 + y) (4.35) holds for all numbers x, y ∈ R + . To that end, we will consider several cases and start with the one where x ≥ y and y ≥ 1. We claim that which follows from s(s + 1) ≤ 2s 2 for s ≥ 1. One also checks that z ≥ 1 implies −1 + z − ln(z) ≥ (8/9)(z − 1) 2 /(z + 1). Together with Eq. (4.36), this gives . (4.37) Next, consider the case x ≥ 2 and y ≤ 1. Here We use the same inequality as above to obtain a lower bound for the right-hand of Eq. (4.38) and find For x ≥ y and x ≤ 2, we again consider the integral representation of S(x, y) from above and replace s by x + y in the denominator of the function under the integral sign. This yields a lower bound of the form . (4.40) It remains to consider the case y ≥ x. Here we argue in the same way as in the previous step: This proves the claimed bound (4.35) for S(x, y) with C ≥ 2/27.
In order to deduce Eq. (4.34) from Eq. (4.35), we follow the argument in the proof of [39,Lemma 7].
Let us start with the first term on the right-hand side of (4.49). We want to show that it is of the order o(N ). To do so, we bound where e i (γ) are the eigenvalues of γ. Since f (t) ≥ 2f (t/4), we have f (e i (γ)) ≤ 2(f (e i (γ)) − f (λ/4)) for e i (γ) > λ. Since f (γ ν,0 ) ≤ f (λ/4), we further have (denoting by ψ i the eigenfunctions of γ) We know from Lemma 4.1 and Eq. (4.32) that Next, consider the second term on the right-hand side of Eq. (4.49). Using the Cauchy-Schwarz inequality and the cyclicity of the trace, we find for 0 ≤ x, y ≤ λ. The fact that max{ γ ν,0 , γ } ≤ λ and an application of Klein's inequality therefore gives Denote by {ψ j } ∞ j=0 the eigenbasis of γ and by {ϕ i } ∞ i=0 the eigenbasis of γ ν,0 . In order to replaceγ by γ, we write . It follows that f (e j (γ)) ≥ 2f (e i (γ ν,0 )) for all i and j such that e j (γ) > λ. Hence we can replace f (1(e j (γ) ≤ λ)e j (γ)) by f (e j (γ)) in Eq.  (4.60) Having Eq. (4.60) and Eq. (4.59) at hand, we combine them with Eq. (4.49) and γ ν,0 1 ≤ N to finally obtain where we inserted j 2 γ N j 2 for γ. To complete the argument it remains to choose the projection P .

The off-diagonal elements of γ N and the final estimate
To complete the proof of Theorem 1.1, it remains to estimate the trace norm of j 2 1 γ N j 2 2 . In combination with Corollary A.2, which shows that the trace norm difference of γ N,0 and γ gc 0 is small, this will allow us to conclude the convergence result (1.25) for the one-particle density matrix in the case N 0 ∼ N and N a N ≥ ǫω −1/2 for some ǫ > 0. Finally, we shall comment on the case where one of these assumptions is not valid.
We defineγ gc 0 = γ gc 0 − N gc 0 |ϕ 0 ϕ 0 | as well as P GP = P GP N 0 a N = |φ GP 1,N 0 a N φ GP 1,N 0 a N | for short. The identity j 1 (x) 2 + j 2 (x) 2 = 1 for all x ∈ R 3 and the triangle inequality allow us to bound gc 0 j 2 1 (4.66) With Eq. (4.17), Eq. (4.65) and j i (x) ≤ 1 for i ∈ {1, 2}, we see that the first two terms on the right-hand side of the above equation are of order o(N ). To derive an estimate for the first term in the second line of Eq. (4.66), we bound j 2 1 γ N j 2 2 1 ≤ P GP j 2 1 γ N j 2 2 1 + (1 − P GP )j 2 1 γ N j 2 2 1 . Since P GP is a rank one projection, we have (recall that · denotes the operator norm) We also estimate Recall that the largest eigenvalue ofγ gc 0 is bounded by a constant times (βω) −1 . Making use of Eq. (4.65), this implies P GP j 2 1 γ N j 2 2 1 ≤ o(N ). On the other hand, where · 2 denotes the Hilbert-Schmidt norm. We apply Eq. (4.17) to see that the trace in the second line of Eq. (4.69) is bounded by o(N ) + N 0 P GP − j 1 P GP j 1 1 . Moreover, the exponential decay of φ GP 1,N 0 a N , see [29, Appendix A], implies for an appropriately chosen constant c > 0. Since R ≫ ω −1/2 we know that the last term on the right-hand side of Eq. (4.70) is of order o(1). Together with Eqs. (4.67)-(4.70), we therefore see that holds.
It remains to give a bound on the last term on the right-hand side of Eq. (4.66). We add and subtract j 2γ gc 0 and use j 2 ≤ 1 to see that holds. The last factor equals (N gc th ) 1/2 . The square of the first factor, on the other hand, can be bounded by To obtain this bound, we used that the support of 0 ≤ 1 − j 2 ≤ 1 is given by B(2R). We claim that γ gc 0 (x, x) β −3/2 holds. This can be seen by an analysis similar to a part of the analysis carried out in Lemma 2.4: We first replace the chemical potential µ 0 by a larger one that guarantees − 3ω 2 − µ ≥ 0 to hold. In fact, we choose µ = − 3ω 2 . Then there exists a constant C > 0 such that holds for all α ≥ 1. This can easily be checked when we realize that the expression on left-hand side of the above equation is monotone decreasing in e α (h) if µ < µ 0 . It therefore suffices to check the inequality with e α (h) replaced by ω. Eq. (4.74) implies Recall that so far we have worked under the assumptions N 0 ∼ N and N a N ≥ ǫω −1/2 for some ǫ > 0. It remains to consider the case where either N 0 = o(N ) and/or a N ≪ ω −1/2 N −1 . In each of these cases, we have P GP − |ϕ 0 ϕ 0 | ≪ 1, and also E GP (N 0 , a N , ω) ≪ ωN . The equivalent of Eq. (4.1) therefore reads and implies With this input, we go through the analysis of Section 4.2. In case µ 0 −ω, we can directly apply Lemma 4.1 and the subsequent estimates, and the equivalent of Eq. (4.65) tells us that γ N − γ gc 0 1 ≤ o(N ). Together with Corollary A.2 this implies the claim. If |µ 0 | ≪ ω, we have to first remove the condensate. Let P = |ϕ 0 ϕ 0 | denote the projection onto the ground state of h, and Q = 1 − P . From the subadditivity of the entropy and Eq. (4.20), we have To conclude the proof of Theorem 1.1 it remains to prove Eq. (1.26). To that end, we write (4.82) The first term on the right-hand side of Eq. (4.82) can be bounded by the trace norm of the same expression, which is bounded by o(N ) according to Eq. (1.25). Similarly, Eq. (A.24) implies that the second term is bounded by a constant times (N ln N ) 1/2 . From the explicit form ofγ gc 0 we deduce that its largest eigenvalue is of the order O((βω) −1 ). Therefore, and the proof of Theorem 1.1 is complete.

A. Some properties of the ideal Bose gas
In this Appendix we collect several statements about the ideal Bose that are needed in the proof of Theorem 1.1 and do not seem to have appeared in the literature before (except for the first part of Proposition A.1). We start by introducing some notation. Fix a nondecreasing sequence {E j } ∞ j=0 of nonnegative real numbers. A vector n = (n 0 , n 1 , . . .) of infinite length is called a configuration if all its entries are nonnegative integers and if only a finite number of them is different from zero. The collection {E j } ∞ j=0 plays the role of the energy levels of a one-particle quantum system with the temperature factor β absorbed, and a configuration n labels an element of the standard basis of the bosonic Fock space. For each configuration n, we define its energy to be E(n) = j≥0 E j n j . The canonical partition function is given by Z(N ) = |n|=N exp (−E(n)) for N ∈ N 0 . Here |n| = ∞ j=0 n j denotes the number of particles in the configuration n. Let A N be the expectation of an operator A in the canonical Gibbs state related to Z(N ). By a * j and a j we denote the bosonic creation and annihilation operator of a particle with the energy E j and we definen j = a * j a j . The free energy of the system is given by We assume that the energy levels {E j } ∞ j=0 and the function f are such that the partition function and the expectations of all f (n j ) are finite. We then have: Proposition A.1. The map N → F (N ) is convex, i.e., F (N + 1) + F (N − 1) ≥ 2F (N ) holds for all N ≥ 1, and for any nonnegative, nondecreasing function f : N 0 → R the map N → f (n j ) N is nondecreasing for each j ∈ N 0 .
The proof of the first statement seems to appear for the first time in [23] and it was later reproven in [42]. The second statement has been shown in [36] for the functions f (x) = x k with k ∈ N and in [42] for f (x) = x. The proof in [42] still works if one replaces the identity by a nonnegative, nondecreasing function f . Proposition A.1 has two consequences that will be of importance for us. The first is naturally formulated in the following setup: Denote by Z gc (µ) = n exp(−(E(n) − µ|n|)) the grand canonical partition function and define λ N,µ = Z(N ) exp(µN )/Z gc (µ). The grand canonical free energy is given by F gc (µ) = N ≥0 λ N,µ F (N ) − S(λ), where S(λ) = − N ≥0 λ N,µ ln(λ N,µ ). The first consequence of Proposition A.1 is the following statement which quantifies the difference between the canonical free energy and its grand canonical counterpart: Corollary A.1. Assume µ is such that N = N ≥0 N λ N,µ ∈ N. Then The second consequence is an estimate of the canonical one/two-particle density in terms of the grand canonical one-particle density. φ j 1 (x)φ j 2 (y)φ j 3 (y)φ j 4 (x) a * j 1 a * j 2 a j 3 a j 4 N the canonical one-particle and two-particle densities, respectively. The grand canonical one-particle density is given by where µ is chosen such that N = N ≥0 N λ N,µ ∈ N. Then Proof. We start with the proof of the inequality for the one-particle density. From the definition of the one-particle densities, Eqs. (A.4) and (A.5), and the monotonicity of the map N → n j N , see Proposition A.1, we find To prove the claim, we need to show that N ≥N λ N,µ ≥ 1.8 40 holds. From Lemma A.1 below we know that This proves the claim for the one-particle densities (assuming the validity of Lemma A.1).
It remains to prove the bound for the two-particle densities. An application of the Cauchy-Schwarz inequality tells us that holds. Here S 2 denotes the group of permutations of two elements. Let us denote by A gc = N ≥0 λ N,µ A N the expectation of an operator A in the grand canonical Gibbs state. We want to derive an upper bound for the expectation value in the second line on the right-hand side of Eq. (A.13). For i = j we have a * i a * i a i a i =n i (n i − 1) withn i = a * i a i . From Proposition A.1 we know that the map N → n j (n j − 1) N is nondecreasing. Using this fact, we argue as in the case of the one-particle density to see that n j (n j − 1) N ≤ 40 1.8 n j (n j − 1) gc (A.14) holds. The right-hand side of Eq. (A.14) can be simplified when we use Wick's Theorem: n j (n j − 1) gc = 2 n j 2 gc . If i = j we have a * i a * j a j a i =n inj and [42] tells us that n inj N ≤ n i N n j N (A. 15) holds. As in the previous case, we use n j N ≤ 40 1.8 n j gc . Combining these estimates with Eq. (A.13), we finally obtain This proves the claim (A.6).
Remark A.1. The first part of the proof shows that n j N ≤ 40 1.8 n j gc for all j, where · gc denotes the corresponding grand canonical state with average particle number N . In particular, N 0 N gc 0 and N th N gc th holds.
The next Lemma provides an estimate of the fourth moment of the particle number in the grand canonical ensemble in terms of the second moment. It is needed in the proof of Proposition A.2.
Lemma A.1. LetN = ∞ j=0n j be the particle number operator and denote by N = N gc the expected number of particles in the grand canonical ensemble. We then have Proof. If we use that ∂ ∂µ Z gc (µ) = Z gc (µ)N holds, a simple computation leads to which together with ∂ ∂µ N = N 2 gc − N 2 allows us to conclude that (A.20) To treat the first term on the right-hand side, we need to do a little computation. It yields ∂ ∂µ where µ < 0 is chosen such that Tr γ gc N = N ∈ N. We shall also introduceγ gc N = γ gc N − N gc 0 |φ 0 φ 0 |, where N gc 0 = (e −µ −1) −1 is the number of particles in the condensate in the grand canonical ensemble, and similarlyγ N = γ N − N 0 |φ 0 φ 0 |, with N 0 = n 0 N the number of particles in the condensate in the canonical ensemble. Proof. Let Γ G N denote the canonical Gibbs state, with one-particle density matrix γ N . With s(γ) defined in Eq. (4.19), we have the entropy inequality (see, e.g., [44, 2.5.14.5]) S(Γ G N ) ≤ s(γ N ). Therefore, Note that γ N and γ gc N have the same eigenfunctions φ j . In particular, they commute. Hence the relative entropy can be written as where S is defined in the proof of Lemma 4.1, and a α resp. b α denote the eigenvalues of γ N and γ gc N , respectively. Using (4.35) and the Cauchy-Schwarz inequality, this implies The claim (A.24) then follows from the fact that N = Tr γ N = Tr γ gc N , which implies that holds.
We shall now apply Lemma A.2 to the case of the harmonic oscillator Hamiltonian in Eq. (1.30). We re-introduce the inverse temperature β, and adjust the notation to the one used in the main text. That is, we denote the one-particle density matrices by γ N,0 and γ gc 0 , respectively. Recall also the definitions N th = N − N 0 and similarly for N gc th . In this case, we have Hence we obtain the following Corollary.
Corollary A.2. Consider the three-dimensional ideal Bose gas in the harmonic oscillator potential, that is, the one-particle Hamiltonian of the system is given by h = −∆ + ω 2 4 x 2 − 3 2 ω. We assume that the chemical potential µ 0 is chosen such that the expected number of particles in the grand canonical ensemble equals N ∈ N. Then