Classical field theory limit of many-body quantum Gibbs states in 2D and 3D

Abstract

We provide a rigorous derivation of nonlinear Gibbs measures in two and three space dimensions, starting from many-body quantum systems in thermal equilibrium. More precisely, we prove that the grand-canonical Gibbs state of a large bosonic quantum system converges to the Gibbs measure of a nonlinear Schrödinger-type classical field theory, in terms of partition functions and reduced density matrices. The Gibbs measure thus describes the behavior of the infinite Bose gas at criticality, that is, close to the phase transition to a Bose–Einstein condensate. The Gibbs measure is concentrated on singular distributions and has to be appropriately renormalized, while the quantum system is well defined without any renormalization. By tuning a single real parameter (the chemical potential), we obtain a counter-term for the diverging repulsive interactions which provides the desired Wick renormalization of the limit classical theory. The proof relies on a new estimate on the entropy relative to quasi-free states and a novel method to control quantum variances.

Appendices

Appendix A: The counter-term problem

In this appendix we discuss the counter-term problem in detail.

1.1 A.1. Hartree versus reduced Hartree energy functional

We recall that to any one-body density matrix $\gamma \geqslant 0$ one can associate a unique Gaussian state

$$\begin{aligned} \Gamma =\frac{e^{-\mathrm {d}\Gamma (a)}}{{{\,\mathrm{\mathrm{Tr}}\,}}_\mathfrak {F}[e^{-\mathrm {d}\Gamma (a)}]} \end{aligned}$$

on the Fock space which has the one-particle density matrix $\Gamma ^{(1)}=\gamma$ [8, 156]. The unique corresponding one-body operator a is given by

$$\begin{aligned} \gamma =\frac{1}{e^a-1}\quad \Longleftrightarrow \quad a=\log \frac{1+\gamma }{\gamma }. \end{aligned}$$

(A.1)

Then its energy terms and entropy can be expressed as in (3.21), resulting the Hartree free energy

$$\begin{aligned} \mathcal {F}^\mathrm{H}[\gamma ]&:={{\,\mathrm{\mathrm{Tr}}\,}}\left[ (-\Delta +V-\nu )\gamma \right] +\frac{\lambda }{2} \iint \gamma (x;x) w(x-y) \gamma (y;y)\, dx\, dy\\&\quad + \frac{\lambda }{2} \iint w(x-y) |\gamma (x;y)|^2 dx dy\\&\quad -T{{\,\mathrm{\mathrm{Tr}}\,}}\left[ (1+\gamma )\log (1+\gamma )-\gamma \log \gamma \right] \end{aligned}$$

Thus if we are interested in equilibrium states minimizing the free energy, in the quasi-free class this leads to the following variational problem

$$\begin{aligned} F^\mathrm{H}_\lambda =\inf _{\gamma =\gamma ^*\geqslant 0}\mathcal {F}^\mathrm{H}[\gamma ]. \end{aligned}$$

(A.2)

When $\widehat{w}\geqslant 0$, the functional $\mathcal {F}^\mathrm{H}[\gamma ]$ turns out to be strictly convex. Hence, with the confining potential V it admits a unique minimizer $\gamma ^\mathrm{H}$, that defines a unique corresponding quasi-free state in Fock space $\Gamma ^\mathrm{H}$ (the proof is the same as that of Lemma 3.2). The optimal density matrix solves the nonlinear equation

$$\begin{aligned} \gamma ^\mathrm{H}=\left\{ \exp \left( \frac{-\Delta +V-\nu +\lambda \rho ^\mathrm{H}*w+\lambda X^\mathrm{H}}{T}\right) -1\right\} ^{-1} \end{aligned}$$

where $\rho ^\mathrm{H}(x)=\gamma ^\mathrm{H}(x;x)$ is the density and $X^\mathrm{H}$ is the exchange operator with integral kernel $X^\mathrm{H}(x;y)=w(x-y)\gamma ^\mathrm{H}(x;y)$.

In the limit $\lambda \rightarrow 0$ with $T=1/\lambda$, the quasi-free state $\Gamma ^\mathrm{H}$ is rather badly behaved. Its density $\rho ^\mathrm{H}$ diverges very fast. However, it turns out that, although $\rho ^\mathrm{H}(x)$ depends on x, its growth as $\lambda \rightarrow 0$ is more or less uniform in x and can be captured by

$$\begin{aligned} \rho ^\mathrm{H}(x) \sim \varrho _0^{\kappa }(\lambda ) := \left[ \frac{1}{e^{\lambda (-\Delta +\kappa )}-1} \right] (x;x) = \frac{1}{(2\pi )^d\lambda ^{\frac{d}{2}}} \int _{k\in {\mathbb {R} }^d} \frac{dk}{e^{|k|^2+\lambda \kappa } -1} \end{aligned}$$

(A.3)

provided that

$$\begin{aligned} \nu =\lambda \widehat{w} (0) \varrho _0^\kappa -\kappa . \end{aligned}$$

(A.4)

Recall that $\lambda \varrho _0^\kappa (\lambda )$ diverges in dimension $d\geqslant 2$ but it does not depend on x by translation invariance of $-\Delta + \kappa$. On the other hand, the exchange term $\lambda X^\mathrm{H}$ typically stays bounded, for instance in the Hilbert-Schmidt norm.

This suggests to simplify things a little bit by removing the exchange term as we did in the paper, that is, to consider the simplified minimization problem associated with the reduced Hartree or, simply, mean-field free energy (3.22) which, we recall, is given by

