Re-examining the Statistical Mechanics of an Interacting Bose Gas

TL;DR

This paper critically re-examines the statistical mechanics of interacting Bose gases, clarifies the canonical versus grand-canonical treatments, and proposes a new framework that allows the number of bosons to vary with temperature.

Contribution

It identifies the misinterpretation in previous treatments as a canonical approach and develops a canonical and grand-canonical formulation where boson number varies with temperature.

Findings

The grand-canonical treatment in previous theories is effectively canonical.

Discontinuity in condensate density at T_c is due to an inappropriate extension of Bogoliubov's prescription.

A new canonical and grand-canonical framework for interacting bosons is proposed.

Abstract

We re-examine the way in which Bogoliubov's theory of a dilute Bose gas at $T=0$ has been extended to describe the statistical mechanics of interacting bosons at finite temperature. We show explicitly that the field-theoretic calculation of the grand partition function in this formulation amounts to a canonical trace over the eigenfunctions of the Bogoliubov Hamiltonian at fixed total number of bosons $N$, and that the additional trace over $N$ that is required in the grand-canonical formalism is never carried out. This implies that what usually passes as the grand-canonical treatment of the Bogoliubov Hamiltonian is not quite grand-canonical, and is in fact a canonical one. We also show that the discontinuity in the condensate density predicted by previous formulations of this theory as the temperature $T$ goes past the critical transition temperature $T_c$ is a direct consequence of an inappropriate generalization of the Bogoliubov prescription to finite temperatures, and that this discontinuity disappears when this prescription is either used as a zero temperature approximation or avoided altogether. Armed with the above findings, we reformulate the statistical mechanics of interacting bosons in the canonical ensemble and derive the thermodynamics of the system. We then show how the canonical treatment can be used to setup a truly grand-canonical description of the statistical mechanics of a weakly interacting Bose gas where the average number of bosons in the system varies with temperature, unlike existing formulations where the total number of bosons $N$ is taken to be a constant that does not depend on $T$. Consequences on the physics of interacting bosons are briefly discussed.

Tables

Table 1. Table 1: Values of σ ~ ~ 𝜎 \tilde{\sigma} and C d subscript 𝐶 𝑑 C_{d} for a few representative values of the parameter n B ​ a 3 subscript 𝑛 𝐵 superscript 𝑎 3 n_{B}a^{3} .

$n_B{a}^3$	$\tilde{\sigma }$	$C_d$
${10}^{-5}$	0.614	0.9967
${10}^{-4}$	0.568	0.9897
${10}^{-3}$	0.452	0.9672
${10}^{-2}$	0.263	0.8990
${10}^{-1}$	0.172	0.7483

Table 2. Table 2: Numerical values of the coefficient C E subscript 𝐶 𝐸 C_{E} in the expression of the T = 0 𝑇 0 T=0 ground-state energy for a few representative values of the parameter n B ​ a 3 subscript 𝑛 𝐵 superscript 𝑎 3 n_{B}a^{3} as obtained from the canonical ensemble description of our variational number-conserving approach.

$n_B{a}^3$	$C_E$
${10}^{-4}$	57.62
${10}^{-3}$	14.64
${10}^{-2}$	2.37

Table 3. Table 3: Critical temperature T c / T c ​ 0 subscript 𝑇 𝑐 subscript 𝑇 𝑐 0 T_{c}/T_{c0} for a few representative values of the parameter n B ​ a 3 subscript 𝑛 𝐵 superscript 𝑎 3 n_{B}a^{3} as obtained from the canonical ensemble description of our variational number-conserving approach.

$n_B{a}^3$	$T_c/T_{c0}$
${10}^{-4}$	1.498
${10}^{-3}$	1.727
${10}^{-2}$	2.022

Equations

$$\Psi (\mathbf{r},\tau )\simeq \sqrt{n_0}+\psi (\mathbf{r},\tau ).$$

$$\frac{\Delta T_c}{T_{c0}}\propto (n_B{a}^3{)}^{\eta },$$

$$\eta =\frac{1}{6}.$$

$$Z_G=\int [d{\Psi }^{\ast }(\mathbf{r},\tau )][d\Psi (\mathbf{r},\tau )]e^{-S[{\Psi }^{\ast },\Psi ]/\hslash },$$

$$\displaystyle S=\int_{0}^{\beta\hbar}d\tau\Big{\{}\int d{\bf r}\,\Psi^{*}({\bf r},\tau)\Big{(}\hbar\partial_{\tau}-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu\Big{)}\Psi({\bf r},\tau)$$

$$\displaystyle+\frac{1}{2}\int d{\bf r}d{\bf r}^{\prime}\,\Psi^{*}({\bf r},\tau)\Psi^{*}({\bf r}^{\prime},\tau)V({\bf r}-{\bf r}^{\prime})\Psi({\bf r}^{\prime},\tau)\Psi({\bf r},\tau)\Big{\}}.$$

$$g=\frac{4\pi {\hslash }^{2}a}{m}.$$

$$\Psi (\mathbf{r},\tau )\simeq \sqrt{n_0}+\psi (\mathbf{r},\tau ),$$

$$\psi (\mathbf{r},\tau )=\sum _{\mathbf{k}\ne 0}\psi (\mathbf{k},\tau ){\phi }_{\mathbf{k}}(\mathbf{r}).$$

$$\psi (\mathbf{k},\tau )=\frac{1}{\sqrt{\beta \hslash }}\sum _{n=-\infty }^{\infty }\psi (\mathbf{k},{\omega }_n)e^{-i{\omega }_n\tau }$$

$$\psi (\mathbf{r},\tau )=\frac{1}{\sqrt{\beta \hslash V}}\sum _{\mathbf{k}\ne 0}\sum _{n=-\infty }^{\infty }\psi (\mathbf{k},{\omega }_n)e^{i(\mathbf{k}\cdot \mathbf{r}-{\omega }_n\tau )}.$$

$$S_0$$

$$\displaystyle+\sum_{\bf k\neq 0}\sum_{n=-\infty}^{\infty}\Big{\{}\psi^{*}({\bf k},\omega_{n})\big{(}-i\hbar\omega_{n}+\varepsilon_{\bf k}-\mu+2gn_{0}\big{)}\psi({\bf k},\omega_{n})$$

$$\displaystyle+\frac{gn_{0}}{2}\big{[}\psi({\bf k},\omega_{n})\psi(-{\bf k},-\omega_{n})+\psi^{*}({\bf k},\omega_{n})\psi^{*}(-{\bf k},-\omega_{n})\big{]}\Big{\}},$$

$$S_1$$

$$\displaystyle+|\psi({\bf r},\tau)|^{4}\Big{\}}.$$

$$S_0=\beta \hslash V(-\mu n_0+\frac{1}{2}g{n}_0^2)$$

$$\displaystyle+\sum_{\bf k\neq 0}\sum_{n=-\infty}^{\infty}\left\{\Big{[}\psi({\bf k},\omega_{n})\,\psi^{*}({\bf k},\omega_{n})\Big{]}G^{-1}({\bf k},\omega_{n})\left[\begin{array}[]{c}\psi({\bf k},\omega_{n})\\ \psi^{*}({\bf k},\omega_{n})\\ \end{array}\right]\right\},$$

$$\displaystyle G^{-1}({\bf k},\omega_{n})=\left(\begin{array}[]{cc}(G^{-1})_{11}&(G^{-1})_{12}\\ (G^{-1})_{21}&(G^{-1})_{22}\\ \end{array}\right),$$

$$(G^{-1}{)}_{11}$$

$$(G^{-1}{)}_{12}$$

$$(G^{-1}{)}_{22}$$

$$\displaystyle G({\bf k},\omega_{n})=\frac{1}{\hbar^{2}\omega_{n}^{2}+E_{\bf k}^{2}}\left(\begin{array}[]{cc}(G^{-1})_{22}&-(G^{-1})_{12}\\ -(G^{-1})_{21}&(G^{-1})_{11}\\ \end{array}\right),$$

$$E_{\bf k}=\sqrt{\big{(}\varepsilon_{\bf k}+2gn_{0}-\mu\big{)}^{2}-(gn_{0})^{2}}.$$

$$\mu =g{n}_0,$$

$$E_{\bf k}=\sqrt{\varepsilon_{\bf k}\big{(}\varepsilon_{\bf k}+2gn_{0}\big{)}}.$$

$$⟨\psi ({\mathbf{k}}_1,{\omega }_n){\psi }^{\ast }({\mathbf{k}}_2,{\omega }_m)⟩=\hslash G_{11}(\mathbf{k},{\omega }_n){\delta }_{{\mathbf{k}}_1,{\mathbf{k}}_2}{\delta }_{n,m},$$

$$⟨\psi ({\mathbf{k}}_1,{\omega }_n)\psi ({\mathbf{k}}_2,{\omega }_m)⟩=\hslash G_{12}(\mathbf{k},{\omega }_n){\delta }_{{\mathbf{k}}_1,{\mathbf{k}}_2}{\delta }_{n,m}.$$

$$N_{\mathbf{k}}$$

$$\displaystyle=\frac{\varepsilon_{\bf k}+gn_{0}}{2E_{\bf k}}\coth\Big{(}\frac{\beta E_{\bf k}}{2}\Big{)}-\frac{1}{2}.$$

$$n_d=\frac{1}{V}\sum _{\mathbf{k}\ne 0}{N}_{\mathbf{k}},$$

$$n_{d}(T)=\int\frac{d\bf k}{(2\pi)^{3}}\left[\frac{\varepsilon_{\bf k}+gn_{0}}{2E_{\bf k}}\coth\Big{(}\frac{\beta E_{\bf k}}{2}\Big{)}-\frac{1}{2}\right].$$

$$n_0(T)=n_B-n_d(T).$$

$$n_{d}(T)=\int_{\bf k}\left[\frac{\varepsilon_{\bf k}+gn_{0}}{2E_{\bf k}}\coth\Big{(}\frac{\beta E_{\bf k}}{2}\Big{)}-\frac{1}{2}\right].$$

$$\displaystyle n_{d}-n_{cr}=\int_{\bf k}\left[\frac{\varepsilon_{\bf k}+gn_{0}}{2E_{\bf k}}\coth\Big{(}\frac{\beta E_{\bf k}}{2}\Big{)}-\frac{1}{2}\coth\Big{(}\frac{\beta\varepsilon_{\bf k}}{2}\Big{)}\right].$$

Full text

I Introduction

Bogoliubov’s description of the quantum-mechanical ground state of dilute gases of interacting bosons Bogoliubov1947 forms the basis of our understanding of these systems. Lee1957a ; LHY1957 ; Lee1957b ; Bruckner1957 ; Beliaev1958a ; Beliaev1958b ; Hugenholtz1959 ; Gavoret1964 ; Popov1965 ; Sawada1959 ; Hohenberg1965 ; Singh1967 ; Cheung1971 ; Szepfalusy1974 ; Wong1974 ; Bijlsma1997 ; Shi1998 ; Zagrebnov2001 ; Capogrosso2010 Yet, despite its widespread acceptance and aura of venerability, it has long been recognized that several aspects of Bogoliubov’s theory were not very well understood.Girardeau1959 ; Misawa1960 ; Gardiner1997 ; Girardeau1998 ; Castin1998 ; Leggett2001 ; Andersen2004 ; Yukalov2006 ; Suto2008 ; LeggettNJP ; Kita2010 ; Bobrov2016

In a recent paper, Ettouhami2012 we analyzed Bogoliubov’s theory in quite some detail, and discussed several problematic aspects of this theory, which originate in the so-called Bogoliubov prescription on one hand, and in the non-conservation of the number of bosons on the other. In particular, we have shown that Bogoliubov’s theory does not model a unique system of $N$ interacting bosons, but rather describes an independent collection of such systems, where in each one of these systems $N$ bosons are allowed to be in either one of the three single particle states with momenta ${\bf 0,\pm k}$. Indeed, in Bogoliubov’s theory, a truncated version $\hat{H}_{B}$ of the total Hamiltonian $\hat{H}$ of the system is considered, and this truncated Hamiltonian $\hat{H}_{B}$ is written as a sum $\hat{H}_{B}=\sum_{\bf k\neq 0}\hat{H}_{\bf k}$, where each $\hat{H}_{\bf k}$ describes a sytem of $N$ bosons that can only be in one of three momentum states: $-{\bf k},{\bf 0}$, and ${\bf k}$. Then, the Hamiltonians $\hat{H}_{\bf k}$ are diagonalized independently from one another, leading to a theory that describes, as we mentioned above, not a unique system of $N$ interacting bosons, where the same particle in the condensed state with ${\bf k}=0$ can be excited to any one of the other single-particle states with ${\bf k}\neq 0$, but a juxtaposition of independent systems, each with its own reservoir of $N$ bosons, and where a boson in the condensate with ${\bf k}={\bf 0}$ can only be excited to one of the single particle states $\pm{\bf k}$. Reformulating the theory so as to describe a unique system of $N$ bosons such that a boson in the condensate with ${\bf k}={\bf 0}$ can now be excited to any one of single-particle states with momenta $\{\pm{\bf k}_{1},\ldots,\pm{\bf k}_{\infty}\}$, and taking care to preserve the conservation of particle number, we have explicitly shownEttouhami2012 that the depletion of bosons is significantly reduced with respect to the Bogoliubov case, and that the elementary excitations of the truncated Hamiltonian $\hat{H}_{B}$ become gapped, in contradiction with the results of the standard Bogoliubov theory.

