# The Dilute Bose Gas Revised

## Abstract

The well–known results concerning a dilute Bose gas with the short–range repulsive interaction should be reconsidered due to a thermodynamic inconsistency of the method being basic to much of the present understanding of this subject. The aim of our paper is to propose a new way of treating the dilute Bose gas with an arbitrary strong interaction. Using the reduced density matrix of the second order and a variational procedure, this way allows us to escape the inconsistency mentioned and operate with singular potentials of the Lennard–Jones type. The derived expansion of the condensate depletion in powers of the boson density $`n=N/V`$ reproduces the familiar result, while the expansion for the mean energy per particle is of the new form:$`\epsilon =2\pi \mathrm{}^2an/m\{1+128/(15\sqrt{\pi })\sqrt{na^3}(15b/8a)+\mathrm{}\}`$, where $`a`$ is the scattering length and $`b0`$ stands for one more characteristic length depending on the shape of the interaction potential (in particular, for the hard spheres $`a=b`$). All the consideration concerns the zero temperature.

It is well–known that to investigate a dilute Bose gas of particles with an arbitrary strong repulsion (the strong–coupling regime), one should go beyond the Bogoliubov approach  (weak–coupling case) and treat the short–range boson correlations in a more accurate way. An ordinary manner of doing so is the use of the Bogoliubov model with the “dressed”, or effective, interaction potential containing “information” on the short–range boson correlations (see Ref. ). Below it is demonstrated that this manner leads to a loss of the thermodynamic consistency. To overcome this trouble, we propose a new way of investigating the strong–coupling regime which concerns the reduced density matrix of the second order (the 2–matrix) and is based on the variational method.

The 2–matrix for the many–body system of spinless bosons can be represented as : $`\rho _2(𝐫_1^{},𝐫_2^{};𝐫_1,𝐫_2)=F_2(𝐫_1,𝐫_2;𝐫_1^{},𝐫_2^{})/\{N(N1)\},`$ where the pair correlation function is given by

$$F_2(𝐫_1,𝐫_2;𝐫_1^{},𝐫_2^{})=\psi ^{}(𝐫_1)\psi ^{}(𝐫_2)\psi (𝐫_2^{})\psi (𝐫_1^{}).$$
(1)   
Here $`\psi (𝐫)`$ and $`\psi ^{}(𝐫)`$ denote the boson field operators. Recently it has been found  that for the uniform system with a small depletion of the zero–momentum state the correlation function (1) can be written in the thermodynamic limit as follows :

$`F_2(𝐫_1,𝐫_2;𝐫_1^{},𝐫_2^{})=n_0^2\phi ^{}(r)\phi (r^{})`$ (2)   
$`+2n_0{\displaystyle \frac{d^3q}{(2\pi )^3}n_q\phi _{𝐪/2}^{}(𝐫)\phi _{𝐪/2}(𝐫^{})\mathrm{exp}\{i𝐪(𝐑^{}𝐑)\}},`$ (3)   
where $`𝐫=𝐫_1𝐫_2,𝐑=(𝐫_1+𝐫_2)/2`$ and similar relations take place for $`𝐫^{}`$ and $`𝐑^{}`$, respectively. In Eq. (3) $`n_0=N_0/V`$ is the density of the particles in the zero–momentum state, $`n_q=a_𝐪^{}a_𝐪`$ stands for the distribution of the noncondensed bosons over momenta. Besides, $`\phi (r)`$ is the wave function of a pair of particles being both condensed. In turn, $`\phi _{𝐪/2}(𝐫)`$ denotes the wave function of the relative motion in a pair of bosons with the total momentum $`\mathrm{}𝐪`$, this pair including one condensed and one noncondensed particles. So, Eq. (3) takes into account the condensate–condensate and supracondensate–condensate pair states and is related to the situation of a small depletion of the zero–momentum one–boson state. For the wave functions $`\phi (r)`$ and $`\phi _𝐩(𝐫)`$ we have