$$\begin{aligned} \mathcal {F}^\mathrm{MF}[\gamma ]&:={{\,\mathrm{\mathrm{Tr}}\,}}\left[ (-\Delta +V-\nu )\gamma \right] +\frac{\lambda }{2} \iint \gamma (x;x) w(x-y) \gamma (y;y)\, dx\, dy\\&\quad -T{{\,\mathrm{\mathrm{Tr}}\,}}\left[ (1+\gamma )\log (1+\gamma )-\gamma \log \gamma \right] \end{aligned}$$

By doing so we will pick as reference state a quasi-free state which is not the absolute minimizer of the true quantum free energy in the quasi-free class. However, manipulating states depending only on a potential simplifies the analysis. Fortunately, it turns out to be sufficient for our purpose, as justified in Lemma 3.2, which we now prove.

1.2 A.2. Proof of Lemma 3.2

Under our assumptions, the first eigenvalue of the Friedrichs realization of $h=-\Delta +V$ is positive. In addition, we have

$$\begin{aligned} {{\,\mathrm{\mathrm{Tr}}\,}}[e^{-h/T}]\leqslant (2\pi )^{-d}\int _{{\mathbb {R} }^d}e^{-p^2/T}\,dp\int _{{\mathbb {R} }^d}e^{-V(x)/T}\,dx<\infty , \end{aligned}$$

by the Golden–Thompson inequality, see [150, Section 8.1]. The same properties hold if we shift V by $\nu \in {\mathbb {R} }$ and only keep the positive part. Then we obtain

$$\begin{aligned}&{{\,\mathrm{\mathrm{Tr}}\,}}\left( \left( -\Delta +(V-\nu )_+\right) \gamma \right) -T{{\,\mathrm{\mathrm{Tr}}\,}}\left( (1+\gamma )\log (1+\gamma )-\gamma \log \gamma \right) \nonumber \\&\quad \geqslant \frac{1}{2}{{\,\mathrm{\mathrm{Tr}}\,}}\left( -\Delta \gamma \right) + T{{\,\mathrm{\mathrm{Tr}}\,}}\left( \log \left( 1-e^{-\frac{-\Delta /2+(V-\nu )_+}{T}}\right) \right) \end{aligned}$$

(A.5)

for all $\gamma \geqslant 0$. The last expression is the minimum free energy associated with $-\Delta /2+(V-\nu )_+$.

In order to prove that the functional $\mathcal {F}^\mathrm{MF}$ is bounded from below, it therefore remains to show that

$$\begin{aligned} -\int _{{\mathbb {R} }^d}\rho (x)(V-\nu )_-(x)\,dx+\frac{\lambda }{2} \int _{{\mathbb {R} }^d}\int _{{\mathbb {R} }^d}w(x-y)\rho (x)\rho (y)\,dx\,dy \end{aligned}$$

is bounded from below, uniformly in $\rho \geqslant 0$, under the assumption that $\widehat{w}\geqslant 0$ and $w\not \equiv 0$ (if $w\equiv 0$ then the Lemma holds true trivially). Since $w\in L^1({\mathbb {R} }^d)$, we can find a $k_0\in {\mathbb {R} }^d$ and a small radius $r>0$ such that the continuous non-negative function $\widehat{w}$ is at least equal to $\widehat{w}(k_0)/2$ on $B(k_0,r)$. We then choose $\varphi$ in the Schwartz class such that $\varphi >0$ and $\widehat{\varphi }\geqslant 0$ with $\mathrm{supp}(\widehat{\varphi })\subset B(k_0,r)$. Since $V\geqslant 0$ and $V\rightarrow +\infty$ at infinity, the function $(V-\nu )_-$ has compact support and is bounded by $\nu _+$. Therefore, we have

$$\begin{aligned} (V-\nu )_-\leqslant C\varphi \end{aligned}$$

for $C= \left| \left| \varphi ^{-1}(V-\nu )_- \right| \right| _{L^\infty ({\mathbb {R} }^d)}<\infty$. After completing the square and using $\widehat{w}\geqslant 0$, we then deduce

$$\begin{aligned}&-\int _{{\mathbb {R} }^d}(V-\nu )_-\rho +\frac{\lambda }{2} \int _{{\mathbb {R} }^d}\int _{{\mathbb {R} }^d}w(x-y)\rho (x)\rho (y)\,dx\,dy\\&\quad \geqslant -C\int _{{\mathbb {R} }^d}\varphi \rho +\frac{\lambda }{2} \int _{{\mathbb {R} }^d}\int _{{\mathbb {R} }^d}w(x-y)\rho (x)\rho (y)\,dx\,dy \\&\quad \geqslant -\frac{C^2}{2\lambda } \int _{{\mathbb {R} }^d}\frac{|\widehat{\varphi }(k)|^2}{\widehat{w}(k)}\,dk\\&\quad \geqslant -\frac{C^2\Vert \widehat{\varphi }\Vert ^2_{L^2}}{\lambda \widehat{w}(k_0)}. \end{aligned}$$

Combining with (A.5) we find as stated that

$$\begin{aligned} \inf _{\gamma =\gamma ^*\geqslant 0}\mathcal {F}^\mathrm{MF}[\gamma ]>-\infty , \end{aligned}$$

for all $\nu \in {\mathbb {R} }$.

Let us now prove the existence of a minimizer. Writing