As has been noted elsewhere, Suto2008 a gapped excitation spectrum does not necessarily indicate that the theory is incorrect or that there is a problem with the underlying variational procedure. Actually, most number-conserving theories tend to predict a gapped excitation spectrum. Suto2008 Of course, experimental evidence dictates that in a complete theory the elementary excitations of a uniform, translationally invariant system of interacting bosons should not have a gap. This, however, should not induce us to turn a blind eye to all the conceptual difficulties of Bogoliubov’s theory and accept it as “the” correct formulation merely because it predicts a gapless excitation spectrum. To do so would violate the integrity of the whole process of scientific inquiry, in which the end does not justify the means, and every approximate theory needs to be conceptually sound before its results can be accepted. So, rather than continuing to use an imperfect formulation which predicts a gapless spectrum using highly questionable intermediate steps, a proper approach would consist in using a theory with sound foundations and no internal inconsistencies even if such a theory predicts elementary excitations that have a finite gap, and work toward improving this theory to eventually reach a refined version that predicts no gap. After all, progress in research is most of the time incremental, and there is a better chance of a conceptually sound gapless theory to emerge from a correctly formulated gapped approach than from an incorrectly formulated gapless one.

In this article, we want to resume the project we started in Ref. Ettouhami2012, and re-examine the way in which Bogoliubov’s theory of a dilute Bose gas at $T=0$ has been extended to describe the statistical mechanics of interacting bosons at finite temperatures. Here again, we find upon close scrutiny that several aspects of the conventional formulation of the statistical mechanics of Bose systems are quite questionable and do not rest on sound foundations. The main manifestation of this is our finding that what usually passes as the grand-canonical treatment of the Bogoliubov Hamiltonian of interacting bosons is not a grand-canonical treatment at all, and is in fact a canonical one. A clue that betrays the canonical nature of the standard Bogoliubov-Beliaev-Popov (BBP) theory at finite temperature is the presence of the density of condensed bosons $n_{0}(T)$ in the expressions of the thermodynamic quantities derived within this approach when in fact in a grand-canonical formulation the number of condensed bosons $N_{0}$ should be traced over and should not even appear in these expressions. Another clue is the fact that the density of bosons $n_{B}$ in these formulations is taken to be a constant that does not vary with temperature, while in a grand-canonical description the average density of bosons is not fixed, but depends on other thermodynamic variables, most notably the temperature. Below, we will show that the main results of the BBP approach can be derived within a purely canonical formulation, and we show how one can extend this theory to the grand-canonical case where the system is in contact with a particle reservoir and the average density of bosons $n_{B}$ is not a constant, but depends on the chemical potential $\mu$, the volume $V$ and the temperature $T$, $n_{B}=n_{B}(\mu,V,T)$, and explicitly derive an expression of this dependence.

The rest of this paper is organized as follows. In Sec. II we review the standard formulation of the statistical mechanics of a dilute Bose gas. This formulation is based on the use of the so-called Bogoliubov prescription, where the imaginary-time boson field $\Psi({\bf r},\tau)$ is approximated as the sum of a condensed part $\sqrt{n_{0}}$, treated as a constant in $({\bf r},\tau)$-space, and a fluctuating field $\psi({\bf r},\tau)$ representing depleted bosons:

[TABLE]

In Ref. Ettouhami2012, , we saw that the zero temperature version of this prescription was the root cause of many shortcomings and inconsistencies of Bogoliubov’s theory. Here, we will show that the finite temperature version of Bogoliubov’s prescription creates one more major inconsistency by artificially leading to the appearance of a jump discontinuity in the density of condensed bosons $n_{0}(T)$ at the critical transition temperature $T_{c}$. A discussion of how this happens and how this problem can be avoided will be presented in Sec. III, where we argue that the quantity $n_{0}$ inside the square root on the rhs of Eq. (1) should not depend on temperature, because it is the result of the ad-hoc prescription (1) and not the result of a thermal averaging process. We furthermore show that interpreting $n_{0}$ on the rhs of Eq. (1) as the density of condensed bosons at zero temperature leads to a critical transition temperature $T_{c}$ which is higher than the transition temperature $T_{c0}$ of an ideal Bose gas, with the difference $\Delta T_{c}=T_{c}-T_{c0}$ varying with the parameter $n_{B}a^{3}$, where $n_{B}$ is the density of bosons and $a$ the $s$-wave scattering length, according to a relation of the form:

[TABLE]

with an exponent:

[TABLE]

We then show in Sec. IV that the formulation of the BBP theory of Sec. II is in fact canonical in nature, and not grand-canonical as is commonly thought. We do so by going back to first principles, and deriving fundamental results of this theory by taking the trace over the orthonormal basis of eigenfunctions of the system at a fixed boson number $N$. When the theory is reformulated in this fashion, where we avoid the use of the sophisticated but not too transparent apparatus of field-theory, we are led to realize that the results of the standard theory can be recovered by taking a canonical trace, and hence previous formulations of the statistical mechanics of interacting bosons are all essentially canonical. In Sec. V, we generalize the variational theory of Ref. Ettouhami2012, to include Fock interactions between depleted bosons, and we find the zero temperature ground state and elementary excitations of the system when these Fock interactions between depleted bosons are taken into account. Armed with the knowledge of the energy levels of the system, in Sec. VI we take canonical traces over these levels to derive the thermodynamics of the interacting Bose gas in the canonical ensemble. In Sec. VII we show how the statistical mechanics of interacting bosons can be formulated properly within the grand-canonical ensemble by doing the extra step of taking the trace over the toal number of bosons $N$, and in Sec. VIII we present our conclusions.

II Bogoliubov prescription for an interacting Bose gas

Let us start our study by reviewing the standard formulation of the statistical mechanics of a condensed Bose gas. The traditional formulation uses imaginary-time equations of motion for Heisenberg field operators to derive a set of Dyson equations for the normal and anomalous Green’s functions. Beliaev1958a ; FetterWalecka ; Shi1998 Modern formulations use the path-integral formalism to derive the same results, and that is the formalism we shall use in the rest of this Section.

Our starting point will be the imaginary-time grand-canonical partition function of the system, which is given by:

[TABLE]

where $S$ is the imaginary-time action:

[TABLE]

In the above equation, $\hbar$ is Planck’s constant $h$ divided by $2\pi$, $m$ is the mass of the bosons and $\mu$ is the chemical potential. On the other hand, $\beta=1/(k_{B}T)$, where $k_{B}$ is Boltzmann’s constant and $T$ is the temperature. In this Section, we shall assume for simplicity that the interaction potential $V(\bf r)$ between bosons is repulsive, and that it can be approximated by a short-range, point-contact interaction of the form $V({\bf r})=g\delta({\bf r})$, where the interaction strength $g$ can expressed in terms of the $s$-wave scattering length $a$ through the relation:

[TABLE]

Below the critical temperature $T_{c}$ the zero momentum state is macroscopically populated, and the boson field $\Psi({\bf r},\tau)$ is rewritten as the sum:

[TABLE]

where $n_{0}$ is the density of the condensate, and where the field $\psi({\bf r},\tau)$ describes excited bosons. If we denote by $\{\varphi_{\bf k}({\bf r})\}$ the set of eigenfunctions of the non-interacting Hamiltonian $\hat{H}_{0}=-\hbar^{2}\nabla^{2}/2m$ such that $\varphi_{\bf k}({\bf r})=e^{i{\bf k}\cdot{\bf r}}/\sqrt{V}$ and $\hat{H}_{0}\varphi_{\bf k}({\bf r})=(\hbar^{2}k^{2}/2m)\varphi_{\bf k}({\bf r})$, then $\psi({\bf r},\tau)$ can be written as a Fourier expansion of the form

[TABLE]

The function $\psi({\bf k},\tau)$ can in turn be expanded in a Matsubara Fourier series,

[TABLE]

with $\omega_{n}=2\pi n/\beta\hbar$, so that $\psi({\bf r},\tau)$ becomes:

[TABLE]

Using the prescription (7) and the decomposition (10) in the expression of the action, and noticing that $\int d{\bf r}\psi({\bf r},\tau)=0$, we find that $S$ can be written as a sum of a part $S_{0}$ that is quadratic in $\psi$ and a part $S_{1}$ that contains cubic and quartic contributions. These are given by:

[TABLE]

At this point, it is convenient to write the quadratic part $S_{0}$ using matrix notation in the form:

[TABLE]

where the inverse Green’s function matrix $G^{-1}({\bf k},\omega_{n})$ is given by:

[TABLE]

with:

[TABLE]

Inverting $G^{-1}$ gives the following expression of the Green’s function matrix $G({\bf k},\omega_{n})$:

[TABLE]

where we defined the enrgy spectrum $E_{\bf k}$ such that:

[TABLE]

Requiring this excitation spectrum to vanish as $k\to 0$ imposes the constraint:

[TABLE]

upon which the expression of $E_{\bf k}$ becomes:

[TABLE]

Using the expression of the quadratic action $S_{0}$ from Eq. (14) and the Green’s function matrix $G({\bf k},\omega_{n})$ from Eq. (21), one can derive the following results for the so called normal and anomalous thermal averages, $\langle\psi({\bf k}_{1},\omega_{n})\psi^{*}({\bf k}_{2},\omega_{m})\rangle$ and $\langle\psi({\bf k}_{1},\omega_{n})\psi({\bf k}_{2},\omega_{m})\rangle$, respectively:

[TABLE]

It then follows that the average number of bosons $N_{\bf k}=\langle\hat{N}_{\bf k}\rangle$ with momentum ${\bf k}$ is given by:

[TABLE]

The density of depleted bosons $n_{d}$ is given by:

[TABLE]

and transforming the sum into an integral, we obtain:

[TABLE]

From the above expession of the density of depleted bosons, one can obtain the density of condensed bosons

[TABLE]

With the knowledge of $n_{0}(T)$, the normal and anomalous Green’s functions $G_{11}$ and $G_{12}$ are now fully specified and can be used to derive the thermodynamic properties of the Bose gas.

Now, in previous literature the thermodynamics of condensed Bose systems has been studied using the method summarized above mainly at temperatures well below the transition temperature $T_{c}$. This is because this formulation predicts a jump discontinuity in the condensate density $n_{0}(T)$ at $T_{c}$. In the following Section, we want to examine the origin of this discontinuity, which we will find can be traced back to the Bogoliubov prescription used in this formulation, and more precisely to the improper way this prescription is used at finite temperatures. This will be done next.

III Origin of the jump discontinuity in the density of condensed bosons at $T=T_{c}$

As mentioned above, one of the salient features of the standard formulation of the statistical mechanics of interacting bosons is the fact that the density of condensed bosons $n_{0}(T)$ has a discontinuityShi1998 at the critical transition temperature $T_{c}$. We now show that this discontinuity is a direct consequence of the inappropriate use of the Bogoliubov prescription, Eq. (7), at finite temperatures. Indeed, when we write the decomposition $\Psi({\bf r},\tau)\simeq\sqrt{n_{0}}+\psi({\bf r},\tau)$ in the expression of the action $S[\Psi,\Psi^{*}]$, $\sqrt{n_{0}}$ is merely the $\Psi({\bf k}=0,\omega_{n}=0)$ component of $\Psi({\bf r},\tau)$. As such, it is not the result of a thermal averaging process, and hence it does not and cannot represent the number of condensed bosons at any finite temperature $T$. In fact, the quantity $n_{0}$ that appears in the decomposition (7) can, at the most, be interpreted as the density of condensed bosons at $T=0$. Interpreting $n_{0}$ as $n_{0}(T)$ in Eq. (7) is not only unjustified, it also has an unphysical and rather annoying consequence in that it causes the density of codensed bosons to have a jump discontinuity at the critical transition temperature $T_{c}$. In the following, we shall discuss how this phenomenon takes places, and how correctly interpreting $n_{0}$ in Eq. (7) as $n_{0}(T=0)$ eliminates this jump discontinuity and leads to a continuous variation of $n_{0}(T)$ across the Bose condensation critical temperature $T_{c}$.

Let us for definiteness rewrite here the density of depleted bosons at temperature $T$ obtained in the previous Section, Eq. (28) – (note that we use the shorthand notation $\int_{\bf k}\equiv\int d{\bf k}/(2\pi)^{3}$):

[TABLE]

Close to the transition temperature $T_{c}$, we shall follow Ref. Shi1998, and calculate the difference:

[TABLE]

Here $n_{cr}(T)$ is the density of excited particles of an ideal Bose gas at a given temperature $T$, and is given by:

[TABLE]