$`\phi (r)=1+\psi (r),\phi _𝐩(𝐫)=\sqrt{2}\mathrm{cos}(\mathrm{𝐩𝐫})+\psi _𝐩(𝐫)(p0)`$ (4)   
with the boundary conditions $`\psi (r)0`$ and $`\psi _𝐩(𝐫)0`$ for $`r\mathrm{}.`$ The functions $`\psi (r)`$ and $`\psi _𝐩(𝐫)`$ can explicitly be expressed in terms of the Bose operators $`a_𝐩^{}`$ and $`a_𝐩`$ . In particular,

$$\stackrel{~}{\psi }(k)=\psi (r)\mathrm{exp}(i\mathrm{𝐤𝐫})d^3r=a_𝐤a_𝐤/n_0.$$
(5)   
Having in our disposal the distribution function $`n_k`$ and the set of the pair wave functions $`\phi (r)`$ and $`\phi _𝐩(𝐫)`$, we are able to calculate the main thermodynamic quantities of the system of interest. In particular, the mean energy per particle is expressed in terms of $`n_k`$ and $`g(r)`$ via the well–known formula

$`\epsilon ={\displaystyle \frac{d^3k}{(2\pi )^3}T_k\frac{n_k}{n}}+{\displaystyle \frac{n}{2}}{\displaystyle g(r)\mathrm{\Phi }(r)d^3r},`$ (6)   
where $`T_k=\mathrm{}^2k^2/2m`$ is the one–particle kinetic energy, $`n=N/V`$ stands for the boson density and the relation

$$g(r)=F_2(𝐫_1,𝐫_2;𝐫_1,𝐫_2)/n^2.$$
(7)   
is valid for the pair distribution function $`g(r)`$.

The starting point of our investigation is the weak–coupling regime which implies weak spatial correlations of particles and, thus, is characterized by the set of the inequalities

$$|\psi (r)|1,|\psi _𝐩(𝐫)|1.$$
(8)   
Specifically, the Bogoliubov model corresponds to the choice 

$$|\psi (r)|1,\psi _𝐩(𝐫)=0.$$
(9)   
Besides, owing to a small depletion of the Bose condensate $`(nn_0)/n`$ we have for the one–particle density matrix $`F_1(r)=\psi ^{}(𝐫_1)\psi (𝐫_2)`$:

$$\left|\frac{F_1(r)}{n}\right|=\left|\frac{d^3k}{(2\pi )^3}\frac{n_k}{n}\mathrm{exp}(i\mathrm{𝐤𝐫})\right|\frac{nn_0}{n}1.$$

So, investigating the Bose gas within the Bogoliubov scheme, we have two small quantities: $`\psi (r)`$ and $`F_1(r)/n`$. This enables us to write Eq. (7) with the help of (3) as follows:

$$g(r)=1+2\psi (r)+\frac{2}{n}\frac{d^3k}{(2\pi )^3}n_k\mathrm{exp}(i\mathrm{𝐤𝐫}),$$
(10)   
where we restricted ourselves to the terms linear in $`\psi (r)`$ and $`F_1(r)/n`$ and put $`\psi ^{}(r)=\psi (r)`$ because the pair wave functions can be chosen as real quantities. Equations for $`\stackrel{~}{\psi }(k)`$ and $`n_k`$ can be found varying the mean energy (6) with (10) taken into account. However, previously one should realize an important point, namely: $`n_k`$ and $`\stackrel{~}{\psi }(k)`$ can not be independent variables. Indeed, when there is no interaction between particles, there are no spatial particle correlations either. So, $`\stackrel{~}{\psi }(k)=0`$ and, since the zero–temperature case is considered, all the bosons are condensed, $`n_k=0`$. While “switching on” the interaction results in appearing the spatial correlations and condensate depletion: $`\stackrel{~}{\psi }(k)0`$ together with $`n_k0`$. In the framework of the Bogoliubov scheme $`\stackrel{~}{\psi }(k)`$ is related to $`n_k`$ by the expression

$$n_k(n_k+1)=n_0^2\stackrel{~}{\psi }^2(k).$$
(11)   
Indeed, the canonical Bogoliubov transformation  implies that

$$a_𝐤=u_k\alpha _𝐤+v_k\alpha _𝐤^{},a_𝐤^{}=u_k\alpha _𝐤^{}+v_k\alpha _𝐤,$$
(12)   
where