$$\begin{aligned} \mathcal {F}^\mathrm{MF}_\nu [\gamma ]=\mathcal {F}^\mathrm{MF}_{\nu +1}[\gamma ]+{{\,\mathrm{\mathrm{Tr}}\,}}\gamma \geqslant \inf _{\gamma '} \mathcal {F}^\mathrm{MF}_{\nu +1}[\gamma ']+{{\,\mathrm{\mathrm{Tr}}\,}}(\gamma ), \end{aligned}$$

where we have displayed the parameter $\nu$ for convenience, we obtain that minimizing sequences $\{\gamma _n\}$ for $\mathcal {F}^\mathrm{MF}_\nu$ are necessarily bounded in the trace-class. In particular, $\Vert \gamma _n\Vert$ is also bounded. In addition, the inequality (A.5) implies that ${{\,\mathrm{\mathrm{Tr}}\,}}(-\Delta )\gamma _n$ is bounded. From the Hoffmann-Ostenhof inequality [88]

$$\begin{aligned} {{\,\mathrm{\mathrm{Tr}}\,}}(-\Delta )\gamma \geqslant \int _{{\mathbb {R} }^d}\left| \nabla \sqrt{\rho _\gamma (x)}\right| ^2\,dx \end{aligned}$$

and the Sobolev inequality, we deduce that $\rho _{\gamma _n}$ is bounded in $L^{p^*/2}({\mathbb {R} }^d)$ where $p^*=+\infty$ in dimension $d=1$, $p^*<\infty$ arbitrarily in dimension $d=2$ and $p^*=2d/(d-2)$ in dimensions $d\geqslant 3$. Hence, up to extraction of a subsequence, we have $\gamma _n\rightharpoonup \gamma$ weakly-$*$ in $\mathfrak {S}^1$ and $\rho _{\gamma _n}\rightharpoonup \rho _{\gamma }$ weakly in $L^1({\mathbb {R} }^d)\cap L^{p^*/2}({\mathbb {R} }^d)$. Using Fatou’s lemma and the concavity (hence weak upper semi-continuity) of the entropy, we obtain that $\gamma$ is a minimizer for $\mathcal {F}_\nu ^\mathrm{MF}$. The nonlinear equation (3.24) follows from classical arguments. Then, according to (A.1), $V_\lambda$ must solve (3.20) because the minimizer $\gamma ^\mathrm{MF}$ is the one-body density matrix of the quasi-free state associated with $\mathrm {d}\Gamma (-\Delta + V_\lambda )$. $\square$

1.3 A.3. Comments on Theorem 3.3

Let us briefly discuss Theorem 3.3. In [60, Section 5], the existence of the solution $V_\lambda$ to (3.20) was established by means of a fixed-point argument (which requires that $d\leqslant 3$ and that $\kappa$ is sufficiently large). The fixed point is performed in the (complete) metric space

$$\begin{aligned} B(V)=\left\{ f \in L^\infty _\mathrm{loc} ({\mathbb {R} }^d): \Vert f\Vert _{B(V)}= \left\| \frac{f}{V} - 1 \right\| _{L^\infty } \leqslant 1/2 \right\} \end{aligned}$$

for the unknown $u=V_\lambda -\kappa$ and provides the Hilbert-Schmidt convergence

$$\begin{aligned} {{\,\mathrm{\mathrm{Tr}}\,}}\Big | (-\Delta +V_\lambda )^{-1} -(-\Delta +V_0)^{-1} \Big |^2 \rightarrow 0. \end{aligned}$$

Our notation here is slightly different from [60] as we shift potentials by a constant. Moreover, since $V_\lambda -\kappa \in B(V)$ we have

$$\begin{aligned} \frac{V}{2} \leqslant V_\lambda -\kappa \leqslant 3\frac{V}{2}. \end{aligned}$$

There remains to discuss the nonlinear equation (3.28) for $V_0$, which we can rewrite in the form

$$\begin{aligned} V_0=w*\left( V+\kappa +\rho \left[ \big (-\Delta +V_0\big )^{-1}-\big (-\Delta +\kappa \big )^{-1}\right] \right) . \end{aligned}$$

(A.6)

Here we just need to pass to the limit in the similar equation at $\lambda >0$

$$\begin{aligned} V_\lambda =w*\left( V+\kappa +\rho \left[ \frac{\lambda }{ e^{\lambda (-\Delta +V_\lambda )}-1}-\frac{\lambda }{ e^{\lambda (-\Delta +\kappa )}-1}\right] \right) . \end{aligned}$$

Since we know that $V_\lambda /V\rightarrow V_0/V$ in $L^\infty$, we have $V_\lambda \rightarrow V_0$ in $L^\infty _\mathrm{loc}$ and it suffices to prove the convergence of the density on the right side, which we denote for simplicity

$$\begin{aligned} \rho _\lambda ^{V_\lambda }(x):=\left[ \frac{\lambda }{ e^{\lambda (-\Delta +V_\lambda )}-1}-\frac{\lambda }{e^{\lambda (-\Delta +\kappa )}-1}\right] (x;x). \end{aligned}$$

In [60, Eq. (5.21)] it is shown that

$$\begin{aligned} |\rho _\lambda ^{V_\lambda }(x)|\leqslant C\kappa ^{d/2-2} \left| \left| \frac{V_\lambda - \kappa }{V} \right| \right| _{L^\infty } V(x). \end{aligned}$$

(A.7)

Hence from the dominated convergence theorem and the assumptions on w, it suffices to prove that

$$\begin{aligned} \rho _\lambda ^{V_\lambda }(x)\rightarrow \left( \big (-\Delta +V_0\big )^{-1}-\big (-\Delta +\kappa \big )^{-1}\right) (x;x) \end{aligned}$$