where $\zeta(x)$ is Riemann’s Zeta function, and $\zeta(3/2)\simeq 2.612$. To calculate the difference in Eq. (30), we note that the main contribution to the integral comes from the region of ${\bf k}$-space where the energy spectrum $E_{\bf k}$ differs significantly from the free-particle energy $\varepsilon_{\bf k}$, which in this case is the region $\varepsilon_{\bf k}\leq gn_{0}$. Assuming that $gn_{0}\ll k_{B}T_{c}$, we can approximate $\coth(\beta E_{\bf k}/2)$ by $(2/\beta E_{\bf k})$ and $\coth(\beta\varepsilon_{\bf k}/2)$ by $(2/\beta\varepsilon_{\bf k})$ to obtain:Shi1998

[TABLE]

We now need to use the relationship between the total number of bosons $n_{B}$, the number of condensed bosons $n_{0}(T)$ and the number of depleted bosons $n_{d}(T)$,

[TABLE]

to eliminate $n_{d}(T)$ from Eq. (32) and rewrite it solely in terms of $n_{B}$ and $n_{0}(T)$, with the result:

[TABLE]

In the standard formulation of the statistical mechanics of an interacting Bose gas, the quantity $n_{0}$ inside the square root, which is the one coming directly from the Bogoliubov prescription of Eq. (7), is taken to be a function of temperature. In this case, Eq. (34) can be seen as a quadratic equation for $\sqrt{n_{0}(T)}$ which has two solutions:

[TABLE]

Here, we followed Ref. Shi1998, and called $\sqrt{n_{g}(T)}$ the coefficient of $\sqrt{n_{0}}$ in Eq. (34), namely:

[TABLE]

Below the transition temperature, $n_{cr}<n_{B}$, and so the solution $n_{0}^{(-)}(T)$ in Eq. (34) is negative and hence unphysical. We then can immediately see that the physical solution

[TABLE]

being the square of a strictly positive number, does not vanish for any value of the temperature $T$. At $T=T_{c}$, $n_{cr}(T_{c})=n_{B}$ and hence $n_{0}(T)$ is discontinuous at the transition:Shi1998

[TABLE]

At this point, we will note that a similar jump discontinuity of $n_{0}(T)$ at $T=T_{c}$ is also found to take place in the formulation of Hartree-Fock theory Zhang2006 that uses Bogoliubov’s prescription with a temperature-dependent $n_{0}$ (more on this below).

Having reviewed how the standard derivation gives rise to a jump discontinutiy in the density of condensed bosons at $T=T_{c}$, we now go back to Eq. (34) and observe that the $n_{0}$ inside the square root of this last equation is the same $n_{0}$ that appears in the expression of the Green’s function $G_{11}({\bf k},\omega_{n})$, and, as such, should be independent of temperature (the other $n_{0}(T)$ on the rhs of Eq. (34) comes from using the relation between $n_{B}$, $n_{0}$ and $n_{1}$, Eq. (33), and has to have a temperature dependence). In other words, we will rewrite the Bogoliubov prescription of Eq. (7) in the form:

[TABLE]

where $\Psi_{0}$ does not depend on temperature because it is the result of the ad-hoc prescription (39) and not the result of a thermal averaging process. Then the action of Eq. (14) and hence the Green’s function matrix of Eq. (21) will all be written in terms of the temperature-independent quantity $\Psi_{0}$ (instead of $\sqrt{n_{0}}$ which one can easily be misled to interpret as a quantity that varies with temperature), and so would the rhs of Eq. (28) giving the density of depleted bosons at temperature $T$. In other words, Eq. (28) should now read:

[TABLE]

with the energy spectrum

[TABLE]

where $|\Psi_{0}|^{2}$ on the rhs of these last two equations does not depend on temperature. Comparing the excitation spectrum (41) with the standard $T=0$ expression $E_{\bf k}=\sqrt{\varepsilon_{\bf k}(\varepsilon_{\bf k}+2gn_{0}(T=0))}$, we see that we can identify $|\Psi_{0}|^{2}$ to be the density of condensed bosons at zero temperature:Remark

[TABLE]

It then follows that a more appropriate way to write Eq. (34) is as follows:

[TABLE]

where we now write the density of condensed bosons inside the square root as $n_{0}(0)$ to emphasize that this is the quantity that is coming directly from the decomposition (39) in the imaginary-time action, and hence is the $T=0$ value of the density of condensed bosons (below in Sec. V we will show that this quantity should in fact be $n_{B}$). Solving for $n_{0}(T)$, we obtain:

[TABLE]

where in going from the first line to the second one we used the expression (31) of $n_{cr}(T)$, and we approximated $n_{0}(0)$ inside the square root with the total density of bosons $n_{B}$. Dividing both sides of the above equation by $n_{B}$, we can write:

[TABLE]

At $T=T_{c}$, we expect $n_{0}(T_{c})=0$. Assuming this to be true, we can write:

[TABLE]

The critical temperature for a free Bose gas $T_{c0}$ can be obtained by letting $g=0$ in the above equation, and satisfies the following relation:

[TABLE]

which gives the well-known result:

[TABLE]

Going back to the condensate fraction of Eq. (45), if we scale the temperature $T$ by $T_{c0}$ we obtain after a few manipulations:

[TABLE]

At $T=T_{c}$, $n_{0}(T)=0$, and so:

[TABLE]

Now, let $\Delta T_{c}=T_{c}-T_{c0}$, so that $T_{c}=T_{c0}(1+\Delta T_{c}/T_{c0})$. Then, if we assume that the difference between critical temperatures is small, i.e. $\Delta T_{c}\ll T_{c0}$, then we can write:

[TABLE]

and Eq. (50) has the approximate solution:

[TABLE]

where we call $c$ the numerical constant:

[TABLE]

We now want to check this result numerically. Going back to Eq. (28), we see that we can write for the condensate fraction $n_{0}(T)/n_{B}=1-n_{d}(T)/n_{B}$ the following expression:

[TABLE]

where we introduced the dimensionless wave-vector $\tilde{\bf k}$ such that:

[TABLE]

with

[TABLE]

In Fig. 1, we plot the condensate fraction as obtained from the above equation for two different values of the parameter $n_{B}a^{3}$. We find that for $n_{B}a^{3}=0.001$, the condensate fraction vanishes at $T_{c}/T_{c0}=1.33$, and for $n_{B}a^{3}=0.01$, the condensate fraction vanishes at $T_{c}/T_{c0}=1.39$.

The expression above for the condensate fraction, Eq. (54), can be used to study the dependence of $T_{c}$ on $n_{B}a^{3}$ numerically by looking for the temperature $T$ that solves the equation $n_{0}(T)/n_{B}=0$. The upper panel of Fig. 2 shows a plot of $\Delta T_{c}/T_{c0}$ vs. $n_{B}a^{3}$ found using this method, where the infinite slope near the origin indicates a power-law relationship with an exponent that is less than unity. In the lower panel, we plot $\ln(\Delta T_{c}/T_{c0})$ vs. $\ln(n_{B}a^{3})$, which shows that there is a linear relationship between these two logarithms with a slope of $0.158\simeq 1/6$, in perfect agreement with the analytical result of Eq. (52).

It is worth taking a pause at this juncture to note the vast body of literature which treated the problem of the shift in the transition temperature due to repulsive interactions between bosons. Beginning with the early work of Lee and Yang, Lee1957b various methods of calculations over the years have produced widely dissimilar results, Glassgold1960 ; Luban1962 ; FetterWalecka ; Kanno1969a ; Kanno1969b ; Toyoda1982 ; Stoof1992 ; Bijlsma1996 ; Huang1999 ; Mueller2001 ; Gruter1997 ; Holzmann1999a ; Holzmann1999b ; Baym1999 ; Baym2000 ; Wilkens2000 ; Arnold2000 ; deSouza2001 ; Kashurnikov2001 ; Arnold2001 ; Holzmann2001 ; Arnold2002 ; Baym2001 ; deSouza2002 ; Kneur2002 ; Braaten2002 ; Kastening2003 with positive and negative shifts in $T_{c}$ that were predicted to be proportional to $a$, $\sqrt{a}$, $a\ln(a)$, etc. It is not our intention here to give support to one of these results or another, our main focus in this paper not being to give an accurate or definitive calculation of the shift in transition temperature, but merely to correct a misconception about the standard Beliaev BBP predicting a jump discontinuity of $n_{0}(T)$ at the transition. It is interesting to note that the positive shift we obtain

[TABLE]

with the exponent

[TABLE]

has been predicted previously, Lee1957b ; Glassgold1960 ; Huang1999 although the author is not aware of this result having been obtained within an approach similar to the one used in the current study.

Now that we got rid of the jump discontinuity in the condensate fraction at the critical temperature $T_{c}$ predicted by the standard, field-theoretic formulation of the statistical mechanics of interacting bosons, one may ask how we should interpret the quantity $n_{B}$ that appears in Eq. (54). Indeed, in the grand-canonical ensemble, where the system is in contact with a particle reservoir, the density of bosons $n_{B}$ should depend on temperature. Hence, the question arises as to whether we should interpret $n_{B}$ in Eq. (54) to be the total density of bosons (i) at zero temperature, (ii) at the critical temperature $T_{c}$, or (iii) at some intermediate temperature between $T=0$ and $T=T_{c}$. The fact is, the plots in Figs. 1 and 2 were drawn with the assumption that $n_{B}$ does not depend on temperature, and even if we are willing to assume a temperature dependence for $n_{B}$, there is simply not enough information in the BBP formulation to find what this temperature dependence exactly is.

In the following, we want to examine the physical content of the BBP formulation by trying to derive some of its major results using first-principles, i.e. using traces over eigenstates of the Bogoliubov Hamiltonian. In doing so, we will establish that what passes as a grand-canonical treatment of the interacting Bose gas is in fact a canonical treatment, not a grand-canonical one as is commonly thought, and hence that we do not have to worry about the temperature dependence of $n_{B}$ in Eq. (54) because this quantity, in the canonical ensemble, is fixed and does not in any way vary with temperature.

IV Statistical mechanics of interacting bosons: canonical formulation

IV.1 Back to basics: thermal averages as traces over eigenvectors of the Bogoliubov Hamiltonian

In this Section, we want to derive fundamental building blocks of the BBP theory, namely the expressions of the normal and anomalous thermal averages $\langle a_{\bf k}^{\dagger}a_{\bf k}\rangle$ and $\langle a_{\bf k}a_{-\bf k}\rangle$ respectively, using a first-principles type of approach. Our starting point will be the canonical partition function of the Bogoliubov model, which for a system of $N$ bosons is given by (note that we are purposely adding a subscript $N$ to the symbol Tr to indicate that we are taking a canonical trace with a fixed number of bosons $N$):

[TABLE]

In the above equation, $\hat{H}_{B}$ is the truncated part of the total Hamiltonian $\hat{H}$ that can be diagonalized using Bogoliubov’s canonical transformations, and is given in Fourier space by the following expression:

[TABLE]

After the above Hamiltonian is diagonalized, it can be written in the form:

[TABLE]

where $E_{0}$ is the ground state energy of the system, and $E_{\bf k}$ is the energy of an excitation of wavevector ${\bf k}$. In the standard Bogoliubov formulation, these two quantities are given by:

[TABLE]

On the other hand, $\alpha_{\bf k}^{\dagger}$ and $\alpha_{\bf k}$ are the creation and annihilation operators of an elementary excitation of momentum ${\bf k}$:

[TABLE]

where the so-called coherence factors $u_{\bf k}$ and $v_{\bf k}$ are given by:

[TABLE]

We shall denote by $|\Psi_{B}(N)\rangle$ the ground state of Bogoliubov’s Hamiltonian $\hat{H}_{B}$, such that $\alpha_{\bf k}|\Psi_{B}(N)\rangle=0$ for all wavevectors ${\bf k}\neq 0$. It then follows that:

[TABLE]

confirming our statement above that $E_{0}$ in Eq. (61) is the ground state energy of $\hat{H}_{B}$. On the other hand, we shall denote by $|\Psi_{\{m_{i}\}}(N)\rangle=|m_{1},m_{2},\ldots,m_{M}\rangle$ the eigenfunction of $\hat{H}_{B}$ corresponding to $m_{i}$ excitations of wavevector ${\bf k}_{i}$, $i=1,\ldots,M$ ($M$ here is the number of momentum modes kept in the Hilbert space, and will eventually be sent to infinity). These eigenfunctions can be written in the form:

[TABLE]

and one can easily verify by direct calculation that $|\Psi_{\{m_{i}\}}(N)\rangle$ satisfies the following equation:

[TABLE]

We now go back to Eq. (59) and rewrite the canonical partition function $Z_{C}(N,V,T)$ in the form:

[TABLE]

In the above expression, the summation extends over all values of $m_{i}$ that are smaller than the dimension of the matrix representation of the Hamiltonian $\hat{H}_{B}$. In the thermodynamic limit $N\to\infty$ and $M\gg 1$, the summations over the $m_{i}$’s can be extended all the way to infinity, with the result:

[TABLE]

where, in going from the second to the third line use has been made of the geometric summation formula:

[TABLE]

Now, the thermal average of the operator $\hat{N}_{{\bf k}_{i}}=a_{{\bf k}_{i}}^{\dagger}a_{{\bf k}_{i}}$ in the canonical ensemble is given by:

[TABLE]