$$u_k^2v_k^2=1.$$
(13)   
At zero temperature $`\alpha _𝐤^{}\alpha _𝐤=0`$ and, using Eqs. (5) and (12) we arrive at

$$n_k=v_k^2,\stackrel{~}{\psi }(k)=u_kv_k/n_0.$$
(14)   
With Eqs. (13) and (14) one can readily obtain (11).

Now, let us show that all the results on the thermodynamics of a weak–coupling Bose gas can be derived for the Bogoliubov scheme with variation of the mean energy (6) under the conditions (10) and (11). Inserting (10) into (6) and, then, varying the obtained expression, we arrive at

$$\delta \epsilon =\frac{d^3k}{(2\pi )^3}\left[\left(T_k+n\stackrel{~}{\mathrm{\Phi }}(k)\right)\frac{\delta n_k}{n}+n\stackrel{~}{\mathrm{\Phi }}(k)\delta \stackrel{~}{\psi }(k)\right].$$
(15)   
Relation (11) connecting $`\stackrel{~}{\psi }(k)`$ with $`n_k`$ results in

$$\delta \stackrel{~}{\psi }(k)=\frac{(2n_k+1)\delta n_k}{2n_0^2\stackrel{~}{\psi }(k)}+\frac{\stackrel{~}{\psi }(k)}{n_0}\frac{d^3q}{(2\pi )^3}\delta n_q,$$
(16)   
where the equality

$$n=n_0+\frac{d^3k}{(2\pi )^3}n_k$$
(17)   
is taken into consideration. Setting $`\delta \epsilon =0`$ and using Eqs. (15) and (16), we derive the following expression:

$`2T_k`$ $`\stackrel{~}{\psi }(k)={\displaystyle \frac{n^2}{n_0^2}}\stackrel{~}{\mathrm{\Phi }}(k)(1+2n_k)`$ (19)   
$`+2n\stackrel{~}{\psi }(k)\left(\stackrel{~}{\mathrm{\Phi }}(k)+{\displaystyle \frac{n}{n_0}}{\displaystyle \frac{d^3q}{(2\pi )^3}\stackrel{~}{\mathrm{\Phi }}(q)\stackrel{~}{\psi }(q)}\right).`$   
Here one should realize that Eq. (19) is able to yield results being accurate only to the leading order in $`(nn_0)/n`$ because the used expression for $`g(r)`$ given by (10) is valid to the next–to–leading order . So, Eq. (19) should be rewritten as

$$2T_k\stackrel{~}{\psi }(k)=\stackrel{~}{\mathrm{\Phi }}(k)(1+2n_k)+2n\stackrel{~}{\psi }(k)\mathrm{\Phi }(k).$$
(20)   
Eq. (20) is an equation of the Bethe–Goldstone type or, in other words, the in–medium Schrödinger equation for the pair wave function. As $`2\stackrel{~}{\mathrm{\Phi }}(k)(n_k+n\stackrel{~}{\psi }(k))`$ is the product of the Fourier transforms of $`\mathrm{\Phi }(r)`$ and $`n(g(r)1)`$, we can rewrite Eq. (20) in the more customary form

$$\frac{\mathrm{}^2}{m}^2\phi (r)=\mathrm{\Phi }(r)+n\mathrm{\Phi }(|𝐫𝐲|)\left(g(y)1\right)d^3y.$$
(21)   
The structure of Eq. (21) is discussed in the papers . Here we only remark that the right–hand side (r.h.s.) of (21) is the in–medium potential of the boson–boson interaction in the weak–coupling approximation. The system of equations (11) and (20) can easily be solved, which leads to the familiar results :