almost everywhere. Applying again [60, Eq. (5.21)] we find that

$$\begin{aligned} \left| \left( \frac{\lambda }{e^{\lambda (-\Delta +V_\lambda )}-1}-\frac{\lambda }{ e^{\lambda (-\Delta +V_0)}-1}\right) (x;x)\right| \leqslant C\kappa ^{d/2-2} \left| \left| \frac{V_\lambda -V_0}{V} \right| \right| _{L^\infty } V(x) \end{aligned}$$

which tends to 0 in $L^\infty _\mathrm{loc}$ since $(V_\lambda -V_0)/V\rightarrow 0$. Hence we can replace $V_\lambda$ by $V_0$ throughout. Next we write, following [60, Eq. (5.16)],

$$\begin{aligned} \rho _\lambda ^{V_0}(x)=-\int _0^1\mathrm{d}s\int _{{\mathbb {R} }^d}\mathrm{d}z\frac{\lambda \,e^{s\lambda (-\Delta +V_0)}}{e^{\lambda (-\Delta +V_0)}-1}(x;z)V_0(z)\frac{\lambda \,e^{(1-s)\lambda (-\Delta +\kappa )}}{e^{\lambda (-\Delta +\kappa )}-1}(z;x). \end{aligned}$$

Using that $V_0\geqslant \kappa$, we have the pointwise bound on the operator kernels

$$\begin{aligned} 0\leqslant \frac{e^{s\lambda (-\Delta +V_0)}}{e^{\lambda (-\Delta +V_0)}-1}(x;z)\leqslant \frac{e^{s\lambda (-\Delta +\kappa )}}{e^{\lambda (-\Delta +\kappa )}-1}(x;z) \end{aligned}$$

by the same argument as in Lemma 11.2 and in [60, Eq. (5.17)]. Using [60, Lemma 5.4] we see that we get a convergent domination. So by the dominated convergence theorem, the strong local convergence of $\rho ^{V_0}_\lambda$ follows from that of the kernels

$$\begin{aligned} \frac{\lambda \, e^{\lambda (-\Delta +V_0)}}{e^{\lambda (-\Delta +V_0)}-1}(x;z)\rightarrow \frac{1}{-\Delta +V_0}(x;z), \end{aligned}$$

$$\begin{aligned} \frac{\lambda \,e^{s\lambda (-\Delta +\kappa )}}{e^{\lambda (-\Delta +\kappa )}-1}(x;z)\rightarrow \frac{1}{-\Delta +\kappa }(x;z). \end{aligned}$$

In fact this convergence is strong in $L^2({\mathbb {R} }^d\times {\mathbb {R} }^d)$ since the corresponding operators converge in the Hilbert-Schmidt norm, by [60, Lemma C.1]. Passing to the limit, this proves the strong local convergence

$$\begin{aligned}&\rho ^{V_0}_\lambda (x)\underset{\lambda \rightarrow 0^+}{\longrightarrow } -\int _0^1\mathrm{d}s\int _{{\mathbb {R} }^d}\mathrm{d}z\frac{1}{-\Delta +V_0}(x;z)V_0(z)\frac{1}{-\Delta +\kappa }(z;x)\\&\quad =\left( \big (-\Delta +V_0\big )^{-1}-\big (-\Delta +\kappa \big )^{-1}\right) (x;x), \end{aligned}$$

where in the last equality we have used the resolvent formula. The uniform bound (A.7) then allows to pass to the limit in the equation for $V_\lambda$ and obtain (A.6). $\square$

Appendix B: Interpretation in terms of the phase transition of the Bose gas

Here we reformulate and discuss our main results (Theorems 3.1 and 3.4) in microscopic variables and clarify the link with the phase transition of the infinite Bose gas.

1.1 B.1. Homogeneous case

We start with the homogeneous case which is the usual setting in which the thermodynamics of the free Bose gas is formulated, see for instance [30, Section 5.2.5], [161, Sections 2.5.19–20] and [165, Chapter 4].

Let us consider the non-interacting Bose gas on the large torus $L\mathbb {T}^d$ of side length L. In the grand canonical setting we choose a chemical potential $\widetilde{\nu }<0$ and a temperature $T>0$. Our system is then represented by the Gaussian quantum state in Fock space, associated with the one-particle operator $(-\Delta _L-\widetilde{\nu })/T$ where $-\Delta _L$ is the L-periodic Laplacian. Its one-body density matrix is³

$$\begin{aligned} \widetilde{\Gamma }^{(1)}_L=\frac{1}{e^{\frac{-\Delta _L-\widetilde{\nu }}{T}}-1} \end{aligned}$$

and the number of particles per unit volume is given by

$$\begin{aligned} \frac{1}{L^d}\sum _{k\in 2\pi {\mathbb {Z} }^d/L}\frac{1}{e^{\frac{|k|^2-\widetilde{\nu }}{T}}-1}\underset{L\rightarrow \infty }{\longrightarrow }\frac{T^{\frac{d}{2}}}{(2\pi )^d}\int _{{\mathbb {R} }^d}\frac{dk}{e^{k^2-\widetilde{\nu }/T}-1}. \end{aligned}$$

The critical density for Bose–Einstein condensation is obtained in the limit $\widetilde{\nu }/T\rightarrow 0^-$ and it equals