Although it is possible to conduct the calculation of the above thermal average in the number-conserving formalism of Ref. Ettouhami2012, , we here for simplicity shall use the number non-conserving formalism in which the relation between the operators $a_{{\bf k}_{i}}$ and $\{\alpha_{{\bf k}_{i}},\alpha_{{\bf k}_{i}}^{\dagger}\}$ is linear, and can be obtained by inverting Eqs. (63), with the result:

[TABLE]

The operator $\hat{N}_{\bf k}=a_{\bf k}^{\dagger}a_{\bf k}$ is therefore given by:

[TABLE]

Let us find the trace of ${\hat{N}}_{\bf k}\exp(-\beta\hat{H}_{B})$ in the base formed by the $|\Psi_{\{m_{i}\}}\rangle$. In order to do that, we shall first use Eq. (67) to calculate the action of this operator on the ket $|\Psi_{\{m_{i}\}}\rangle$, with the result:

[TABLE]

Hence:

[TABLE]

Using the expression of $\hat{N}_{{\bf k}_{j}}$ in terms of the $\alpha_{\bf k}$’s, Eq. (73), in the above expression, we obtain:

[TABLE]

We now use the fact that $\alpha_{{\bf k}_{j}}^{\dagger}\alpha_{{\bf k}_{j}}$ simply counts the number of excitations of wavevector ${\bf k}_{j}$ in the excited state $|\Psi_{\{m_{i}\}}(N)\rangle$, which is nothing but $m_{j}$, to write:

[TABLE]

Also, from the canonical commutation relations of the $\alpha_{\bf k}$’s, we see that $\alpha_{-{\bf k}_{j}}\alpha_{-{\bf k}_{j}}^{\dagger}=1+\alpha_{-{\bf k}_{j}}^{\dagger}\alpha_{-{\bf k}_{j}}$, and hence:

[TABLE]

where we denote by $m_{-j}$ the number of excitations of wavevector $-{\bf k}_{j}$. Now, if we use Eqs. (77) and (78) in Eq. (76), this last equation becomes:

[TABLE]

Replacing this result back into Eq. (75) and taking the trace, we obtain:

[TABLE]

Dividing by ${Z_{C}}$, as in Eq. (71), we obtain that the thermal average of $\hat{N}_{\bf k}$ in the canonical ensemble is given by:

[TABLE]

Recalling that $E_{-{\bf k}_{j}}=E_{{\bf k}_{j}}$, this last equation becomes:

[TABLE]

Using Eq. (70) and the fact that

[TABLE]

we finally obtain:

[TABLE]

Now, using the expression of the coherence factors $u_{{\bf k}_{i}}^{2}$ and $v_{{\bf k}_{i}}^{2}$ from Eq. (64), we can rewrite $\langle\hat{N}_{{\bf k}_{j}}\rangle$ in the form:

[TABLE]

This is the exact same expression that was obtained previously from the field-theortic formulation of the statistical mechanics of an interacting Bose gas, see Eq. (26b). The fact that we obtained this expression within a strictly canonical formalism, where we traced over the basis of eigenstates $|\Psi_{\{m_{i}\}}(N)\rangle$ corresponding to a fixed number of bosons $N$, is quite remarkable. In fact, and as we already mentioned in the beginning of this section, a quite conspicuous clue that the usual treatment is canonical is the presence of the density of bosons $n_{B}$ on the rhs of Eq. (85), for in a correctly formulated grand-canonical treatment the number of bosons $N$ would be traced over, and the rhs of this last equation would no longer contain the density of bosons $n_{B}=N/V$. (Strictly speaking, in the standard Bogoliubov derivation $n_{0}$ appears on the rhs of Eq. (85) rather than $n_{B}$, and the same argument can be made that in a grand-canonical formulation $n_{0}$ should be traced over and hence should not appear on the rhs of this last equation.)

We now want to calculate the anomalous thermal average $\langle a_{{\bf k}_{j}}a_{-{\bf k}_{j}}\rangle$. Following the analysis we made in Ref. Ettouhami2012, of the meaning of this anomalous average at zero temperature, where we showed that it can be interpreted as the expectation value $\langle\Psi(N-2)|a_{{\bf k}_{j}}a_{-{\bf k}_{j}}|\Psi(N)\rangle$, $|\Psi(N)\rangle$ and $|\Psi(N-2)\rangle$ being the Bogoliubov ground states of a system of $N$ and $(N-2)$ bosons respectively, we here will define the canonical trace in the following way:

[TABLE]

Using Eq. (67), we can write:

[TABLE]

We now use Eq. (72) to express $a_{\bf k}a_{-\bf k}$ in terms of the $\alpha_{\bf k}$’s, with the result:

[TABLE]

Taking the quantum expectation value, we obtain:

[TABLE]

Replacing this result in Eq. (87), and using expression (69c) of the canonical partition function $Z_{C}$, we can write:

[TABLE]

Performing the summations with the help of Eqs. (70) and (83), we finally obtain:

[TABLE]

This is again the exact same expression as the one obtained from the standard Bogoliubov theory,Shi1998 ; FetterWalecka except that the total density of bosons $n_{B}$ appears on the rhs instead of $n_{0}$. We emphasize that this is so because we used $n_{B}$ instead of $n_{0}$ in the expressions of the coherence factors $u_{\bf k}$ and $v_{\bf k}$, Eq. (64).

The fact that the thermal averages $\langle a_{\bf k}^{\dagger}a_{\bf k}\rangle$ and $\langle a_{\bf k}a_{-\bf k}\rangle$ that are derived in the BBP approach can be derived using canonical traces is quite remarkable. It tells us that thermal averages of physical observables (one or two-body operators that can be expressed as combinations of $a_{\bf k}^{\dagger}a_{\bf k}$ and $a_{\bf k}a_{-\bf k}$ and their hermitian conjugates) can be obtained in the canonical ensemble, and hence that the whole BBP formulation is canonical in nature, and not grand canonical. In the following subsection, we want to analyze the various tracing schemes that one can use to obtain partition functions for interacting bosons, and try to understand how in the BBP, starting from a supposedly grand-canonical partition function we end up with results that can be derived in the canonical ensemble.

IV.2 Analyzing the various possible tracing schemes

We now are in a good position to discuss what exactly is not quite right about the standard field-theoretic formulations of the statistical mechanics of interacting bosons. These formulations are so elegant and so mathematically sophisticated that it is hard to imagine that they actually do anything other than what they purport to be doing. Unfortunately, mathematical sophistication is no guarantee of correctness, and below we will show that in the standard formulation of BBP’s theory, the number of condensed bosons $N_{0}$ is never traced over, and this incomplete tracing operation effectively reduces what is supposed to be a grand-canonical trace to a canonical one instead.

Let us for example have a detailed look at the manipulations that led us to Eq. (21). When we write the action $S_{0}$ as a quadratic form in the fields $\psi$ and $\psi^{*}$, as we did in Eq. (14), and take the inverse of the matrix $G^{-1}$ to find the Green’s function matrix $G$, these two quick steps are thought to be equivalent to the following ones in first-quantized methods:

The eigenvalues and eigenfunctions of the Bogoliubov Hamiltonian are found for a system of $N$ bosons; 2. 2.

A canonical trace of the Boltzmann weight $e^{-\beta(\hat{H}-\mu\hat{N})}$ is taken over these eigenfunctions at constant number of bosons $N$; 3. 3.

A second trace is taken over the total number of bosons $N$ from $N=0$ to $N=\infty$ to find the grand-canonical partition function.

The fact that the above three steps do actually take place behind the scenes is never doubted and is generally taken for granted. It turns out, unfortunately, that taking the inverse of the matrix $G^{-1}$ to go from Eq. (18) to Eq. (21) succeeds in completing steps 1 and 2 of the above program, but fails to complete step 3. In other words, finding the Green’s function matrix $G$ is equivalent to taking the trace over the eigenfunctions $|\Psi_{\{m_{i}\}}(N)\rangle$ of the Bogoliubov Hamiltonian $\hat{H}_{B}$ at fixed boson number $N$, which is a canonical trace, but is not equivalent to taking the extra step of doing a trace over $N$ that would make the approach grand-canonical.

We now want to discuss the various kinds of trace that one may use to evaluate partition functions in the canonical and grand-canonical ensembles.

IV.2.1 Canonical ensemble: tracing over the basis of single-particle states

Let us first consider the relatively simple situation of a gas of $N$ bosons in the canonical ensemble. The canonical partition function ${Z_{C}}=\mbox{Tr}_{N}\left(e^{-\beta\hat{H}}\right)$ can be written as the following trace over boson single-particle states:

[TABLE]

where the ket $|N_{0},\ldots,N_{M}\rangle$ represents the state having $N_{0}$ bosons with momentum ${\bf k}_{0}={\bf 0}$, $N_{1}$ bosons with momentum ${\bf k}_{1}$, and so on, i.e.:

[TABLE]

Note how, in Eq. (92), we chose to take the trace of the Hamiltonian $\hat{H}_{B}$ over the complete orthonormal basis composed of the eigenstates of the non-interacting Hamiltonian $\hat{H}_{0}=\sum_{\bf k}\varepsilon_{\bf k}a_{\bf k}^{\dagger}a_{\bf k}$ at fixed number of bosons $N$, this last constraint being enforced by the Kronecker delta on the rhs of this last equation, which ensures that the total number of particles in all the basis kets $|N_{0},N_{1},\ldots,N_{M}\rangle$ is equal to $N$.

IV.2.2 Canonical ensemble: tracing over the basis of elementary excitations

Now, after we diagonalize the Hamiltonian $\hat{H}_{B}$ using a number-conserving formalism such as the one developed in Ref. Ettouhami2012, , we end up with another basis of eigenfunctions describing elementary excitations of the system of $N$ bosons, with one such eigenfunction $|\Psi_{\{m_{i}\}}(N)\rangle=|m_{1},m_{2},\ldots,m_{M}\rangle$ describing a state having $m_{i}$ excitations in the state of momentum ${\bf k}_{i}$, with $i=1,\ldots,M$, see Eq. (66). The trace being independent of the orthonormal basis chosen, we observe that there is a second, equivalent way to calculate the partition function $Z_{C}$ which consists in taking a trace over the basis of eigenfunctions $\{|\Psi_{\{m_{i}\}}(N)\rangle\}$ at fixed total number of bosons $N$:

[TABLE]

In the above equation, the $m_{i}$’s take integer values $m_{i}=0,1,2,\ldots,D(N,M)$, corresponding to the number of excitations of wavevector ${\bf k}_{i}$. The upper limit of the summation over any single quantum number $m_{i}$ is the dimension $D(N,M)$ of the matrix representation of the Hamiltonian $\hat{H}_{B}$ and depends on the number of bosons and the number of momentum modes considered. Assuming that our system has $N$ bosons and $M$ momentum modes, any given state of the Hilbert space can be written in the form $|N_{0},N_{1},N_{2},\ldots,N_{M}\rangle$, and hence the dimension of the Hilbert space is the number of ways we can solve the equation:

[TABLE]

This is a well-known combinatorial problem, which can be mapped onto the problem of finding the number $D(N,M)$ of arrangements of $N$ balls and $M-1$ dividers (all indistinguishable), and is given by:

[TABLE]

It can be shown that $D(N,M)$ becomes extremely large for values of $N$ and $M$ such that $N+M\gg N,M$, and hence for all practical purposes we can extend the summations in Eq. (94) all the way to infinity:

[TABLE]

Note the fundamental difference between Eq. (92) and Eq. (97). For while in the former there is a need for a Kronecker delta on the rhs to ensure that the constraint (95) is satisfied, no such constraint is imposed on the $m_{i}$’s. This has the important consequence that while it is notoriously difficult to calculate the canonical trace directly from Eq. (92), precisely because of the constraint $\sum_{i}n_{i}=N$ on the rhs of this equation, in the basis of Bogoliubov eigenfunctions there is no such constraint and the sums on the rhs of Eq. (94) can be decoupled into a product of geometric series and hence can be easily evaluated.

IV.2.3 Grand-canonical ensemble: which trace is used in the BBP formalism ?

We now want to discuss the grand-canonical ensemble, where the grand-canonical partition function can be obtained from the canonical one using the following relation:

[TABLE]

Inserting the expression of ${Z_{C}}$ from Eq. (92) into Eq. (98), one can show that the summation over all values of $N$ from $N=0$ to $\infty$ is equivalent to extending the summations over the $N_{i}$’s from [math] to $\infty$ in Eq. (92), with the result:Remark-ideal

[TABLE]

Alternatively, if we start from the expression (94) of the canonical partition function, where the trace is evaluated in the basis of the $|\Psi_{\{m_{i}\}}(N)\rangle$’s, the grand partition function takes the form:

[TABLE]

The two tracing schemes in Eqs. (99) and (100) are mathematicall equivalent and have to produce the same result for the grand partition function $Z_{G}$. Notice from Eq. (99) that all single particle occupation numbers $N_{0},\ldots,N_{M}$ are traced over and should not appear in the expression of $Z_{G}$. And, since the total number of bosons $N$ can be decomposed as the sum $N=N_{0}+\cdots+N_{M}$, $N$ itself should also not appear in the expression of $Z_{G}$. Similarly, since Eq. (100) should give the exact same result at Eq. (99), we conclude that the rhs of Eq. (100) should not depend on any one of the single particle occupation numbers $N_{0},\ldots,N_{M}$, and should also not depend on the total number of bosons $N$.