$`n_k={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{T_k+n\stackrel{~}{\mathrm{\Phi }}(k)}{\sqrt{T_k^2+2nT_k\stackrel{~}{\mathrm{\Phi }}(k)}}}1\right),`$ (22)   
$`\stackrel{~}{\psi }(k)={\displaystyle \frac{\stackrel{~}{\mathrm{\Phi }}(k)}{2\sqrt{T_k^2+2nT_k\stackrel{~}{\mathrm{\Phi }}(k)}}}.`$ (23)   
Now we are able to demonstrate that the investigation of the strong–coupling case based on the Bogoliubov model with the effective boson–boson interaction, results in a loss of the thermodynamic consistency. Indeed, as it was shown in the previous paragraph, any calculating scheme using the basic relations of the Bogoliubov model (10), (11) conclusively leads to Eqs. (20)-(23) provided this scheme does yield the minimum of the mean energy. In this case Eqs. (20)-(23) certainly includes the quantity $`\mathrm{\Phi }(r)`$ which is the “bare” interaction potential appearing in (6). The use of the Bogoliubov model with the effective interaction potential substituted for $`\mathrm{\Phi }(r)`$ can in no way disturb the relations given by (10) and (11). And Eq. (6) is the same in both the weak– and strong–coupling regimes. Thus, any attempts of replacing $`\mathrm{\Phi }(r)`$ by the effective “dressed” potential without modifications of (10) and (11) results in a calculating procedure which does not really provide the minimum of the mean energy. It is nothing else but a loss of the thermodynamic consistency. Remark that we do not mean, of course, that the t–matrix approach or the pseudopotential method can not be applied in the quantum scattering problem. It is only stated that the usual way of combining the ladder diagrams with the random phase approximation faces the trouble mentioned above. Though our present investigation is limited by the consideration of the Bose systems, the derived result gives a hint that the similar situation is likely to take place in the Fermi case, too. In this connection it is worth noting the problem associated with the lack of self–consistency of the standard method of treating the dilute Fermi gas .

The strong–coupling regime is characterized by significant spatial correlations. So, Eq. (9) resulting in (10) is not relevant for an arbitrary strong repulsion between bosons at small separations when we have $`\psi (0)=1,\psi _𝐩(0)=\sqrt{2}`$ (see Refs. ). Therefore, to investigate the strong–coupling regime, Eq. (10) should be abandoned in favor of (3). Expression (3) is accurate to the next–to–leading order in $`(nn_0)/n`$. So, using (3) and (7), we can write

$$g(r)=\phi ^2(r)+\frac{2}{n}\frac{d^3q}{(2\pi )^3}n_q\left(\phi _{𝐪/2}^2(𝐫)\phi ^2(r)\right).$$
(24)   
Let us now perturb $`\stackrel{~}{\psi }(k)`$ and $`n(k)`$. Working to the first order in the perturbation and keeping in mind conditions (11) and (24), from (6) we derive:

$$2T_k\stackrel{~}{\psi }(k)=\stackrel{~}{U}(k)(1+2n_k)+2n\stackrel{~}{\psi }(k)\stackrel{~}{U}^{}(k)$$
(25)   
with

$$\stackrel{~}{U}(k)=\phi (r)\mathrm{\Phi }(r)\mathrm{exp}(i\mathrm{𝐤𝐫})d^3r$$
(26)   
and

$`\stackrel{~}{U}^{}(k)={\displaystyle \left(\phi _{𝐤/2}^2(𝐫)\phi ^2(r)\right)\mathrm{\Phi }(r)d^3r}.`$ (27)   
Using Eqs. (26), (27) as well as the relation $`\psi _𝐤(𝐫)\sqrt{2}\psi (r)(k0)`$ (see the boundary conditions (4)), we obtain $`\stackrel{~}{U}(0)\stackrel{~}{U}^{}(0).`$ This implies that the system of Eqs. (11) and (25) is not able to yield the relation $`n_k1/k(k0)`$ following from the “$`1/k^2`$” theorem of Bogoliubov for the zero temperature . Indeed, let us assume $`n_k\mathrm{}`$ for $`k0.`$ Then, from Eq. (11) at $`n=n_0`$ we find $`n|\stackrel{~}{\psi }(k)|/n_k1`$ when $`k0.`$ On the contrary, Eq. (25) gives $`n|\stackrel{~}{\psi }(k)|/n_k\stackrel{~}{U}(0)/\stackrel{~}{U}^{}(0)1`$ for $`k0.`$ So, consideration of the Bose gas based on Eqs. (3) and (11) does not produce satisfactory results. Nevertheless, it is worth noting that Eq. (25) has an important peculiarity which differentiate it from Eq. (20) in an advantageous way. The point is that in both the limits $`n0`$ and $`k\mathrm{}`$ Eq. (25) is reduced to