$$\begin{aligned} \rho _c(T)=\frac{T^{\frac{d}{2}}}{(2\pi )^d} \int _{{\mathbb {R} }^d}\frac{dk}{e^{k^2}-1} ={\left\{ \begin{array}{ll} +\infty &{}\quad \text {for } d=1,2,\\ \frac{T^{\frac{d}{2}}\zeta \left( \frac{d}{2}\right) }{2^d\pi ^{\frac{d}{2}}}<\infty &{}\quad \text {for } d\geqslant 3. \end{array}\right. } \end{aligned}$$

(B.1)

The grand-canonical model in infinite space stops to exist at $\widetilde{\nu }=0$, where the one-body density matrix converges weakly to

$$\begin{aligned} \widetilde{\Gamma }_\infty ^{(1)}=\frac{1}{e^{\frac{-\Delta }{T}}-1}. \end{aligned}$$

(B.2)

The corresponding infinite Gaussian state (properly defined over the $C^*$-algebra of Canonical Commutation Relations [30]) has the number of particles per unit volume equal to $\rho _c(T)$.

The phenomenon of Bose–Einstein Condensation (BEC) is better understood in the canonical setting with N particles, going back to the thermodynamic limit. In dimensions $d\geqslant 3$, fixing the density $N/L^d=\rho >\rho _c(T)$ one obtains a density matrix which has a two-scale behavior [30, Sec. 5.2.5]. It contains a rank-one part with the diverging eigenvalue $L^d(\rho -\rho _c(T))$ and constant eigenfunction $f_0(x)=L^{-d/2}$, living at the macroscopic scale (the Bose–Einstein condensate). When this rank-one operator is removed, the rest of the density matrix converges weakly to $\widetilde{\Gamma }_\infty ^{(1)}$ in (B.2).

To understand the behavior of the system just before the phase transition, we have to look at the simultaneous limit $L\rightarrow \infty$ with $\widetilde{\nu }\rightarrow 0^-$, see [161, Sec. 2.5.20.3]. At the macroscopic scale (that is, after rescaling length by L), the (rescaled) one-particle density matrix equals

$$\begin{aligned} \Gamma ^{(1)}_L=\frac{1}{e^{\frac{-\Delta -L^2\widetilde{\nu }}{L^2T}}-1}\quad \text {in } L^2(\mathbb {T}^d) \end{aligned}$$

and thus we see that the natural scaling for $\widetilde{\nu }$ is

$$\begin{aligned} \widetilde{\nu }(L)=-\frac{\kappa }{L^2} \end{aligned}$$

(B.3)

with a fixed $\kappa >0$, in all dimensions $d\geqslant 1$. We are therefore exactly in the setting studied in this paper and in [104] with the choice

$$\begin{aligned} \boxed {\lambda =\frac{1}{TL^2}\underset{L\rightarrow \infty }{\longrightarrow }0^+.} \end{aligned}$$

The density matrix converges to

$$\begin{aligned} \frac{\Gamma _L^{(1)}}{TL^2}\underset{L\rightarrow \infty }{\longrightarrow }\frac{1}{-\Delta +\kappa }=\int |u\rangle \langle u|\,d\mu _\kappa (u) \end{aligned}$$

(B.4)

strongly in the Schatten space $\mathfrak {S}^p(L^2(\mathbb {T}^d))$ for all $p>d/2$ ($p\geqslant 1$ if $d=1$), where $\mu _\kappa$ is the classical Gaussian measure with covariance $(-\Delta +\kappa )^{-1}$. Equivalently,

$$\begin{aligned} \frac{\Gamma _L^{(1)}}{L^2}\underset{L\rightarrow \infty }{\longrightarrow }\frac{T}{-\Delta +\kappa }=\int |u\rangle \langle u|\,d\mu _{\kappa ,T}(u) \end{aligned}$$

where $\mu _{\kappa ,T}$ is the Gaussian measure with covariance $T(-\Delta +\kappa )^{-1}$. Similar properties hold for higher density matrices.

Our conclusion is that, close to its phase transition, the free Bose gas is properly described by the classical Gaussian measure $\mu _{\kappa ,T}$ on $\mathbb {T}^d$ (macroscopic scale) where $\kappa$ describes the speed at which the chemical potential $\widetilde{\nu }$ approaches 0 via (B.3) or, equivalently, at which the corresponding density approaches the critical density $\rho _c(T)$. This elementary fact is rarely mentioned in textbooks on Bose gases.

The speed at which the density approaches $\rho _c(T)$ is computed by using the following

Lemma B.1

(Particle number of the free Bose gas) In dimensions $d\geqslant 1$, we have for every fixed $\kappa >0$

$$\begin{aligned} \sum _{k\in 2\pi {\mathbb {Z} }^d}\frac{1}{e^{\lambda (|k|^2+\kappa )}-1}=\frac{\lambda ^{-\frac{d}{2}}}{(2\pi )^d}\int _{{\mathbb {R} }^d}\frac{dk}{e^{k^2+\lambda \kappa }-1}+\frac{\varphi _d(\kappa )}{\lambda }+ o\left( \lambda ^{-1}\right) _{\lambda \rightarrow 0^+} \end{aligned}$$

(B.5)

with the positive decreasing function

$$\begin{aligned} \varphi _d(\kappa )=\sum _{\ell \in {\mathbb {Z} }^d{\setminus }\{0\}}\int _0^\infty \frac{e^{-t\kappa }}{(4\pi t)^{\frac{d}{2}}} e^{-\frac{|\ell |^2}{4t}}\,dt ={\left\{ \begin{array}{ll} \displaystyle \frac{1}{2}\sum _{\ell \in {\mathbb {Z} }{\setminus }\{0\}} e^{-\kappa |\ell |}&{}\quad \text {for } d=1,\\ \displaystyle \frac{1}{4\pi }\sum _{\ell \in {\mathbb {Z} }^3{\setminus }\{0\}}\frac{e^{-\sqrt{\kappa }|\ell |}}{|\ell |}&{}\quad \text {for } d=3. \end{array}\right. } \end{aligned}$$

(B.6)

The proof is provided below in “Appendix B.3”. The integral in the first equality of (B.6) is the Fourier transform of $(2\pi )^{-d/2}(|k|^2+\kappa )^{-1}$, that is, the Klein-Gordon Green function. In dimensions $d\leqslant 3$, we can also write by Poisson’s formula

$$\begin{aligned} \varphi _d(\kappa )=\sum _{k\in 2\pi {\mathbb {Z} }^d}\left( \frac{1}{|k|^2+\kappa }-\frac{1}{(2\pi )^d}\int _{(-\pi ,\pi )^d}\frac{dp}{|k+p|^2+\kappa }\right) \end{aligned}$$

but the sum is not absolutely convergent in dimensions $d\geqslant 4$. Since we have the expansions

$$\begin{aligned} \frac{1}{(2\pi )^d}\int _{{\mathbb {R} }^d}\frac{dk}{e^{|k|^2+a}-1}&={\left\{ \begin{array}{ll} \frac{1}{2\sqrt{a}}+\frac{\zeta (\frac{1}{2})}{2\sqrt{\pi }}+o(1)_{a\rightarrow 0^+}&{}\quad \text {for } d=1,\\ \frac{-\log (a)}{4\pi }+\frac{a}{8\pi }+o(a)_{a\rightarrow 0^+}&{}\quad \text {for } d=2, \end{array}\right. }\nonumber \\&=\frac{\zeta \left( \frac{d}{2}\right) }{2^d\pi ^{\frac{d}{2}}}-{\left\{ \begin{array}{ll} \frac{\sqrt{a}}{4\pi }+O(a)_{a\rightarrow 0^+}&{}\quad \text {for } d=3,\\ \frac{(\log (1/a)+1)a}{16\pi ^2}+o(a)_{a\rightarrow 0^+}&{}\quad \text {for } d=4,\\ \frac{\zeta \left( \frac{d}{2}-1\right) }{2^d\pi ^{\frac{d}{2}}}a+o(a)_{a\rightarrow 0^+}&{}\quad \text {for } d\geqslant 5. \end{array}\right. } \end{aligned}$$

(B.7)

we find from (B.5) that the density approaches the critical density $\rho _c(T)$ from below as

$$\begin{aligned}&\frac{{{\,\mathrm{\mathrm{Tr}}\,}}[\Gamma ^{(1)}_L]}{L^d} ={\left\{ \begin{array}{ll} \left( \frac{1}{2\sqrt{\kappa }}+\varphi _1(\kappa )\right) TL+o(L)&{}\quad \text {for } d=1,\\ \frac{T}{2\pi }\log (L)+\left( \varphi _2(\kappa )-\frac{\log (\kappa /T)}{4\pi }\right) T+o(1)&{}\quad \text {for } d=2, \end{array}\right. }\nonumber \\&\quad =\frac{T^{\frac{d}{2}}\zeta \left( \frac{d}{2}\right) }{2^d\pi ^{\frac{d}{2}}}-{\left\{ \begin{array}{ll} \left( \frac{\sqrt{\kappa }}{4\pi }-\varphi _3(\kappa )\right) \frac{T}{L}+o(L^{-1})&{}\quad \text {for } d=3,\\ \frac{\kappa T}{8\pi ^2}\frac{\log L}{L^2}+\left( \frac{1+\log \frac{T}{\kappa }}{16\pi ^2}\kappa -\varphi _4(\kappa )\right) \frac{T}{L^2}+o(L^{-2})&{}\quad \text {for } d=4,\\ \frac{\zeta \left( \frac{d}{2}-1\right) }{2^d\pi ^{\frac{d}{2}}}\frac{T^{\frac{d}{2}-1}\kappa }{L^2}+o(L^{-2})&{}\quad \text {for } d\geqslant 5. \end{array}\right. } \end{aligned}$$

(B.8)

Note that $\varphi _3(\kappa )\leqslant \sqrt{\kappa }/(4\pi )$ which is its behavior when $\kappa \rightarrow 0^+$.

So far our discussion applies to any dimension $d\geqslant 1$. Next we discuss the inclusion of interactions for $d\leqslant 3$. In our work the potential w is introduced at the macroscopic level. Re-expressed in microscopic variables and taking $T=1$ for simplicity, we obtain the microscopic n-particle Hamiltonian

$$\begin{aligned} \widetilde{H}_{n,L}=\sum _{j=1}^n(-\Delta _L-\widetilde{\nu })_{x_j}+\frac{1}{L^4}\sum _{1\leqslant j<k\leqslant n}w\left( \frac{x_j-x_k}{L}\right) . \end{aligned}$$

(B.9)

The interaction has the very small intensity $L^{-4}$ but varies on length scales comparable with the size of the box. The form of the microscopic Hamiltonian (B.9) is the same in all dimensions. The thermodynamic limit $L\rightarrow \infty$ of this model at fixed $\widetilde{\nu }<0$ is the same as the non-interacting case. This is because we have the lower bound

$$\begin{aligned} \widetilde{H}_{n,L}\geqslant \sum _{j=1}^n-(\Delta _L)_{x_j}-\left( \widetilde{\nu }+\frac{w(0)}{2L^4}\right) n \end{aligned}$$

due to the fact that $\widehat{w}\geqslant 0$ (for an upper bound on the free energy, use the non-interacting state). We conclude that w does not, to leading order, change the phase diagram as compared to the non-interacting case. The effect of w is only visible when zooming just before the phase transition. From Theorem 3.1, when the chemical potential goes to zero as

$$\begin{aligned} \widetilde{\nu }(L)=\frac{\nu (L^{-2})}{L^2}={\left\{ \begin{array}{ll} \displaystyle \frac{\widehat{w} (0)\,\log (L)}{2\pi L^2} -\frac{\nu _0}{L^2}+o(1)_{\lambda \rightarrow 0^+}&{}\quad \text {for } d=2,\\ \displaystyle \frac{\widehat{w} (0)\zeta \left( \frac{3}{2}\right) }{8\pi ^{\frac{3}{2}}L}-\frac{\nu _0}{L^2}+o(1)_{\lambda \rightarrow 0^+}&{}\quad \text {for } d=3, \end{array}\right. } \end{aligned}$$

then the behavior close to the transition is described by the nonlinear Gibbs measure $\mu$ at the macroscopic scale, which depends on w and $\kappa$ solving

$$\begin{aligned} \nu _0={\left\{ \begin{array}{ll} \kappa +\widehat{w}(0)\frac{\log (\kappa )}{4\pi }-\widehat{w}(0)\,\varphi _2(\kappa )&{}\quad \text {for } d=2,\\ \kappa +\widehat{w}(0)\frac{\sqrt{\kappa }}{4\pi }-\widehat{w}(0)\,\varphi _3(\kappa )&{}\quad \text {for } d=3, \end{array}\right. } \end{aligned}$$

(B.10)

More physical interactions are much bigger and have a much shorter range. Although a universal behavior can still be expected at the phase transition, the phase diagram depends on w at leading order and a mathematical treatment seems out of reach with the present techniques. A simpler behavior is however expected in the dilute regime $\rho \rightarrow 0$ with $\rho \sim T^{d/2}$ (Gross-Pitaevskii regime [48]). In dimension $d=3$ and at our macroscopic scale, the Gross-Pitaevskii limit corresponds to replacing $\lambda w$ by $\lambda w_\lambda$ with $w_\lambda (x)=\lambda ^{-3} w(x/\lambda )$ in our many-particle Hamiltonian. In this case one would expect the phase transition to be described by the (appropriately renormalized) nonlinear Gibbs measure $\mu$ over the torus $\mathbb {T}^3$, with w replaced by the Dirac delta $8\pi a \delta _0$ where a is the scattering length of w [48, 115, 116]. Proving such a result seems a formidable task.

1.2 B.2. Trapped gases

The theory of Bose–Einstein condensation for trapped gases is analogous to the homogeneous case, but the formulas are slightly different, see for instance [9, 161, Sec. 2.5.15] and [49]. Here we only discuss the case $V(x)=|x|^s$ for simplicity. At the microscopic scale the one-body Hamiltonian takes the form

$$\begin{aligned} \widetilde{h}_L:=-\Delta +\frac{|x|^s}{L^{2+s}}-\widetilde{\nu } \end{aligned}$$

where L is now a parameter used to open the trap whenever the number of particles grows. At fixed $\widetilde{\nu }<0$, the number of particles in the non-interacting Gaussian state is given by

$$\begin{aligned}&{{\,\mathrm{\mathrm{Tr}}\,}}\left( \frac{1}{e^{T^{-1}(-\Delta +L^{-2-s}|x|^s-\widetilde{\nu })}-1}\right) \nonumber \\&\quad ={{\,\mathrm{\mathrm{Tr}}\,}}\left( \frac{1}{e^{-T^{-1}\left( L^{-\frac{4+2s}{s}}\Delta +|x|^s-\widetilde{\nu }\right) }-1}\right) \nonumber \\&\quad \underset{L\rightarrow \infty }{\sim } \frac{(L^2T)^{d\left( \frac{1}{2}+\frac{1}{s}\right) }}{(2\pi )^d}\iint _{{\mathbb {R} }^d\times {\mathbb {R} }^d}\frac{dx\,dk}{e^{|k|^2+|x|^s-\widetilde{\nu }/T}-1}. \end{aligned}$$

(B.11)

In the first equality we have rescaled lengths by the factor $L^{(2+s)/s}$, which places the system in a conventional semi-classical limit with effective parameter $\hbar =L^{-(2+s)/s}\rightarrow 0$, hence the second limit. Computing in the same manner the average against |x| one sees that the gas is extended at the length scale

$$\begin{aligned} \ell _\mathrm{gas}\sim (L^2T)^{\left( \frac{1}{2}+\frac{1}{s}\right) }T^{-\frac{1}{2}}. \end{aligned}$$

(B.12)

Dividing by the effective volume $(\ell _\mathrm{gas})^d$, the density obtained in the thermodynamic limit is therefore proportional to

$$\begin{aligned} \frac{T^{\frac{d}{2}}}{(2\pi )^d}\iint _{{\mathbb {R} }^d\times {\mathbb {R} }^d}\frac{dx\,dk}{e^{|k|^2+|x|^s-\widetilde{\nu }/T}-1}. \end{aligned}$$

Bose–Einstein condensation is obtained as before when $\widetilde{\nu }/T\rightarrow 0^-$, with the critical density

$$\begin{aligned} \rho '_c(T)=\frac{T^{\frac{d}{2}}}{(2\pi )^d}\iint _{{\mathbb {R} }^d\times {\mathbb {R} }^d}\frac{dx\,dk}{e^{|k|^2+|x|^s}-1}. \end{aligned}$$

This is finite for all $s>0$ in dimensions $d\geqslant 2$. In dimension $d=1$, $\rho '_c(T)$ is finite for $s<2$ and infinite otherwise.

As before the Bose–Einstein condensate emerges in the limit $L\rightarrow \infty$ in the canonical setting, when $N>\rho '_c(T)(\ell _\mathrm{gas})^d$. The corresponding condensate wavefunction is the first eigenfunction of $-\Delta +L^{-2-s}|x|^s$. This is nothing but that of $-\Delta +|x|^s$ dilated to the scale L. Therefore, in this system the BEC length scale is L and it is always smaller than the natural extension length $\ell _\mathrm{gas}$ of the cloud in (B.12), at a fixed temperature $T>0$. The two coincide only in a sharp container ($s=+\infty$).

The Gaussian Gibbs measure based on $h=-\Delta +|x|^s$ emerges at the BEC length scale L, whenever the chemical potential is chosen as

$$\begin{aligned} \widetilde{\nu }(L)=-\frac{\kappa }{L^2}. \end{aligned}$$

At this scale the one-particle Hamiltonian just becomes $(-\Delta +|x|^s+\kappa )/L^2$ so that $\lambda =1/(TL^2)$ like in the homogeneous case. The arguments from [60] and “Appendix A.3” in the non-interacting case give

$$\begin{aligned}&\rho \left[ \frac{1}{e^{\frac{-\Delta +|x|^s+\kappa }{TL^2}}-1}\right] (x)=\frac{T^{\frac{d}{2}}L^d}{(2\pi )^d}\int _{{\mathbb {R} }^d}\frac{dk}{e^{|k|^2+\frac{\kappa }{TL^2}}-1}\nonumber \\&\quad +TL^2\;\rho \left[ \frac{1}{-\Delta +|x|^s+\kappa }-\frac{1}{-\Delta +\kappa }\right] (x)+o(TL^2)_{L\rightarrow \infty } \end{aligned}$$

(B.13)

for every $x\in {\mathbb {R} }^d$. Here the density on the second line plays the role of $\varphi _d$ in the homogeneous case. In particular we obtain all the same formulas as in (B.8) with $\varphi _d(\kappa )$ replaced by

$$\begin{aligned} \varphi _d'(\kappa ):={{\,\mathrm{\mathrm{Tr}}\,}}\left[ \frac{1}{-\Delta +|x|^s+\kappa }-\frac{1}{-\Delta +\kappa }\right] . \end{aligned}$$

In the inhomogeneous case, interactions were introduced at the scale L of the BEC and the interpretation is the same as in the homogeneous case.

1.3 B.3. Proof of Lemma B.1

We write

$$\begin{aligned} \frac{\lambda }{e^{\lambda h}-1}=\lambda \sum _{n\geqslant 1}e^{-n\lambda h} \end{aligned}$$

and obtain

$$\begin{aligned}&\lambda \left( \sum _{k\in 2\pi {\mathbb {Z} }^d}\frac{1}{e^{\lambda (|k|^2+\kappa )}-1}-\frac{1}{(2\pi )^d}\int _{{\mathbb {R} }^d}\frac{dp}{e^{\lambda (|p|^2+\kappa )}-1}\right) \\&\quad =\lambda \sum _{n\geqslant 1}e^{-n\lambda \kappa }\sum _{k\in 2\pi {\mathbb {Z} }^d}\left( e^{-n\lambda |k|^2}-\frac{1}{(2\pi )^d}\int _{(-\pi ,\pi )^d}e^{-n\lambda |k-p|^2}\,dp\right) . \end{aligned}$$

We use Poisson’s formula

$$\begin{aligned} \sum _{k\in 2\pi {\mathbb {Z} }^d}\widehat{f}(k)=\frac{1}{(2\pi )^{\frac{d}{2}}}\sum _{\ell \in {\mathbb {Z} }^d}f(\ell ) \end{aligned}$$

for $f(x)=e^{-n\lambda |x|^2}-e^{-n\lambda |\cdot |^2}*\chi (x)$ where $\chi =(2\pi )^{-d}{\mathbb {1} }_{(-\pi ,\pi )^d}$ and find

$$\begin{aligned}&\lambda \left( \sum _{k\in 2\pi {\mathbb {Z} }^d}\frac{1}{e^{\lambda (|k|^2+\kappa )}-1}-\frac{1}{(2\pi )^d}\int _{{\mathbb {R} }^d}\frac{dp}{e^{\lambda (|p|^2+\kappa )}-1}\right) \\&\quad =\frac{\lambda }{(4\pi )^{\frac{d}{2}}}\sum _{n\geqslant 1}\frac{e^{-n\lambda \kappa }}{(4\pi n\lambda )^{\frac{d}{2}}}\sum _{\ell \in {\mathbb {Z} }^d{\setminus }\{0\}} e^{-\frac{|\ell |^2}{4n\lambda }}. \end{aligned}$$

This is a Riemann sum which converges to

$$\begin{aligned} \sum _{\ell \in {\mathbb {Z} }^d{\setminus }\{0\}}\int _0^\infty \frac{e^{-t\kappa }}{(4\pi t)^{\frac{d}{2}}} e^{-\frac{|\ell |^2}{4t}}\,dt \end{aligned}$$

where the right side is the Fourier transform of $k\mapsto (2\pi )^{-d/2}(|k|^2+\kappa )^{-1}$, see [114]. $\square$