We now can make an attempt at explaining the reason why the result we obtained in Sec. II for the thermal average \big{\langle}a_{\bf k}^{\dagger}a_{\bf k}\big{\rangle}, which was supposedly derived in the grand-canonical ensemble, coincides with the result of Eq. (85), which was derived in the canonical ensemble (provided $n_{0}$ is replaced with $n_{B}$). This reason has to do with the Bogoliubov prescription, which replaces Ettouhami2012 the density of bosons $n_{B}$ in the expression of the ground state and excitation energies with $n_{0}$ (we shall indeed see in Sec. V below that when the Bogoliubov prescription is avoided, the condensate density $n_{0}$ does not appear in the expressions of the ground state and excitation energies of the system, and that only $n_{B}$ appears in these expressions), and with the fact that the occupation number of the condensate $N_{0}$ is kept as an immutable constant and never traced over. In other words, the BBP formulation is equivalent to the canonical trace in Eq. (97), as the final trace over $N$ in Eq. (100) is never performed in this approach. The fact that the trace is taken in an incomplete way leads to a canonical result instead of a grand-canonical one for the average N_{\bf k}=\big{\langle}a_{\bf k}^{\dagger}a_{\bf k}\big{\rangle} of the number of bosons in the single-particle state of momentum ${\bf k}$.

The realization that the statistical treatments of the interacting Bose gas are canonical in nature and not grand-canonical as is commonly thought is the main result of this paper. In what follows, we want to extend the work we did in Ref. Ettouhami2012, and generalize the variational number-conserving theory we developed in that reference to finite temperatures. However, before we can do that, we need to find a way to take Fock interactions between depleted bosons into account. This will be the subject of the following Section.

V Ground state energy and excitation spectrum of interacting bosons: taking Fock interactions between depleted bosons into account

In Sec. VI below, we will want to examine the canonical formulation of the statistical mechanics of interacting bosons by extending the variational approach developed in Ref. Ettouhami2012, to finite temperatures. The motivation behind our desire to use the variational approach resides in the fact that it will allow us to avoid the use of Bogoliubov’s prescription, and hence will allow us to seamlessly avoid issues such as the jump discontinuity in the condensed density $n_{0}(T)$ at the critical temperature $T_{c}$, and unambiguously demonstrate that $n_{0}$ which appears in the expression of the density of depleted bosons should in fact be replaced by $n_{B}$.

Because we will be interested in the thermodynamics of the Bose gas in the whole range of temperatures between $T=0$ and the critical transition temperature $T_{c}$, we will need to take into account the Fock interactions between depleted bosons, which are neglected in the $T=0$ diagonalization of Bogoliubov’s Hamiltonian. This can be done in two different ways. One way consists in treating these Fock interactions as a perturbation to the quadratic action $S_{0}$, and performing a perturbative expansion of the partition function in terms of this perturbation. Another, more compact way to treat the Fock terms, which avoids the difficulties of perturbation theory and mirrors more closely the inner workings of an exact solution, consists in trying to generalize the variational approach of Ref. Ettouhami2012, to find the ground state and excitation energies in presence of the Fock interactions between depleted bosons, and is the way we are going to follow in this Section. To this end, let us consider the following Hamiltonian:

[TABLE]

Here $\hat{H}_{B}$ is Bogoliubov’s Hamiltonian of Eq. (60), and $\hat{H}_{F}^{\prime}$ is the part of the total Hamiltonian that represents the Fock interaction between depleted bosons,

[TABLE]

Henceforth, we shall write the Hamiltonian $\hat{H}$ in the form:

[TABLE]

where the single-mode Hamiltonian $\hat{H}_{\bf k}$ is given by:

[TABLE]

We shall first find the variational, number-conserving ground state of the total Hamiltonian $\hat{H}$ in Subsection V.1 before deriving the excitation spectrum of this Hamiltonian in Subsection V.2.

V.1 Variational ground state energy

As we mentioned above, our aim in this Subsection is to generalize the variational approach of Ref. Ettouhami2012, to the Hamiltonian $\hat{H}$ of Eqs. (103)-(104), which includes the Fock interaction between depleted bosons. To this end, we shall use for $|\Psi(N)\rangle$ an expression of the form (we again here denote by $M$ the total number of momentum modes kept in the calculation, which will eventually be sent to infinity):

[TABLE]

where the normalized basis wavefunctions are given by:

[TABLE]

As we explained in detail in Ref. Ettouhami2012, , the ground state wavefunction in Eq. (105) is not a simple product of ground state wavefunctions of the single-mode Hamiltonians $\hat{H}_{\bf k}$, as was the case in previous approaches to this problemLee1957a ; HuangBook ; Leggett2001 where the kets $|N-2\sum_{i=1}^{M}n_{i};n_{1},n_{1};\ldots;n_{M},n_{M}\rangle$ on the rhs of Eq. (105) were replaced with $|N_{0};n_{1},n_{1};\ldots;n_{M},n_{M}\rangle$, i.e. the number of bosons in the single-particle state of momentum ${\bf k}=0$ was kept as an immutable constant $N_{0}$ regardless of the number of bosons $n_{1},n_{2},\ldots,n_{M}$ in the other single-particle states ${\bf k}_{i}\neq 0$. Here by contrast, the presence of the quantity \big{[}N-2\sum_{i=1}^{M}n_{i}\big{]} in the various kets that appear on the rhs of Eq. (105) acts like an implicit and rather nontrivial coupling between all the single-mode Hamiltonians $\{\hat{H}_{\bf k}\}$, and goes a step beyond wavefunctions studied in previous literature.

We now will follow the same steps as in Ref. Ettouhami2012, , and use the following ansatz for the coefficients $C_{n_{i}}$:

[TABLE]

It can then be shown that the expectation value of the Hamiltonian $\hat{H}_{{\bf k}_{j}}$ in the state $|\Psi(N)\rangle$ of Eq. (105) is given by the following expression (see Appendix A):

[TABLE]

where $\bar{v}(k_{j})$ and $\tilde{v}(k_{j})$ are given by:

[TABLE]

Minimization of the expectation value in Eq. (108) with respect to the variational constants $\{c_{{\bf k}_{j}}\}$ leads to the following equation:

[TABLE]

In the above equation, $\tilde{\cal E}_{\bf k_{j}}$ denotes the quantity:

[TABLE]

Transforming the sum over momentum modes in Eq. (111b) into an integral,Ettouhami2012 we can rewrite $\sigma_{{\bf k}_{j}}$ in the form:

[TABLE]

Solving Eq. (110) for $c_{\bf k}$, we obtain:

[TABLE]

where the sign of the second term has been chosen so that $0<c_{\bf k}<1$. In the particular case where $v({\bf k})=g$, corresponding to $v({\bf r})=g\delta({\bf r})$ in real space, the ratio $[v({\bf k}-{\bf k}^{\prime})/v({\bf k})]=1$, and the expression of $\sigma_{\bf k}$ can be written in the form:

[TABLE]

where $n_{d}$ is the density of depleted bosons:

[TABLE]

As can be seen, the quantity on the rhs of Eq. (114b) does not depend on ${\bf k}$, and we shall henceforth drop the subscript ${\bf k}$ from $\sigma_{\bf k}$, and rewrite Eqs. (111) in the form:

[TABLE]

In Ref. Ettouhami2012, , we used a Gaussian form for the interaction potential between bosons of the form:

[TABLE]

where $\lambda$ is a positive quantity having the dimensions of length, so as to avoid an ultraviolet divergence in the integral on the rhs of Eq. (116b) in three dimensions, which occurs when $v(k)=g$ is a constant because in that case $c_{\bf k}$ from Eq. (113) behaves like $1/k^{2}$ as $k\to\infty$. In Fourier space, the interaction potential of Eq. (117) is given by:

[TABLE]

Let us now introduce dimensionless units, where energies are measured in units of $n_{B}v({\bf 0})=gn_{B}$, and wavevectors are measured in units of $k_{0}=\sqrt{2mgn_{B}}/\hbar$. Then we can write for $v({\bf k})$ the following expression:

[TABLE]

where in the notation of Eq. (55) $\tilde{\bf k}$ is the dimensionless wavevector $\tilde{\bf k}={\bf k}/{k_{0}}$. As we mentioned above, the reason behind using an interaction potential with Gaussian form in Ref. Ettouhami2012, was to avoid the ultraviolet divergence that results from summing over momentum modes in the integral on the rhs of Eq. (116b). It is not difficult to see from Eq. (119) that such a Gaussian interaction potential $v({\bf k})$ introduces a momentum cut-off around values of $\tilde{k}$ such that $8\pi n_{B}a\lambda^{2}\tilde{k}^{2}\sim 1$. Hence, for reasons of numerical efficiency we’ll find it expedient to replace the above interaction potential with the following approximation:

[TABLE]

where the ultraviolet momentum cut-off $\Lambda$ is given by:

[TABLE]

For simplicity, we shall take the characteristic length scale $\lambda$ which governs the range of the interaction potential $v({\bf r})$ to be the scattering length $a$. The expression of $\Lambda$ becomes:

[TABLE]

We can then write for $c_{\bf k}$ the following expression:

[TABLE]

where we denote by $\tilde{\sigma}$ the dimensionless quantity:

[TABLE]

and where $C_{d}$ is given by:Ettouhami2012

[TABLE]

Here, the fraction of depleted bosons $N_{d}/N$ is given by:Ettouhami2012

[TABLE]

where we used the shorthand notation:

[TABLE]

If we use the expression (123) of $c_{\bf k}$ in Eq. (116b), and change the variable of integration from ${\bf k}$ to the dimensionless wavevector $\tilde{\bf k}={\bf k}/k_{0}$ (see Eq. (56)), we can write for the dimensionless quantity $\tilde{\sigma}$ the following self-consistent equation:

[TABLE]

A numerical solution to the above equation for $\tilde{\sigma}$ can be found by iteration in the following way. First, one starts with an initial guess for $\tilde{\sigma}$. Using this initial guess, one computes the ratio $N_{d}/N$ using Eq. (126) above, which allows us to find the “depletion constant” $C_{d}=1-N_{d}/N$. This value of $C_{d}$ is then used to solve equation (128) for $\tilde{\sigma}$. This computed value of $\tilde{\sigma}$ is then used again as an input to find a better estimate of the ratio $N_{d}/N$ using Eq. (126), and the process is repeated until convergence and a stable solution for $C_{d}$ and $\tilde{\sigma}$ is found.

In the following, we shall be mostly interested in dilute Bose gases, for which $n_{B}a^{3}\ll 1$. Under these circumstances, a numerical solution of the self-consistency equation (128) for a few representative values of $n_{B}a^{3}$ between $10^{-5}$ and $10^{-1}$ yields the values for the quantities $\tilde{\sigma}$ and $C_{d}$ shown in Table 1. The variation of $\tilde{\sigma}$ vs. $n_{B}a^{3}$ over the whole range $10^{-6}\leq n_{B}a^{3}\leq 0.1$ is shown in Fig. 3. As can be seen, $\tilde{\sigma}$ takes values that are close to $0.6$ for $n_{B}a^{3}<10^{-4}$, and decreases to about $0.2$ for $n_{B}a^{3}$ between $0.02$ and $0.1$.

We are now in a position to find the ground state energy $E_{0}=\langle\Psi(N)|\hat{H}|\Psi(N)\rangle$ of the system, which is given by Eq. (108). For simplicity, in Eq. (109b) we shall approximate the interaction potential as $v({\bf r})=g\delta({\bf r})$, which gives us:

[TABLE]

It then follows that Eq. (108) can be written in the form:

[TABLE]

Using the expression (123) of the coefficients $c_{\bf k}$ into this last equation and transforming the sum into an integral, we obtain (in three dimensions):

[TABLE]

The physical interpretation of the above equation is quite simple: it is just the ground state energy derived in Ref. Ettouhami2012, , plus an extra term

[TABLE]

which represents the Fock interaction between depleted bosons.

Numerical evaluation of the integral on the rhs of Eq. (131) for $n_{B}a^{3}=10^{-4}$ yields, for $\tilde{\sigma}=0.568$ and $C_{d}=0.9897$:

[TABLE]

On the other hand, for $n_{B}a^{3}=10^{-3}$, using $\tilde{\sigma}=0.452$ and $C_{d}=0.9672$, we have:

[TABLE]

Finally, for $n_{B}a^{3}=10^{-2}$ we let $\tilde{\sigma}=0.263$ and $C_{d}=0.8990$, and we obtain:

[TABLE]

At this point we will introduce some notation, and write the ground state energy $E_{0}$ at $T=0$ in the form:

[TABLE]