$$\frac{\mathrm{}^2}{m}^2\phi (r)+\mathrm{\Phi }(r)\phi (r)=0.$$
(28)   
As it is seen, this is the exact “bare” (not in–medium) Schrödinger equation, other than its Born approximation following from (21). Thus, we can expect the line of our investigation to be right.

As it was shown in the previous paragraph, an approach adequate for a dilute Bose gas with an arbitrary strong interaction can not be constructed without modifications of Eq. (11). This is also in agreement with a consequence of the relation

$$|a_𝐤a_𝐤|^2a_𝐤a_𝐤^{}a_𝐤^{}a_𝐤$$
(29)   
resulting from the inequality of Cauchy–Schwarz–Bogoliubov  $`|\widehat{A}\widehat{B}|^2\widehat{A}\widehat{A}^{}\widehat{B}^{}\widehat{B}.`$ With (5) and (29) one can easily derive $`n_0^2\stackrel{~}{\psi }^2(k)n_k(n_k+1)`$. Thus, it is reasonable to assume that Eq. (11) takes into account only the condensate–condensate channel and ignores the supracondensate–condensate ones. Now the question arises how to find corrections to the r.h.s. of Eq. (11). At present we have no regular procedure allowing us to do this in any order of $`(nn_0)/n`$. However, there exists an argument which makes it possible to realize the first step in this direction. The matter is that the alterations needed have to produce the equation for $`\stackrel{~}{\psi }_𝐩(𝐤)`$ which is reduced to the equation for $`\stackrel{~}{\psi }(k)`$ in the limit $`p0.`$ Though this requirement does not uniquely determine the corrections to Eq. (11), it turns out to be significantly restrictive. In particular, even the simplest variant of correcting Eq. (11) in this way, leads to promising results. Indeed, this variant is specified by the expression

$`n_k(n_k+1)=n_0^2\stackrel{~}{\psi }^2(k)+2n_0{\displaystyle \frac{d^3q}{(2\pi )^3}n_q\stackrel{~}{\psi }_{𝐪/2}^2(𝐤)}.`$ (30)   
Eq. (30) is valid to the next–to–leading order in $`(nn_0)/n`$. So, we may rewrite it as

$$n_k(n_k+1)=n^2\stackrel{~}{\psi }^2(k)+2n\frac{d^3q}{(2\pi )^3}n_q(\stackrel{~}{\psi }_{𝐪/2}^2(𝐤)\stackrel{~}{\psi }^2(k)).$$
(31)   
Perturbing $`\stackrel{~}{\psi }(k)`$ and $`n_k`$ and bearing in mind conditions (24) and (31), (6) gives Eq. (25) again. However, now $`\stackrel{~}{U}^{}(k)`$ obeys the new relation

$`\stackrel{~}{U}^{}(k)`$ $`=`$ $`{\displaystyle \left(\phi _{𝐤/2}^2(𝐫)\phi ^2(r)\right)\mathrm{\Phi }(r)d^3r}`$ (33)   
$`{\displaystyle \frac{d^3q}{(2\pi )^3}\frac{\stackrel{~}{U}(q)\left(\stackrel{~}{\psi }_{𝐤/2}^2(𝐪)\stackrel{~}{\psi }^2(q)\right)}{\stackrel{~}{\psi }(q)}}`$   
which significantly differs from (27). Indeed, the choice of the pair wave functions as real quantities implies that operating with integrands in (26) and (33), one can exploit $`\psi _𝐩(𝐫)\sqrt{2}\psi (r)p^2`$ at small $`p`$ . For $`k0`$ this provides $`\stackrel{~}{U}^{}(k)\stackrel{~}{U}(k)=t_k=ck^4+\mathrm{}.`$ Similar to Eq. (20), Eq. (25) can yields results correct only to the leading order in $`(nn_0)/n`$. So, it has to be solved together with (11) where $`n_0^2`$ should be replaced by $`n^2`$, rather than with (31). This leads to the following relation:

$`n_k={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\stackrel{~}{T}_k+n\stackrel{~}{U}(k)}{\sqrt{\stackrel{~}{T}_k^2+2n\stackrel{~}{T}_k\stackrel{~}{U}(k)}}}1\right),`$ (34)   
$`\stackrel{~}{\psi }(k)={\displaystyle \frac{\stackrel{~}{U}(k)}{2\sqrt{\stackrel{~}{T}_k^2+2n\stackrel{~}{T}_k\stackrel{~}{U}(k)}}},`$ (35)   
where $`\stackrel{~}{T}_k=T_k+nt_k`$. In the limit $`k0`$ Eq. (35) gives $`n_k(\sqrt{nm\stackrel{~}{U}(0)}/\mathrm{}k1)/2`$, which is fully consistent with the “$`1/k^2`$” theorem of Bogoliubov for the zero temperature . Eqs. (26) and (35) should be solved in a self–consistent manner. So, for $`n0`$ one can derive

$$\stackrel{~}{U}(k)=\stackrel{~}{U}^{(0)}(k)(1+8\sqrt{na^3}/\sqrt{\pi }).$$
(36)   
Here $`\stackrel{~}{U}^{(0)}(k)=\phi ^{(0)}(r)\mathrm{\Phi }(r)\mathrm{exp}(i\mathrm{𝐤𝐫})d^3r`$, where $`\phi ^{(0)}(r)`$ obeys Eq. (28). Further, substituting $`k=\sqrt{n}y`$ in the integral for the condensate depletion $`(nn_0)/n=1/(2\pi )^3_0^+\mathrm{}𝑑k\mathrm{\hspace{0.17em}4}\pi k^2n_k/n`$, we obtain the familiar result

$$(nn_0)/n=8\sqrt{na^3}/(3\sqrt{\pi })+\mathrm{},$$
(37)   
$`a`$ being the scattering length. Inserting (24), (34) and (35) into Eq. (6) and using (36), in a similar manner we derive

$$\epsilon =\frac{2\pi \mathrm{}^2an}{m}\left\{1+\frac{128}{15\sqrt{\pi }}\sqrt{na^3}\left(1\frac{5}{8}\frac{b}{a}\right)+\mathrm{}\right\},$$
(38)   
where $`b0`$ is one more characteristic length defined as

$$b=\frac{1}{4\pi }(\phi ^{(0)}(r))^2d^3r.$$
(39)   
As it is seen, the well–known result of papers  can be derived from (38) with the choice $`b=0`$. However, this approximation is rather crude because the case of the hard–sphere interaction ($`\mathrm{\Phi }(r)=0(r>a)`$ and $`\mathrm{\Phi }(r)\mathrm{}(r<a)`$) is specified by $`b=a`$:

$$\epsilon =\frac{2\pi \mathrm{}^2an}{m}\left\{1+\frac{16}{5\sqrt{\pi }}\sqrt{na^3}+\mathrm{}\right\}.$$
(40)   
In the general case and, in particular, for the singular potentials of the Lennard–Jones type we have $`ab`$. Remark that the last term in the r.h.s. of (30) does not make any contribution into the results given by (37) and (38). However, the next orders in the expansions of the energy and depletion depend on its contribution essentially.

Concluding let us take notice of the important points of this Letter once more. It was demonstrated that thermodynamically consistent calculations based on (10) and (11) conclusively result in Eqs. (20)-(23). Therefore, using the Bogoliubov model with the “dressed” interaction does not provide the satisfactory solution of the problem of the strong–coupling Bose gas. As it was shown, when investigating this subject, one should go beyond the Bogoliubov scheme. To do this, we developed the approach reduced to the system of Eqs. (26), (33), (34) and (35). This equations reproduce the well–known result (37) for the condensate depletion and yields the new expansion (38) in powers of $`n`$ for the energy, (40) being the particular case of the hard spheres. One can expect alterations for the excitation spectrum, too.

This work was supported by the RFBR Grant No. 97-02-16705.