where $C_{E}$ is a positive dimensionless quantity whose precise value depends on the parameter $n_{B}a^{3}$. Just as we discussed in Ref. Ettouhami2012, , we will not introduce any “renormalization” of the “bare” interaction strength $g$ to convert the negative ground state energy $E_{0}$ into a positive quantity. In this author’s view, such “renormalizations” where divergent integrals are substracted from $g$ and where the ground state energy is converted from a negative value to a positive one while leaving the coefficients of the wavefunction intact,FetterWalecka ; LeggettNJP are improper and mathematically unjustified and hence will be avoided in the rest of this paper. The negative $E_{0}$ in Eq. (136) with $C_{E}>0$ is what we will use to derive the thermodynamics of the dilute Bose gas in the canonical ensemble in Sec. IV below. For the reader’s convenience, the three representative values of $C_{E}$ listed in Eq. (133)-(135) are summarized in Table 2. The variation of $C_{E}$ over a wider range of the parameter $n_{B}a^{3}$ is shown in Fig. 5.

V.2 Excitation energies

Having found the ground state energy of the system of interacting bosons in presence of Fock interactions between depleted bosons, we now want to find the excitation spectrum of the system within our variational approach. But, before we do so, we shall digress a little bit and uncover yet another flaw in the standard Bogoliubov theory by finding the excitation energies in presence of Fock interactions between depleted bosons in that theory and showing that these excitation energies diverge in the limit of small momenta. This will be done in the following paragraph.

V.2.1 Digression: excitation energies in presence of Fock interactions between depleted bosons in the standard Bogoliubov formulation

Let us find the amount of energy $E_{\bf k}$ required to create a single excitation of wavevector ${\bf k}$ above the ground state in the standard Bogoliubov theory when Fock interactions between depleted bosons is taken into account. This is the quantity:

[TABLE]

Using the fact that $\hat{H}=\hat{H}_{B}+\hat{H}_{F}^{\prime}$, and the fact that the quantity $\langle\Psi_{B}|\alpha_{\bf k}\hat{H}_{B}\alpha_{\bf k}^{\dagger}|\Psi_{B}\rangle-\langle\Psi_{B}|\hat{H}_{B}|\Psi_{B}\rangle$ is nothing but the Bogoliubov excitation energy $E_{\bf k}^{B}=\sqrt{\varepsilon_{\bf k}(\varepsilon_{\bf k}+2gn_{0})}$, we obtain:

[TABLE]

Now, noticing that:

[TABLE]

and using the commutation relation $\alpha_{\bf k}\alpha_{\bf k}^{\dagger}=1+\alpha_{\bf k}^{\dagger}\alpha_{\bf k}$, we can write:

[TABLE]

For the purpose of calculating the commutator $[\hat{H}_{F}^{\prime},\alpha_{\bf k}^{\dagger}]$, we shall for simplicity use the linearized approximation to $\alpha_{\bf k}$ and $\alpha_{\bf k}^{\dagger}$ given in Eq. (64). Tedious but straightforward algebra leads to the result:

[TABLE]

where $n_{d}=\sum_{\bf k\neq 0}\langle\Psi_{B}|a_{\bf k}^{\dagger}a_{\bf k}|\Psi_{B}\rangle/V$ is the density of depleted bosons. In this last equation, the second term becomes negligibly small compared to the first term in the thermodynamic limit $V\to\infty$, and so we obtain for the excitation energy:

[TABLE]

As can be seen, apart from the familiar term $E_{\bf k}^{B}$ the excitation energy acquires an addtional contribution due to the addition of the Hamiltonian $\hat{H}_{F}^{\prime}$. This additional contribution does not vanish in the limit $k\to 0$, but instead diverges like $1/k$, which is an erroneous feature of the standard, number non-conserving Bogoliubov theory, to be added to the list of other erroneous features we already discussed in Ref. Ettouhami2012, . To remedy this divergence, we shall evaluate the excitation energies directly within the variational approach we developed in this last reference. This will be done in the next paragraph.

V.2.2 Excitation energies in presence of Fock interactions between depleted bosons: variational approach

We will not give here all the steps regarding the derivation of the excitation energies, and refer the interested reader to the calculation presented in Ref. Ettouhami2012, , which remains unchanged, except for the additonal $gn_{d}/2$ term on the rhs of Eq. (116b) or equivalently the additonal $n_{d}/2n_{B}$ on the rhs of Eq. (128), with the result for the excitation energy $E_{\bf k}$ of wavevector ${\bf k}$ given by:

[TABLE]

Due to the fact that $Q^{2}\to\tilde{\sigma}/C_{d}$ as $k\to 0$, we see that the excitation spectrum $E_{\bf k}$ has a finite gap at $k=0$:

[TABLE]

Fig. 6 shows how the gap varies with the parameter $n_{B}a^{3}$. The gap is largest for very small values $n_{B}a^{3}<10^{-3}$, where $E_{k\to 0}\simeq n_{B}v({\bf 0})$, while for $0.01<n_{B}a^{3}<0.1$ the gap is a little smaller, of order $0.8n_{B}v({\bf 0})$.

In anticipation of the customary criticism that the finite gap in Eq. (144) violates Goldstone’s theorem, we here shall repeat a counter-example from Ref. Ettouhami2012, proving that this theorem does not even apply to this kind of microscopic Hamiltonian. Indeed, if Goldstone’s theorem applied to Bogoliubov’s Hamiltonian, it would equally apply to the Hartree-Fock Hamiltonian given by the following expression:

[TABLE]

in which case the excitation spectrum of $\hat{H}_{HF}$ should not have a gap. Unfortunately, this turns out to be not true, as the excitation spectrum of $\hat{H}_{HF}$, which is known exactly,Ettouhami2012 does have a gap. We therefore conclude that Goldstone’s theorem does not apply to the Hartree-Fock Hamiltonian $\hat{H}_{HF}$, and in a similar fashion does not apply to the Hamiltonian $\hat{H}$ studied in this Section. We refer the reader to section 8.4 of Ref. Ettouhami2012, where we discuss questions related to the apparent violation of the Goldstone and Hugenholtz-Pine theorems in quite some detail, and show explicitly that the gap in Eq. (144) does not violate either theorem.

Before we close this Section, we want to attract the reader’s attention to the fact that neither the ground state energy in Eq. (131) nor the excitation spectrum in Eq. (143a) depend on the density of condensed bosons $n_{0}$. Hence the canonical partition function $Z_{C}$ in Eq. (94) will only depend on $(N,V,T)$ as it should, and not also on $N_{0}$. When in Sec. VII we will want to derive the grand partition function $Z_{G}$, we will just have to perform an additional trace over $N$ as in Eq. (100), and the resulting expression of $Z_{G}$ will not depend on $N_{0}$.

With the knowledge of the ground state and excitation energies, we are now equipped to perform traces over eigenfunctions corresponding to these energies and obtain partition functions to study the thermodynamics of interacting bosons. In the following Section, we will study the statistical mechanics of interacting bosons in the canonical formalism using the ground state and excitation energies derived in the current Section. A discussion of the grand-canonical formalism will be presented in Sec. VII.

VI Statistical mechanics of interacting bosons: canonical number-conserving formulation

We now want to discuss the thermodynamics of an interacting Bose gas in the canonical ensemble using the number conserving formalism developed in Ref. Ettouhami2012, which avoid the use of Bogoliubov’s prescription. We will start by investigating the condensate fraction and the shift in transition temperature $T_{c}$ in Subsection VI.1 before discussing the thermodynamic properties in Subsection VI.2.

VI.1 Condensate fraction and shift in transition temperature

In Sec. IV, we already derived the expectation value $\langle\hat{N}_{\bf k}\rangle$ of the number of bosons of wavevector ${\bf k}$ as a canonical trace at fixed number of bosons $N$, and the main steps of that derivation remain valid in the variational approach of Sec. V provided that we replace the coherence factors $u_{\bf k}$ and $v_{\bf k}$ we used in that section by the appropriate expressions for the number-conserving theory:

[TABLE]

where $\bar{v}({\bf k})$ was defined in Eq. (109a). Using Eqs. (146) into Eq. (84), we obtain after a few manipulations:

[TABLE]

where in going from the first to the second line we introduced the dimensionless quantity $Q^{2}$ defined in Eq. (127). On the other hand, it is not difficult to show that the quantity $\beta E_{\bf k}$ inside the hyperbolic cotangent can be written in dimensionless form, where temperature $T$ is measured in units of the critical temperature of a non-interacting gas $T_{c0}$, see Eq. (48), as follows:

[TABLE]

Before going any further, we pause a moment to note the differences between the expression of $N_{\bf k}$ given in Eq. (147a) and the similar expression derived within the number non-conserving approach, Eq. (26b). Apart from the appearance of the quantity $\sigma$ on the rhs of Eq. (147a), we also note that $n_{0}$ does not appear on the rhs of that equation, and instead $n_{0}$ is replaced with the total density of bosons $n_{B}$, which we remind the reader is fixed and does not depend on temperature. We thus see that in the variational number-conserving approach the confusion about whether $n_{0}$ on the rhs of Eq. (26b) depends or not on temperature does not arise because $n_{0}$ does not appear in the expression of $N_{\bf k}$ in the first place, hence confirming our earlier claim in Sec. IV that $n_{0}$ on the rhs of Eq. (26b) should in fact be replaced by $n_{B}$.

Now, using Eq. (148) in Eq. (147b), we obtain after a few manipulations that the condensate fraction in the canonical ensemble using the variational number-conserving approach of Sec. V is given by:

[TABLE]

In Fig. 7 we plot the condensate fraction $n_{0}(T)/n_{B}$ vs. temperature $T$ for $n_{B}a^{3}=10^{-4}$ (upper panel) and $10^{-3}$ (lower panel). The solid curve is the result of the variational number-conserving approach, and the dashed curve is the result obtained from the standard, number non-conserving approximation. It is seen that the curves obtained from the number-conserving approach lie higher than the curves predicted by the number non-conserving approximation. One can see that it takes higher temperatures to deplete the condensate completely in the number-conserving approach, hence this approach leads to higher critical temperatures than what is predicted by the standard Bogoliubov formulation.

Unlike the case of the number non-conserving approximation where we were able to derive an analytical expression for $\Delta T_{c}/T_{c0}$ in terms of $n_{B}a^{3}$ for small values of $n_{B}a^{3}$, in the variational number-conserving approach such a derivation is made difficult by the fact that $\tilde{\sigma}$ depends on $n_{B}a^{3}$, and the fact that this dependence is unknown analytically. Here we shall investigate this dependence numerically, and to this end in the upper panel of Fig. 8 we plot $\Delta T_{c}/T_{c0}$ vs. $n_{B}a^{3}$ in the variational number-conserving approximation (see also Table 3). This plot shows the infinite slope behaviour near the origin which is charateristic of a dependence $\Delta T_{c}/T_{c0}\propto(n_{B}a^{3})^{\eta}$ with $\eta<1$. In the lower panel we plot $\ln(\Delta T_{c}/T_{c0})$ vs. $\ln(n_{B}a^{3})$, and it turns out that this plot can be well approximated by a straight line of slope $0.16357$, which seems to point to a dependence of the form:

[TABLE]

which is reminiscent of the number non-conserving result of Eq. (52). The only difference between the variational and standard approach is that while in the latter the curve of $\Delta T_{c}/T_{c0}$ vs. $n_{B}a^{3}$ has a maximum around $n_{B}a^{3}\approx 0.01$ and decreases monotonically for values $n_{B}a^{3}\geq 0.01$, no such maximum appears in the number-conserving case, and $\Delta T_{c}/T_{c0}$ seems to be following the variation in Eq. (150) all the way across the region $0\leq n_{B}a^{3}\leq 0.1$.

VI.2 Thermodynamics in the canonical ensemble

We now turn our attention to the thermodynamic functions of the system. The result derived in Sec. IV for the canonical partition function $Z_{C}(N,T,V)$, Eq. (69c), remains valid in the number-conserving variational formalism, provided that we use for the excitation energies $E_{\bf k}$ the variational expression from Eq. (143a). It then follows that the Helmholtz free energy

[TABLE]

is given by (we here assume that $N\gg 1$ so that $(N-1)\simeq N$):

[TABLE]

The internal energy of the gas $U=\langle\hat{H}_{B}\rangle$ can be obtained using the thermodynamic identity $U=F+\beta(\partial F/\partial\beta)$, valid in the canonical ensemble, with the result:

[TABLE]

If we call excess internal energy $U_{exc}(N,V,T)$ the temperature-dependent term on the rhs of Eq. (153), i.e.

[TABLE]

then in dimensionless units this quantity is given by:

[TABLE]

The heat capacity at constant volume $C_{v}=(\partial U/\partial T)_{N,V}$ is given by:

[TABLE]

and in dimensionless units can be written in the form:

[TABLE]

where we remind the reader that the quantity $\beta E_{\bf k}$ is given in dimensionless units by the rhs of Eq. (148). In Fig. 9, we show plots of the excess internal energy $U_{exc}(T)$ and specific heat $C_{v}(T)$ as obtained in the canonical formalism based on our variational number-conserving theory, Eqs. (155) and (157). In both plots, the dashed vertical line at the tip of the curve indicates the location of the critical temperature $T_{c}$ for the corresponding value of $n_{B}a^{3}$. As can be seen, at any given value of temperature, the excess internal energy $U_{exc}$ and heat capacity $C_{v}$ decrease as $n_{B}a^{3}$ increases, reflecting the fact that higher repulsive interactions tend to reduce the number of depleted bosons, hence reducing the values of $U_{exc}$ and $C_{v}$.

For completeness, let us examine what the above expressions of $U_{exc}$ and $C_{v}$ become in the naive Bogoliubov theory. In that case, $\tilde{\sigma}=0$, $C_{d}=1$ and $Q^{2}=\tilde{k}^{2}$, and we can write:

[TABLE]

In Fig. 10 we plot the excess internal energy (upper panel) and specific heat (lower panel) in the naive Bogoliubov theory for $n_{B}a^{3}=10^{-4},10^{-3}$, and $10^{-2}$. In both plots, the dashed vertical line at the tip of each curve indicates the location of the critical temperature $T_{c}$ for the corresponding value of $n_{B}a^{3}$, with $T_{c}=1.246T_{c0}$ for $n_{B}a^{3}=10^{-4}$, $T_{c}=1.332T_{c0}$ for $n_{B}a^{3}=10^{-3}$ and $T_{c}=1.396T_{c0}$ for $n_{B}a^{3}=10^{-2}$. These plots are qualitatively similar to the plots obtained from the variational approach shown in Fig. 9. There is a quantitative difference coming from the fact that the critical temperature $T_{c}$ is a little bit higher in the variational approach than in the naive Bogoliubov approach for any given value of the parameter $n_{B}a^{3}$ owing to the fact that when the conservation of the number of bosons is enforced, the ground state depletion becomes extremely small, and hence higher temperatures are needed in order to completely deplete the state ${\bf k}=0$ of all its bosons.

VII Statistical mechanics of interacting bosons: grand-canonical formulation

Let us recapitulate what we have done so far. After having realized that the field-theoretic calculation of the grand partition function based on Bogoliubov’s prescription involved an incomplete trace where the number of bosons in the condensate $N_{0}$ is not traced over, we generalized the variational approach of Ref. Ettouhami2012, to include Fock interactions between depleted bosons, and derived the thermodynamics of a dilute Bose gas in the canonical ensemble within this method. The advantage of using the variational method lies in the fact that this approach allows us to avoid Bogoliubov’s prescription, which on one hand eliminates the spurious temperature dependence introduced by this prescription in the BBP approach, and on the other hand allows us to evaluate traces in the canonical ensemble cleanly, with all the relevant variables being traced over. The result is a canonical partition function $Z_{C}(N,V,T)$ that does not depend on $N_{0}$ at all, but depends on the total number of bosons $N$ instead, as it should. In this section, we want to explore how one can perform the extra trace over $N$ to derive the grand canonical partition function of the system. This will be done in the following Subsection.

VII.1 Using saddle-point techniques to evaluate the grand partition function

In this Subsection, we want to examine the question of how to formulate the grand-canonical description of an interacting Bose gas using the canonical formulation of Bogoliubov’s theory from Sec. V as a starting point. In this Section, we shall use the expression of the ground state energy $E_{0}$ given by the standard Bogoliubov theory, in spite of the grave reservations Ettouhami2012 this author has about the way this expression is usually derivedFetterWalecka and in particular about the sign of the constant $C_{E}$, which in this approach is given by:

[TABLE]

which results in an overall positive correction to the Gross-Pitaevskii result $\frac{1}{2}gn_{B}^{2}$. Our goal behind using the standard Bogoliubov theory in this section is to show for the record how to derive a grand-canonical formulation for the naive Bogoliubov approach. On a more practical level, using Bogoliubov’s theory with constant $C_{E}$ allows us to ignore the complication arising from the variation of $C_{E}$ with the expansion parameter $n_{B}a^{3}$, which in turn allows us to avoid taking derivatives of $C_{E}$ with respect to $n_{B}$ at various steps of the calculation.

We shall start by rewriting the expression of the grand-canonical partition function $Z_{G}(\mu,V,T)$, which is given by:

[TABLE]

where ${Z_{C}}$ is the canonical partition function and $F_{N}$ is the canonical free energy for a system of $N$ bosons, and where in going from the first to the second line use has been made of the fact that ${Z_{C}}=\exp(-\beta F_{N})$. Now, in Sec. V we saw that the canonical free energy $F_{N}$ is given by:

[TABLE]

Dividing both sides by $V$ and transforming the sum to an integral, we obtain:

[TABLE]

where we approximated $(N-1)\simeq N$ and used the shorthand notation $\int_{\bf k}=\int d{\bf k}/(2\pi)^{3}$. Introducing the dimensionless wavevector $\tilde{k}=k/k_{0}$, we obtain:

[TABLE]

where the excitation spectrum $E_{\tilde{k}}$ is given by:

[TABLE]

Now, using Eq. (136) in Eq. (163), we finally obtain:

[TABLE]

where we denote by $f_{1}(x)$ the function:

[TABLE]

A plot of this function is shown in Fig. 11. It is seen that $f_{1}(x)$ takes large negative values for $0<x<1$, and goes to zero negatively for $x>1$.

Replacing the expression of $F_{N}$ from Eq. (165) into Eq. (160b), we find that the grand partition function

[TABLE]

can be written in the form (note that $n_{B}=N/V$ here is not an actual density of bosons but merely a dummy variable which is being summed over):

[TABLE]

Now, the quantity on the rhs can be interpreted as a Riemann sum, which in the limit $V\to\infty$ converges exactly to an integral over the variable $n_{B}=N/V$, with the result:

[TABLE]

where $g(n_{B})$ is the following function:

[TABLE]

The integral in Eq. (169) can be calculated to a very good approximation when $V\to\infty$ using saddle-point methods. To this end, we start by finding the value of $n_{B}$ for which the argument of the exponential is maximal. Setting the derivative $g^{\prime}(n_{B})$ to zero, we find that the value of $n_{B}$ maximizing the exponential satisfies the following equation:

[TABLE]

where we defined:

[TABLE]

For the rest of this section we shall place ourselves in the particular case where the interaction potential is a delta function in real space, $v({\bf r})=g\delta({\bf r})$. Then $v({\bf k})=v({\bf 0})=g$, and the above expression of $f_{2}$ reduces to (note that $f_{2}$ in this case is simply the derivative of $f_{1}$, $f_{2}(x)=f_{1}^{\prime}(x)$):

[TABLE]

A plot of this function is shown in Fig. 12. It is seen that $f_{2}(x)$ takes large positive values for $0<x<1$ and tends to zero positively for $x>1$.

In order to solve Eq. (171) for $n_{B}$ anlytically, it is necessary to make an assumption about the value of the quantity $n_{B}a^{3}$. If we denote by $n_{B}^{*}$ the value of $n_{B}$ that solves Eq. (171), we shall assume that $\sqrt{n_{B}^{*}a^{3}}\ll 1$. This will allow us organize the solution as an expansion in powers of the small parameter $\sqrt{n_{B}^{*}a^{3}}$. Then, to zero-th order in this small parameter, we simply have:

[TABLE]

The next order approximation to the value of $n_{B}^{*}$ is obtained by letting $n_{B}=\mu/v({\bf 0})$ in the $(n_{B}a^{3})^{\frac{1}{2}}$ terms in Eq. (171), with the result:

[TABLE]

As mentioned above, the functions $f_{1}(x)$ and $f_{2}(x)$ assume very large values at small arguments, and become small for values of $x$ that are greater than unity. Therefore, the expansion in powers of $\sqrt{n_{B}^{*}a^{3}}$ in Eq. (175) is only valid for values of $\beta\mu$ such that $\beta\mu>1$, i.e. $k_{B}T<\mu$; otherwise, the expansion (175) is no longer valid, and a numerical method would be needed to solve Eq. (171). Assuming the condition $\beta\mu>1$ to be true, we now want to perform a Taylor expansion of $g(n_{B})$ around $n_{B}^{*}$, writing:

[TABLE]

where in going from the first equality to the second we used the fact that $g^{\prime}(n_{B}^{*})=0$ and we assumed that $g^{\prime\prime}(n_{B}^{*})<0$, i.e. that $n_{B}^{*}$ is a maximum of $g(n_{B})$. The grand-canonical partition function now becomes:

[TABLE]

where in going from the first to the second line we extended the lower limit of integration to $-\infty$ because the exponential decays extremely rapidly away from $n_{B}^{*}$ when $V\to\infty$. At this point the integral can be easily calculated, and we obtain the following result for the grand-partition function $Z_{G}$:

[TABLE]

To lowest order in $(n_{B}^{*}a^{3})^{\frac{1}{2}}$, the second derivative $g^{\prime\prime}(n_{B}^{*})$ is given by:

[TABLE]

and is indeed negative if $v(\bf 0)>0$. On the other hand, upon using the value of $n_{B}^{*}$ from Eq. (175) into the expression of $g(n_{B}^{*})$ we obtain:

[TABLE]

Given this result for the grand-canonical partition function, one can immediately write the expression of the grand-potential $\Omega=-k_{B}T\ln{Z_{G}}$, which is given by:

[TABLE]

The important thing to note in this result is that the rhs only depends on the thermodynamic variables $(\mu,T,V)$, as should be the case in a truly grand-canonical formulation:

[TABLE]

In particular, previous formulations almost invariably express $\Omega$ as a function of $\mu$, $T$, $V$, but also of $n_{0}(T)$ which, in an incomplete tracing procedure is never traced over (sometimes, both $n_{0}$ and the density of bosons $n_{B}$, which is considered to be temperature-independent, appear in the expression of $\Omega$). Here on the contrary, both $n_{0}$ and $n_{B}$ have been traced over and do not appear on the rhs of the above equation for $\Omega(\mu,T,V)$. Moreover, we are going to show that $n_{B}$ itself depends on temperature, and cannot be taken to be a temperature-independent quantity as in previous formulations.

VII.2 Finding the average number of bosons \mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}\hat{N}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}} in the system

Let us now find the value of the average number of bosons \mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}\hat{N}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}} for given values of $\mu$, $T$ and $V$. In the rest of this section, we will use double angular brackets to denote the grand-canonical average. We do so in order to emphasize that this average consists of a double trace procedure, by contrast with the canonical average which involves a single trace. We have:

[TABLE]

Performing the calculation, we find after a few steps:

[TABLE]

Note that this value of \mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}\hat{N}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}} agrees with the value $n_{B}^{*}$ of $n_{B}$ which maximizes the argument $g(n_{B})$ of the exponential, see Eq. (175). More importantly, the expression on the rhs of Eq. (185) shows that \mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}\hat{N}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}} varies with the thermodynamic variables $\mu$, $T$ and $V$:

[TABLE]

as it should in a grand-canonical description where the system is in contact with a particle reservoir and the average number of particles in the system varies, in particular, when temperature is varied. As $T\to 0$, using the fact that $\lim_{x\to\infty}f_{1}(x)=\lim_{x\to\infty}f_{2}(x)=0$, we obtain that the average number of bosons \mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}\hat{N}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}} approaches the limiting value:

[TABLE]

In order to be able to draw plots, we now need to introduce dimensionless units. By analogy with the canonical case where temperatures were measured in units of the critical temperature of an ideal Bose gas $T_{c0}$, we here shall introduce the characteristic temperature $T_{0}$ such that:

[TABLE]

We shall start by writing:

[TABLE]

Then, using the definition of $T_{0}$, Eq. (188), and the fact that $v({\bf 0})=g=4\pi a\hbar^{2}/m$, we obtain after a few manipulations:

[TABLE]

On the other hand, going back to Eq. (189), we find that we can write $\beta\mu$ in dimensionless units in the form:

[TABLE]

We now can write the dimensionless quantity $a^{3}g(n_{B})$ in the form:

[TABLE]

In Fig. 13, we plot the dimensionless function $a^{3}g(n_{B})$ for $(\mu a^{3}/v({\bf 0}))=10^{-4}$ and $T/T_{0}=0.2$, where $v({\bf 0})/(k_{B}Ta^{3})$ on the rhs of Eq. (192) is expressed in dimensionless units using Eq. (190). It is seen that $g(n_{B})$ has a pronounced maximum around $n_{B}a^{3}=10^{-4}$, in agreement with the zeroth-order approximation of Eq. (174). On the other hand, in Fig. 14 we plot the ratio \mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}\hat{N}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}}/\big{[}\mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}\hat{N}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}}|_{T=0}\big{]} as obtained by taking the ratio of Eqs. (185) and (187) for a couple of values of the dimensionless parameter $(\mu a^{3}/v({\bf 0}))$, where on the rhs of Eq. (185) the dimensionless expression of $\beta\mu$ given in Eq. (191) has been used. The two curves in the figure correspond to $(\mu a^{3}/v({\bf 0}))=10^{-3}$ and $10^{-4}$ from top to bottom. The average density of bosons decreases quite substantially with temperature for both values of the parameter $(\mu a^{3}/v({\bf 0}))$, and appear to go to zero for temperatures in the range of $T_{0}$. This shows that the variation of \mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}\hat{N}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}} with temperature is not a small perturbative effect that can be neglected, and has to be taken into account in any grand canonical formulation of the statistical mechanics of an interacting Bose gas.

VII.3 Condensate fraction in the grand-canonical ensemble

Having found the average number of bosons \mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}\hat{N}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}}(\mu,V,T), we now want to study the condensate fraction

[TABLE]

and in particular how this quantity varies with temperature $T$. In order to do so, we shall start by deriving the average number of depleted bosons in the single-particle state of momentum ${\bf k}$, \mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}\hat{N}_{\bf k}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}}=\mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}a_{\bf k}^{\dagger}a_{\bf k}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}}. We have:

[TABLE]

In Eq. (194a), we denote by $\mbox{Tr}_{G}$ the grand-canonical trace: it consists of a trace over all eigenvectors of the Hamiltonian $\hat{H}$ at a given number of bosons $N$, plus a trace over all values of $N$ from [math] to $+\infty$. The fact that the grand-canonical trace consists of two separate traces is explicitly spelled out in Eq. (194b), where we separated out the canonical trace, denoted by $\mbox{Tr}_{N}$ (the index $N$ emphasizing that the trace is performed at a fixed number of bosons), and the summation over all values of the total number of bosons $N$. Now, since $\hat{H}$ commutes with $\hat{N}$, $\exp[-\beta(\hat{H}-\mu\hat{N})]=\exp(\beta\mu\hat{N})\exp(-\beta\hat{H})$, and we can write:

[TABLE]

where in the last equation we used the definition of the average number of bosons in the single-particle state of momentum ${\bf k}$ in the canonical ensemble (we again use the subscript $N$ to emphasize that the average is taken at fixed number of bosons $N$):

[TABLE]

Now, in Eq. (195) we use the fact ${Z_{C}}=\exp(-\beta F_{N})$ to write:

[TABLE]

Transforming the sum into an integral over the dummy variable $n_{B}=N/V$, we can write:

[TABLE]

Using again the Taylor expansion of $g(n_{B})$ around $n_{B}^{*}$, Eq. (176b), we obtain (note that we are extending the lower limit of integration to $-\infty$):

[TABLE]

Now, in the limit $V\to\infty$, the exponential is so sharply peaked around $n_{B}^{*}$ that a very good approximation to the integral can be obtained by using the value $n_{B}^{*}$ that maximizes $g(n_{B})$ in the expression of $\langle a_{\bf k}^{\dagger}a_{\bf k}\rangle_{N}$. Doing that, and taking this quantity outside the integral, we find:

[TABLE]

Given the expression of $Z_{G}$ found in Eq. (179), we finally can write (note that the single angular brackets average on the rhs of Eq. (201a) is a canonical average taken at $n_{B}=n_{B}^{*}$):

[TABLE]

We therefore see that in the grand-canonical ensemble the average number of bosons in the single-particle state of momentum $\bf k$ is given by the same expression as in the canonical ensemble, Eq. (85), except that now the total density of bosons $n_{B}$ is replaced by the average density $n_{B}^{*}(\mu,V,T)$.

Taking the sum over all vectors ${\bf k}$, we find that the condensate fraction in the grand-canonical ensemble is given by (note that this result is identical to the one in Eq. (54), except for the fact that $n_{B}$ in this last equation is now replaced with $n_{B}^{*}(T)$):

[TABLE]

Rearraging terms, the above expression can be written in the form:

[TABLE]

A plot of the condensed fraction \mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}n_{0}(T)\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}}/\mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}n_{B}(0)\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}} vs. temperature $T$ in the grand canonical ensemble is shown in Fig. 15, with the solid curve corresponding to $(\mu a^{3}/v({\bf 0}))=10^{-4}$ and the dashed curve to $(\mu a^{3}/v({\bf 0}))=10^{-3}$. We again note that higher values of $(\mu a^{3}/v({\bf 0}))$ lead to higher values of the critical temperature $T_{c}$, which is qualitatively the same behaviour of $T_{c}(a)$ we found previously in the canonical ensemble.

This concludes our discussion of the grand-canonical formulation of the statistical mechanics of interacting bosons. Our take-away from this section is that any such formulation should necessarily take into account the fact that the total number of bosons in the system is not fixed, but is a function of the grand-canonical variables $(\mu,V,T)$.

VIII Discussion and Conclusions

In this paper, we have discussed many aspects of the statistical mechanics of a dilute gas of interacting bosons. In doing so, we have covered a lot of ground and brought forth a number of new ideas, and so we will offer the following thoughts by way of conclusion:

The standard BBP formulation of the statistical mechanics of interacting bosons performs an incomplete trace over the occupation numbers $N_{i}$ of single particle states of momentum ${\bf k}_{i}$, where the occupation number $N_{0}$ of the condensate is not traced over. The result of this partial tracing procedure is equivalent to a canonical trace at fixed boson number $N$. The realization that previous statistical treatments of interacting Bose gases are canonical in nature and not grand-canonical as commonly thought is the main result of this paper, a result that sets the record straight on several decades of research on the statistical mechanics of interacting bosons using BBP type of methods. 2. 2.

Another important result of this paper is the realization that the discontinuous jump in the condensate density $n_{0}(T)$ at the critical transition temperature $T_{c}$ that is commonly obtained by these BBP type of theories originates from the inappropriate generalization of the Bogoliubov prescription $\Psi({\bf r},\tau)\simeq\sqrt{n_{0}}+\psi({\bf r},\tau)$ to finite temperatures. Indeed, not only is the ${\bf k}=0$ Fourier component of $\Psi({\bf r},\tau)$, which is a flucuating field in $\tau$-space, replaced with a scalar $\sqrt{n_{0}}$, but furthermore this scalar is erroneously assumed to have a temperature dependence even though no thermal averaging process is involved. This state of great confusion that is directly caused by Bogoliubov’s prescription undescores the importance of using methods that circumvent this type of ad-hoc approximations. Within the number-conserving variational approach used in this paper, where Bogoliubov’s prescription is avoided altogether, no such issues arise, as the condensate fraction $n_{0}(T)/n_{B}$ vanishes continuously at $T=T_{c}$ and the jump discontinuity discussed above does not occur. 3. 3.

A hallmark of grand-canonical formulations is that the total number of particles in the system is not constant, but depends on the thermodynamic variables $(\mu,V,T)$. To the best of the author’s knowledge, none of the previous formulations of the statistical mechanics of interacting bosons made any attempt (let alone succeeded) at reproducing this important property, as in all these theories the density of bosons in the system $n_{B}$ is commonly taken to be a constant. A prominent example of this is the common use of the relation $n_{B}=n_{0}(T)+n_{1}(T)$ with a temperature-independent $n_{B}$ (most notably to find the condensate fraction $n_{0}(T)/n_{B}$), which is only correct to do in the canonical ensemble. In Sec. VII, we have shown how one can devise a formulation which is truly grand-canonical by tracing over the number of bosons in the system using saddle-point techniques. When this is done, the average number of bosons \mathopen{\hbox{{\langle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\langle}}}\hat{N}\mathclose{\hbox{{\rangle}}\mkern 2.0mu\kern-3.49998pt\leavevmode\hbox{{\rangle}}} that is obtained does vary with temperature, and we have shown that this variation is not small, and may therefore not be neglected. 4. 4.

Since previous field-theoretic formulations of the BBP approach were essentially canonical in nature, and not grand-canonical in the least, the validity of many higher order perturbation theory results Shi1998 for the spectral density function, damping effects and the lifetime of quasiparticles becomes highly questionable. On one hand, the perturbative methods used to derive these results assume the underlying statistical ensemble to be grand-canonical, while these studies are essentially canonical. On the other hand, the Bogoliubov prescription $\Psi({\bf r},\tau)\simeq\sqrt{n_{0}(T)}+\psi({\bf r},\tau)$ introduces a spurious temperature dependence in the expression of the normal and anomalous Green’s functions $G_{11}({\bf k},\omega)$ and $G_{12}({\bf k},\omega)$ making any results derived within these theories, one is sorry to say, of highly dubious character. To be clear, this is not a statement that this author is taking lightly, nor is there any derogatory or polemical intent in it. The fact is, if we can all agree that Bogoliubov’s prescription introduces an erroneous temperature dependence in the expression of $n_{0}$ on the rhs of this prescription (which at the most should be regarded as $n_{0}(T=0)$ and not $n_{0}(T)$), and that this erroneous temperature dependence trickles down to the expression of the Green’s functions, then one is automatically led to the conclusion that any results based on such erroneous Green’s functions cannot be taken at face value and need to be thoroughly and carefully re-examined. 5. 5.

In addition to the four items above, which constitute the most important results of this paper, we have derived other, auxiliary rsults along the way which, while not as important, should not be completely overlooked. An example of this would be our prediction that the shift in critical temperature $\Delta T_{c}$ due to the repulsive interactions between bosons varies with the dimensionless parameter $n_{B}a^{3}$ like $\Delta T_{c}\propto(n_{B}a^{3})^{\eta}$, with the exponent $\eta=1/6$ appearing to be quite robust, as it is obtained both in the naive Bogoliubov theory and in our number-conserving variational approximation. 6. 6.

Finally, the variation of the “gap paramter” $\tilde{\sigma}$ and the prefactor $C_{E}$ in the expression of the ground state energy of the dilute Bose gas with the interaction parameter $(n_{B}a^{3})$ are also interesting results in their own right. Having been derived somewhat hastily here, these results deserve a more thorough examination, perhaps in a future contribution. The author understands that the variational theory studied in Ref. Ettouhami2012, and in the present paper, having two variational parameters to be determined self-consistently, is definitely more cumbersome that the sleek-looking and elegant Bogoliubov method. However, as we mentioned in the introduction, this should not be an excuse for us to continue using Bogoliubov’s theory totally unabated, after having witnessed all the shortcomings of this theory whether it be at $T=0$ or at finite temperatures, just because this theory predicts a gapless spectrum. If anything, this study should give us a fresh impetus to look for a more satisfying theory of Bose systems: one that rigorously conserves the number of bosons, and at the same time produces a gapless excitation spectrum. But then again, such a theory may well turn out to be cumbersome and not be as sleek and elegant as the number non-conserving Bogoliubov theory we currently have.

An important lesson of this paper is that there is indeniable value in non field-theoretic, good old first principles calculations. As shown in Sec. IV, it is only when we “looked under the hood,” in a figurative sense, and analyzed the thermal averaging process using a back to basics, no-nonsense approach based on tracing over actual eigenvalues of the Bogoliubov Hamiltonian that we discovered the very important fact that field-theoretic formulations of the statistical mechanics of Bose systems correspond to an incomplete tracing procedure and are canonical in nature, not grand-canonical as they purport to be. Hence, first principles calculations based on actual eigenvalues and eigenfunctions of the Hamiltonian should be used, whenever it is possible, as a sanity check to convoluted field-theoretic calculations where the heavy formalism, sophisticated though it is, can conceal seemingly insignificant details which have great physical importance. It is to be hoped that the ideas presented in this study will help correct common misconceptions about the theories of interacting bosons, hence helping us formulate improved theories for these systems in the future.

Appendix A Expectation value of the Hamiltonian in the number-conserving approach

In this Appendix, we briefly describe how we calculate the expectation value of the total Hamiltonian $\hat{H}=\sum_{\bf k_{j}\neq 0}\hat{H}_{{\bf k}_{j}}$ where

[TABLE]

in the variational ground state $|\Psi(N)\rangle$ of Eq. (105), namely (note that $Z$ here is a normalization constant, chosen so that $\langle\Psi(N)|\Psi(N)\rangle=1$, not a partition function):

[TABLE]

If we denote by $\hat{H}_{\bf k_{j}}^{(0)}$ the part consisting of the first six terms of $\hat{H}_{\bf k_{j}}$ and $\hat{H}_{\bf k_{j}}^{(1)}$ the part consisting of the last term on the rhs of Eq. (204a), then the expectation value $\langle\Psi(N)|\hat{H}_{{\bf k}_{j}}^{(0)}|\Psi(N)\rangle$ was already calculated in Appendix B of Ref. Ettouhami2012, , with the result:

[TABLE]

Let us now calculate the expectation value $\langle\Psi(N)|\hat{H}_{\bf k_{j}}^{(1)}|\Psi(N)\rangle$. We have:

[TABLE]

Hence:

[TABLE]

Using the above result, we obtain that the expectation value of the Hamiltonian $\hat{H}_{\bf k_{j}}^{(1)}$ is given by:

[TABLE]

Rearranging terms, we obtain for the total Hamiltonian $\hat{H}=\sum_{j=1}^{M}\hat{H}_{{\bf k}_{j}}$ the following expression:

[TABLE]

This is Eq. (108) from the text, with $\bar{v}({\bf k})$ and $\tilde{v}({\bf k})$ defined as in Eqs. (109a) and (109b), respectively.

Taking the partial derivative of $\langle\hat{H}\rangle=\langle\Psi(N)|\hat{H}|\Psi(N)\rangle$ with respect to $c_{{\bf k}_{j}}$, using steps that are similar to those in Appendix C of Ref. Ettouhami2012, , we obtain:

[TABLE]

where $\tilde{\cal E}_{{\bf k}_{j}}$ is given by:

[TABLE]

where $\sigma_{\bf k_{j}}$ is defined in Eq. (111b). Now, if we set the rhs of Eq. (211) to zero, we obtain:

[TABLE]

which is nothing but Eq. (110) of the text.