text:
# Extension of Bogoliubov theory to quasi-condensates

## I Construction of a discrete phase-density representation

### I.1 Why discretize the real space ?

In previous studies of quasi-condensates the basic tools of the theory are an operator

$$\widehat{\rho }(𝐫)=\widehat{\psi }^{}(𝐫)\widehat{\psi }(𝐫)$$
(1)
giving the density in $`𝐫`$ and an operator $`\widehat{\theta }(𝐫)`$ giving the phase of $`\widehat{\psi }(𝐫)`$, the field operator in $`𝐫`$, the position $`𝐫`$ being a continuous variable griffin . A small parameter of the theory characterizing the regime of quasi-condensates is then that the density fluctuations, that is the fluctuations of $`\widehat{\rho }(𝐫)`$, are small in relative values:

$$\text{Var}(\widehat{\rho }(𝐫))\widehat{\rho }(𝐫)^2\widehat{\rho }(𝐫)^2\widehat{\rho }(𝐫)^2.$$
(2)
However one finds that the expectation value of $`\widehat{\rho }(𝐫)^2`$ is infinite in every point with a non-vanishing mean density $`\rho (𝐫)=\widehat{\rho }(𝐫)`$:

$$\widehat{\rho }(𝐫)^2=\delta (\mathrm{𝟎})\rho (𝐫)+\widehat{\psi }^{}(𝐫)\widehat{\psi }^{}(𝐫)\widehat{\psi }(𝐫)\widehat{\psi }(𝐫)$$
(3)
where $`\delta (\mathrm{𝟎})`$, the value of the Dirac distribution at the origin, is infinite, and the second term in the right hand side, giving the probability density of finding two atoms in the same point of space, is finite in any realistic model. Mathematically this divergence is due to the use of the bosonic commutation relations of the field operators $`\widehat{\psi }(𝐫)`$ and $`\widehat{\psi }^{}(𝐫)`$ in the same point of space to put the atomic field product in normal order.

In order to have small, and therefore finite, density fluctuations, one is forced to discretize the space, that is to collect the particles in little boxes at the nodes of a spatial grid. Each little box has equal length $`l`$ along each dimension of space and is parameterized by the position $`𝐫`$ of its center. The field operator $`\widehat{\psi }(𝐫)`$ has the effect of removing a particle in the box at position $`𝐫`$ and it now satisfies the bosonic commutation relations:

$$[\widehat{\psi }(𝐫),\widehat{\psi }^{}(𝐫^{})]=\frac{\delta _{𝐫,𝐫^{}}}{l^D}$$
(4)
where $`\delta _{𝐫,𝐫^{}}`$ is the discrete Kronecker delta and $`D`$ is the dimension of space. The variance of the operator giving the density is now finite:

$`\text{Var}(\widehat{\rho }(𝐫))`$ $`=`$ $`\widehat{\psi }^{}(𝐫)\widehat{\psi }^{}(𝐫)\widehat{\psi }(𝐫)\widehat{\psi }(𝐫)\rho ^2(𝐫)`$ (5)
$`+`$ $`{\displaystyle \frac{\rho (𝐫)}{l^D}}.`$
In the validity domain of the theoretical approach of this paper, this variance will be much smaller than $`\rho ^2(𝐫)`$ because both the term on the first line of the right hand side and the term on the second line are small:

$`|\widehat{\psi }^{}(𝐫)\widehat{\psi }^{}(𝐫)\widehat{\psi }(𝐫)\widehat{\psi }(𝐫)\rho ^2(𝐫)|\rho ^2(𝐫)`$ (6)
$`\rho (𝐫)l^D1.`$ (7)
### I.2 The phase operator

In usual continuous space theories an hermitian field phase operator $`\widehat{\theta }(𝐫)`$ is introduced subject to the following commutation relation with the operator giving the density:

$$[\widehat{\rho }(𝐫),\widehat{\theta }(𝐫^{})]=i\delta (𝐫𝐫^{}).$$
(8)
In our discrete model the desired commutation relation is modified into

$$[\widehat{\rho }(𝐫),\widehat{\theta }(𝐫^{})]=i\frac{\delta _{𝐫,𝐫^{}}}{l^D}.$$
(9)
First we recall briefly that there exists actually no hermitian operator $`\widehat{\theta }(𝐫)`$ satisfying strictly the above commutation relation. From the identity (9) one can indeed show that the operator

$$T(\alpha )e^{i\alpha \widehat{\theta }(𝐫)},$$
(10)
where $`\alpha `$ is any real number, is a translation operator for the density hint :

$$T(\alpha )^{}\widehat{\rho }(𝐫)T(\alpha )=\widehat{\rho }(𝐫)+\frac{\alpha }{l^D}.$$
(11)
This identity contradicts two fundamental properties of $`\widehat{\rho }(𝐫)`$, the positiveness and the discreteness of its spectrum continuous .

We now proceed with the construction of a phase operator $`\widehat{\theta }(𝐫)`$ approximately satisfying the commutation relation (9). The key ingredients allowing such an approximate construction are (i) to be in the limit of a large occupation number of the considered box of the lattice, and (ii) to construct the operator $`e^{i\widehat{\theta }}`$ first, which, according to (11) taken with $`\alpha =1`$, simply reduces the number of particles in the considered box by one.

In each spatial box we introduce the basis of Fock states $`|n,𝐫`$ with exactly $`n`$ particles in the box. In this basis the field operators have the following matrix elements:

$`\widehat{\psi }(𝐫^{})|n,𝐫`$ $`=`$ $`{\displaystyle \frac{\delta _{𝐫,𝐫^{}}}{l^{D/2}}}\sqrt{n}|n1,𝐫,`$
$`\widehat{\psi }^{}(𝐫^{})|n,𝐫`$ $`=`$ $`{\displaystyle \frac{\delta _{𝐫,𝐫^{}}}{l^{D/2}}}\sqrt{n+1}|n+1,𝐫`$ (12)
as a consequence of the commutation relation (4). The atomic density $`\widehat{\rho }`$ defined by $`\widehat{\rho }(𝐫)=\widehat{\psi }^{}(𝐫)\widehat{\psi }(𝐫)`$ is diagonal in the Fock state basis:

$$\widehat{\rho }(𝐫^{})|n,𝐫=\delta _{𝐫,𝐫^{}}\frac{n}{l^D}|n,𝐫.$$
(13)
We then introduce the operator $`\widehat{A}`$ defined by:

$$\widehat{\psi }(𝐫)\widehat{A}(𝐫)\sqrt{\widehat{\rho }(𝐫)}.$$
(14)
In the Fock space $`\widehat{A}(𝐫)`$ reduces by one the number of particles $`n`$ in the box $`𝐫`$:

$$\widehat{A}(𝐫^{})|n,𝐫=\left(1\delta _{n,0}\right)\delta _{𝐫,𝐫^{}}|n1,𝐫.$$
(15)
Note that its action on the vacuum state of the box gives zero. For each box $`𝐫`$, the definition of $`\widehat{A}`$ leads to the exact relations:

$`\widehat{A}\widehat{A}^{}=I,\widehat{A}^{}\widehat{A}=I|00|\text{and }[\widehat{A},\widehat{A}^{}]=|00|,`$ (16)
where $`I`$ is the identity operator and $`|0`$ is the zero-particle state or vacuum state in the box of center $`𝐫`$. We find that the operator $`\widehat{A}`$ is almost unitary, i.e. it is effectively unitary for a physical state of the system with a negligible probability of having an empty box. In what follows, we assume that this condition is satisfied, so that the projector $`|00|`$ can be neglected:

$$\text{occupation probability of}|n=0,𝐫1.$$
(17)
In this case, we write the approximately unitary operator $`\widehat{A}`$ as:

$$\widehat{A}(𝐫)e^{i\widehat{\theta }(𝐫)}\text{with}\widehat{\theta }^{}(𝐫)\widehat{\theta }(𝐫)$$
(18)
which amounts to writing the field operator as:

$$\widehat{\psi }(𝐫)e^{i\widehat{\theta }(𝐫)}\sqrt{\widehat{\rho }(𝐫)}.$$
(19)
This should be understood as a formal writing, allowing for example to recover the matrix elements of $`\widehat{A}`$ and therefore of the field operator $`\widehat{\psi }`$ from the commutation relation (9). We summarize below all the commutation relations of our phase-density representation:

$`[\widehat{\rho }(𝐫),\widehat{\theta }(𝐫^{})]{\displaystyle \frac{i\delta _{𝐫,𝐫^{}}}{l^D}}[\widehat{\rho }(𝐫),\widehat{\rho }(𝐫^{})]=0[\widehat{\theta }(𝐫),\widehat{\theta }(𝐫^{})]0`$ (20)
We come back to the constraint (17) at the basis of the construction of $`\mathrm{exp}(i\widehat{\theta })`$. A sufficient condition to have a low probability for a zero particle occupation in a box is obtained for a large mean number of particles in the box and with small relative particle number fluctuations. This is the regime that we wish to consider in this paper. We are therefore back to the discussion of the previous subsection and to the conditions (6,7) for weak density fluctuations. In particular, the construction of the operator $`\mathrm{exp}(i\widehat{\theta })`$ becomes problematic in the limit $`l0`$, that is in the continuous model.

### I.3 How to choose the grid spacing $`l`$

Working on a grid can be also seen as performing a coarse-grain average over all physical quantities on a scale $`l`$. This averaging suppresses the short wavelength modes (shorter than $`l`$) and thus introduces an energy cutoff:

$$E_{\mathrm{cut}}\frac{\mathrm{}^2}{ml^2}.$$
(21)
This cutoff is of no physical consequence if all characteristic energies ($`\mu `$, $`k_BT`$) are smaller, i.e., that $`l`$ is smaller than the corresponding characteristic lengths. This leads, for instance, to the following restrictions for $`l`$:

$`l<\xi \text{and}l<\lambda `$ (22)
where

$$\xi =\frac{\mathrm{}}{\sqrt{m\mu }}$$
(23)
is the healing length, and

$$\lambda =\sqrt{\frac{2\pi \mathrm{}^2}{mk_BT}}$$
(24)
is the thermal de Broglie wavelength. These two restrictions, combined with Eq. (7), impose:

$`\rho \lambda ^D1,`$ (25)
$`\rho \xi ^D1.`$ (26)
These are conditions of validity for our discrete model.

The first one, Eq. (25), is the quantum degeneracy regime occurring at sufficiently low temperatures. The second restriction, Eq. (26), corresponds to the regime of weakly interacting systems. Its dependence on the density varies according to the dimension of space. In 1D and 3D, the mean field prediction for the chemical potential is $`\mu g\rho `$ where $`g`$ is a constant characterizing the interaction potential, the so-called coupling constant. In 1D, Eq.(26) is the high density limit where a mean field theory is valid; we recall that the small density limit $`\rho \xi =\mathrm{}\sqrt{\rho /mg}1`$ corresponds to the strongly interacting (or strongly correlated) Tonks gas regime. In 3D, the effective coupling constant $`g`$ is related to the $`s`$-wave scattering length $`a`$ of the interaction potential, $`g=4\pi \mathrm{}^2a/m`$, so that $`\rho \xi ^31/\sqrt{\rho a^3}1`$: one recovers the usual small gaseous parameter $`\sqrt{\rho a^3}`$. In 2D, the chemical potential scales as $`\mathrm{}^2\rho /(m\mathrm{ln}(1/\rho a^2))`$ where $`a`$ is the scattering length of the 2D interaction potential so that the condition $`\rho \xi ^21`$ results in a low density condition, $`\mathrm{ln}(1/\rho a^2)1`$.

## II Perturbative treatment of a model Hamiltonian

### II.1 Model Hamiltonian

In our lattice model, we represent the binary interaction potential among the particles by a discrete delta potential:

$$V(𝐫_1𝐫_2)=\frac{g_0}{l^D}\delta _{𝐫_1,𝐫_2}$$
(27)
where $`g_0`$ is the bare coupling constant. Note that $`g_0`$ in general differs from the effective coupling constant $`g`$, and we shall come back to this point in section III.1. With this model potential, the grand canonical Hamiltonian is:

$$H=\underset{𝐫}{}l^D\left[\frac{\mathrm{}^2}{2m}\widehat{\psi }^{}(𝐫)\mathrm{\Delta }\widehat{\psi }(𝐫)+(U(𝐫)\mu )\widehat{\psi }^{}(𝐫)\widehat{\psi }(𝐫)+\frac{g_0}{2}\widehat{\psi }^{}(𝐫)\widehat{\psi }^{}(𝐫)\widehat{\psi }(𝐫)\widehat{\psi }(𝐫)\right]$$
(28)
where $`U(𝐫)`$ is an external trapping potential and where the Laplacian is a symmetric operator coupling the different neighboring boxes:

$$\mathrm{\Delta }f(𝐫)=\underset{j}{}\frac{f(𝐫+l𝐞_𝐣)+f(𝐫l𝐞_𝐣)2f(𝐫)}{l^2}.$$
(29)
The $`𝐞_𝐣`$ are the unitary vectors and $`j`$ the different orthogonal space directions (for example, $`j=x,y,z`$ in 3D). As usual we take periodic boundary conditions inside a rectangular box with lengths integer multiple of $`l`$.

We now rewrite the Hamiltonian in the density-phase representation, that is in term of the operators $`\widehat{\rho }`$ and $`\widehat{\theta }`$ giving the density and the phase as defined in the previous section. The contributions of the trapping potential and of the interaction potential to the Hamiltonian are local in real space and therefore include the operator $`\widehat{\rho }`$ only:

$$H_{\mathrm{pot}}=\underset{𝐫}{}l^D\widehat{\rho }(𝐫)\left[U(𝐫)\mu +\frac{g_0}{2}\left(\widehat{\rho }(𝐫)\frac{1}{l^D}\right)\right]$$
(30)
where we have used the bosonic commutation relation (4) to exchange one of the $`\widehat{\psi }^{}`$ with $`\widehat{\psi }`$ in the interaction term. The kinetic energy term involves explicitly the phase operator:

$$H_{\mathrm{kin}}=\frac{\mathrm{}^2}{2ml^2}\underset{𝐫}{}l^D\underset{j}{}\sqrt{\widehat{\rho }}\left(e^{i(\widehat{\theta }_{+j}\widehat{\theta })}\sqrt{\widehat{\rho }_{+j}}+e^{i(\widehat{\theta }_j\widehat{\theta })}\sqrt{\widehat{\rho }_j}2\sqrt{\widehat{\rho }}\right)$$
(31)
where we have introduced the notations $`\widehat{\theta }_{\pm j}=\widehat{\theta }(𝐫\pm l𝐞_𝐣)`$ and $`\widehat{\rho }_{\pm j}=\widehat{\rho }(𝐫\pm l𝐞_𝐣)`$. A remarkable property of this writing, to be used below, is that it involves only the difference of two phase operators between two neighboring points of the lattice.

### II.2 Hamiltonian quadratisation and cubisation

We now expand the Hamiltonian to third order in powers of two small parameters. As already discussed in subsection I.1 the regime of quasi-condensates that we are interested in corresponds to small relative fluctuations $`\delta \widehat{\rho }`$ of the density. In the zeroth order approach totally neglecting the density fluctuations, the density is set to a deterministic value $`\rho _0`$, as we shall see. The second order expansion allows us to describe the density fluctuations:

$$\widehat{\rho }(𝐫)=\rho _0(𝐫)+\delta \widehat{\rho }(𝐫).$$
(32)
The third order expansion allows us to calculate the mean value of $`\delta \widehat{\rho }(𝐫)`$.

The first small parameter of the systematic expansion used in this paper is therefore given by

$$ϵ_1=\frac{|\delta \widehat{\rho }|}{\rho _0}1.$$
(33)
where $`|\delta \widehat{\rho }|`$ is the typical value of the operator $`\delta \widehat{\rho }`$ in the physical state of the system. Mathematically this allows to expand $`\sqrt{\widehat{\rho }}`$ as

$$\widehat{\rho }^{1/2}=\rho _0^{1/2}+\frac{1}{2}\frac{\delta \widehat{\rho }}{\rho _0^{1/2}}\frac{1}{8}\frac{\delta \widehat{\rho }^2}{\rho _0^{3/2}}+\frac{1}{16}\frac{\delta \widehat{\rho }^3}{\rho _0^{5/2}}\mathrm{}$$
(34)
The second small parameter of the expansion is given by

$$ϵ_2=|l\widehat{\theta }|1.$$
(35)
Here $``$ represents the gradient on the lattice:

$$f(𝐫)=\underset{j}{}\frac{f(𝐫+l𝐞_𝐣)f(𝐫l𝐞_𝐣)}{2l}𝐞_𝐣$$
(36)
where $`f`$ is an arbitrary function. Physically the existence of the small parameter $`ϵ_2`$ is reasonable: it is at the basis of our hypothesis that the continuous quantum field problem can be well approximated by a discrete lattice model, provided that $`l`$ is small enough, see subsection I.3. Mathematically this second small parameter allows to expand in (31) the exponentials of the phase differences:

$$e^{i(\widehat{\theta }_{+j}\widehat{\theta })}=1+i(\widehat{\theta }_{+j}\widehat{\theta })\frac{1}{2}(\widehat{\theta }_{+j}\widehat{\theta })^2\mathrm{}$$
(37)
From the fact that the discretisation length $`l`$ is on the order of the smallest of the two macroscopic length scales $`\xi `$ and $`\lambda `$, see (22), it will be checked later that the parameters $`ϵ_1`$ and $`ϵ_2`$, though of an apparently different physical origin, can be chosen to be of the same order of magnitude,

$$ϵ_1ϵ_2\frac{1}{\sqrt{\rho _0l^D}}$$
(38)
and can therefore be treated mathematically as infinitesimals of the same order. The mathematical details of the expansion

$$H=H_0+H_1+H_2+H_3+\mathrm{}$$
(39)
are given in the appendix A, we present here only the results:

$`H_0`$ $`=`$ $`{\displaystyle \underset{𝐫}{}}l^D\left[{\displaystyle \frac{\mathrm{}^2}{2m}}\sqrt{\rho _0}\mathrm{\Delta }\sqrt{\rho _0}+{\displaystyle \frac{g_0}{2}}\rho _0^2+(U(𝐫)\mu )\rho _0\right]`$
$`H_1`$ $`=`$ $`{\displaystyle \underset{𝐫}{}}l^D\left[{\displaystyle \frac{\mathrm{}^2}{2m\sqrt{\rho _0}}}\mathrm{\Delta }\sqrt{\rho _0}+U(𝐫)\mu +g_0\rho _0\right]\delta \widehat{\rho }`$
$`H_2`$ $`=`$ $`E_2[\rho _0]+{\displaystyle \underset{𝐫}{}}l^D[{\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{\delta \widehat{\rho }}{2\sqrt{\rho _0}}}\mathrm{\Delta }\left({\displaystyle \frac{\delta \widehat{\rho }}{2\sqrt{\rho _0}}}\right)+{\displaystyle \frac{\mathrm{}^2\delta \widehat{\rho }^2}{8m\rho _0^{3/2}}}\mathrm{\Delta }\sqrt{\rho _0}+{\displaystyle \frac{g_0}{2}}\delta \widehat{\rho }^2`$ (40)
$`+`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \underset{j}{}}\sqrt{\rho _0(𝐫)\rho _0(𝐫+l𝐞_𝐣)}{\displaystyle \frac{(\widehat{\theta }(𝐫+l𝐞_𝐣)\widehat{\theta }(𝐫))^2}{l^2}}]`$
$`H_3`$ $`=`$ $`{\displaystyle \frac{g_0}{2}}{\displaystyle \underset{𝐫}{}}\delta \widehat{\rho }+{\displaystyle \frac{\mathrm{}^2}{4ml^2}}{\displaystyle \underset{𝐫,j}{}}l^D(\widehat{\theta }_{+j}\widehat{\theta })\left({\displaystyle \frac{\rho _{0,+j}^{1/2}}{\rho _0^{1/2}}}\delta \widehat{\rho }+{\displaystyle \frac{\rho _0^{1/2}}{\rho _{0,+j}^{1/2}}}\delta \widehat{\rho }_{+j}\right)(\widehat{\theta }_{+j}\widehat{\theta })`$
$`+`$ $`{\displaystyle \frac{\mathrm{}^2}{8m}}{\displaystyle \underset{𝐫}{}}{\displaystyle \frac{\delta \widehat{\rho }}{\rho _0}}\left(\rho _0^{1/2}\mathrm{\Delta }\rho _0^{1/2}\rho _0^{1/2}\mathrm{\Delta }\rho _0^{1/2}\right){\displaystyle \frac{\mathrm{}^2}{16m}}{\displaystyle \underset{𝐫}{}}l^D\left[{\displaystyle \frac{\delta \widehat{\rho }^3}{\rho _0^{5/2}}}\mathrm{\Delta }\sqrt{\rho _0}{\displaystyle \frac{\delta \widehat{\rho }^2}{\rho _0^{3/2}}}\mathrm{\Delta }\left({\displaystyle \frac{\delta \widehat{\rho }}{\sqrt{\rho _0}}}\right)\right].`$
The quantity $`E_2`$ in $`H_2`$ is a c-number functional of the density $`\rho _0`$, given in the appendix A, which has therefore no contribution to the dynamics of the quantum field.

### II.3 Iterative solution for the quadratic Hamiltonian

We now solve perturbatively, order by order, the Hamiltonian problems defined by $`H_0`$, $`H_0+H_1`$ and $`H_0+H_1+H_2`$. To zeroth order in $`ϵ_{1,2}`$, the Hamiltonian is a c-number. As the chemical potential is fixed in our approach, $`H_0`$ is minimized for a density profile $`\rho _0(𝐫)`$ such that $`\sqrt{\rho _0}`$ solves the discrete version of the Gross-Pitaevskii equation:

$$\left[\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U(𝐫)\mu +g_0\rho _0\right]\sqrt{\rho _0}=0.$$
(41)
This density profile constitutes the zeroth order approximation to the density $`\rho `$. It contains a number of particles that we call $`N_0`$:

$$N_0\underset{𝐫}{}l^D\rho _0(𝐫).$$
(42)
Note that $`N_0`$ coincides with the mean total number of particles $`N`$ only to lowest order in the theory. Equation (41) defines $`\rho _0`$ and therefore $`N_0`$ as function of the chemical potential $`\mu `$. It will however turn out to be more convenient to parameterize the theory in terms of $`N_0`$ rather than in terms of $`\mu `$. We will therefore consider $`\mu `$ and $`\rho _0`$ as functions of $`N_0`$:

$`\mu `$ $`=`$ $`\mu _0(N_0)`$ (43)
$`\rho _0(𝐫)`$ $`=`$ $`\rho _0(𝐫;N_0).`$ (44)
$`\mu _0`$ is therefore the Gross-Pitaevskii prediction for the chemical potential of a gas of $`N_0`$ particles.

For the choice of density profile (41), the first order correction $`H_1`$ to the Hamiltonian vanishes. We therefore have now to solve the Hamiltonian problem defined by $`H_2`$, in order to determine the lowest order approximation to the density fluctuation $`\delta \widehat{\rho }`$ and the phase $`\widehat{\theta }`$. It is instructive to write the corresponding Heisenberg equations of motion, which are linear (and therefore trivially solvable) since $`H_2`$ is quadratic. As $`\widehat{\theta }`$ and $`\delta \widehat{\rho }`$ are two canonically conjugate variables, the equations of motion are:

$$\begin{array}{ccccc}\hfill \mathrm{}_t\widehat{\theta }& & \frac{1}{l^D}\frac{H_2}{\delta \widehat{\rho }(𝐫)}\hfill & =& \frac{\mathrm{}^2}{2m\sqrt{\rho _0}}\left[\mathrm{\Delta }\left(\frac{\delta \widehat{\rho }}{2\sqrt{\rho _0}}\right)\delta \widehat{\rho }\frac{\mathrm{\Delta }\sqrt{\rho _0}}{2\rho _0}\right]g_0\delta \widehat{\rho }\hfill \\ \hfill \mathrm{}_t\delta \widehat{\rho }(𝐫)& & \frac{1}{l^D}\frac{H_2}{\widehat{\theta }(𝐫)}\hfill & =& \frac{\mathrm{}^2}{m}\sqrt{\rho _0}\left[\mathrm{\Delta }\left(\widehat{\theta }\sqrt{\rho _0}\right)\widehat{\theta }\mathrm{\Delta }\sqrt{\rho _0}\right]\hfill \end{array}$$
(45)
An important difference of these equations with the so-called quantum hydrodynamics equations for $`\delta \widehat{\rho }`$ and $`\widehat{\theta }`$ is that our formalism keeps the so-called quantum pressure term for $`_t\widehat{\theta }`$ whereas it is usually neglected in the literature griffin . This allows our treatment to have a cut-off energy larger than $`\mu `$, whereas the usual treatment is restricted to energy modes much below $`\mu `$.

Furthermore one can simplify these equations using the Gross-Pitaevskii equation (41) to eliminate $`\mathrm{\Delta }\sqrt{\rho _0}`$:

$`\mathrm{}_t\widehat{\theta }`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{\rho _0}}}\left[{\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{\Delta }+U+3g_0\rho _0\mu \right]\left({\displaystyle \frac{\delta \widehat{\rho }}{\sqrt{\rho _0}}}\right)`$ (46)
$`\mathrm{}_t\delta \widehat{\rho }(𝐫)`$ $`=`$ $`2\sqrt{\rho _0}\left[{\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{\Delta }+U+g_0\rho _0\mu \right]\left(\widehat{\theta }\sqrt{\rho _0}\right).`$ (47)
This gives the idea of a very simple canonical transformation which, remarkably, maps our equations for a quasi-condensate (46,47) into the equations for the Bogoliubov modes of a condensate: the field

$$\widehat{B}=\frac{\delta \widehat{\rho }}{2\sqrt{\rho _0}}+i\sqrt{\rho _0}\widehat{\theta }$$
(48)
has bosonic commutation relations

$$[\widehat{B}(𝐫),\widehat{B}^{}(𝐫^{})]=\frac{\delta _{𝐫,𝐫^{}}}{l^D}$$
(49)
and it obeys the standard Bogoliubov equations

$$i\mathrm{}_t\left(\begin{array}{c}\widehat{B}\\ \widehat{B}^{}\end{array}\right)=_{\mathrm{GP}}\left(\begin{array}{c}\widehat{B}\\ \widehat{B}^{}\end{array}\right)\left(\begin{array}{cc}\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U\mu +2g_0\rho _0& g_0\rho _0\\ g_0\rho _0& \left(\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U\mu +2g_0\rho _0\right)\end{array}\right)\left(\begin{array}{c}\widehat{B}\\ \widehat{B}^{}\end{array}\right).$$
(50)
This mapping can be readily extended to the Hamiltonian $`H_2`$, which is expected to be canonically equivalent to the Bogoliubov Hamiltonian:

$$H_2=l^D\underset{𝐫}{}\widehat{B}^{}\left(\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U+g_0\rho _0\mu \right)\widehat{B}+g_0\rho _0\left[\widehat{B}^{}\widehat{B}+\frac{1}{2}\left(\widehat{B}^2+\widehat{B}^2\right)\right].$$
(51)
We have checked that the identity (51) indeed holds by replacing $`\widehat{B}`$ by its expression (48) in terms of $`\delta \widehat{\rho }`$ and $`\widehat{\theta }`$, and by using the value of the commutators (20) and the fact that $`\sqrt{\rho _0}`$ solves the Gross-Pitaevskii equation. Remarkably the energy functional $`E_2[\rho _0]`$ is exactly compensated by the contribution of the commutators.

This mapping therefore allows to reuse the standard diagonalisation of the Bogoliubov Hamiltonian. We recall here briefly the procedure described in lewenstein ; yvan2 . One introduces the normal eigenmodes $`(u_s,v_s)`$ of the Bogoliubov operator $`_{\mathrm{GP}}`$ with an energy $`ϵ_s`$, normalisable as

$$\underset{𝐫}{}l^D[|u_s(𝐫)|^2|v_s(𝐫)|^2]=1.$$
(52)
Then $`(v_s^{},u_s^{})`$ is an eigenmode of $`_{\mathrm{GP}}`$ with the energy $`ϵ_s`$. To form a complete family of vectors one has to further introduce the zero energy mode of $`_{\mathrm{GP}}`$, given by $`(\varphi _0,\varphi _0)`$, and the anomalous mode $`(\varphi _a,\varphi _a)`$ with

$`\varphi _0=\sqrt{\rho _0/N_0}\text{and}\varphi _a=\sqrt{N_0}_{N_0}\sqrt{\rho _0}.`$ (53)
The corresponding normalization of the anomalous mode is such that the scalar product of $`\varphi _0`$ and $`\varphi _a`$ is 1/2. With these definitions, one introduces the components of $`(\widehat{B},\widehat{B}^{})`$ on the zero energy mode, on the anomalous mode and on the regular $`(u_s,v_s)`$ modes:

$$\left(\begin{array}{c}\widehat{B}\\ \widehat{B}^{}\end{array}\right)=i\sqrt{N_0}\widehat{Q}\left(\begin{array}{c}\varphi _0\\ \varphi _0\end{array}\right)+\frac{\widehat{P}}{\sqrt{N_0}}\left(\begin{array}{c}\varphi _a\\ \varphi _a\end{array}\right)+\underset{s}{}\widehat{b}_s\left(\begin{array}{c}u_s\\ v_s\end{array}\right)+\widehat{b}_s^{}\left(\begin{array}{c}v_s^{}\\ u_s^{}\end{array}\right).$$
(54)
$`\widehat{Q}`$ is a collective coordinate representing the quantum phase of the field and $`\widehat{P}`$ is its conjugate momentum

$$[\widehat{P},\widehat{Q}]=i.$$
(55)
Physically $`\widehat{P}`$ corresponds to fluctuations in the total number of particles, as expected, and as shown in more details later, see (67). The operators $`\widehat{b}_s`$ are bosonic annihilation operators with the usual commutation relations $`[\widehat{b}_s,\widehat{b}_s^{}]=\delta _{s,s^{}}`$. They commute with $`\widehat{P}`$ and $`\widehat{Q}`$. The inverse formulas giving $`\widehat{b}_s,\widehat{Q}`$ and $`\widehat{P}`$ in terms of $`\widehat{B}`$ can be found for example in yvan2 . Equation (54) results in the following modal expansion for the density fluctuations and the phase operators:

$$\begin{array}{ccc}\hfill \widehat{\theta }(𝐫)& =\hfill & \underset{s}{}\theta _s(𝐫)\widehat{b}_s+\theta _s^{}(𝐫)\widehat{b}_s^{}\widehat{Q}\\ \hfill \delta \widehat{\rho }(𝐫)& =\hfill & \underset{s}{}\delta \rho _s(𝐫)\widehat{b}_s+\delta \rho _s^{}(𝐫)\widehat{b}_s^{}+\widehat{P}_{N_0}\rho _0\end{array}$$
(56)
where

$`\theta _s(𝐫)`$ $`=`$ $`{\displaystyle \frac{u_s(𝐫)v_s(𝐫)}{2i\sqrt{\rho _0(𝐫)}}},`$
$`\delta \rho _s(𝐫)`$ $`=`$ $`\sqrt{\rho _0(𝐫)}(u_s(𝐫)+v_s(𝐫)).`$ (57)
By construction, this modal expansion, when inserted into the quadratic Hamiltonian $`H_2`$, results in

$$H=\underset{s}{}ϵ_s\widehat{b}_s^{}\widehat{b}_s+\frac{1}{2}\widehat{P}^2\mu _0^{}+\stackrel{~}{E}_2[\rho _0],$$
(58)
where $`\mu _0^{}=d\mu _0/dN_0`$. This is the sum of uncoupled harmonic oscillators, plus a massive free degree of freedom corresponding to the unbound phase variable $`\widehat{Q}`$. The effective mass of the phase variable is given by $`1/\mu _0^{}`$. The energy functional $`\stackrel{~}{E}_2[\rho _0]`$ will be calculated in III.2, where it will be shown that it leads to exactly the same ground state energy as the number conserving Bogoliubov theory. This shows that the Bogoliubov theory can be used to calculate the ground state energy even for e.g. 1D quasi-condensates, a fact commonly used in the literature lieb ; gaudin  but which looks rather heuristic in the absence of justification.

### II.4 Effect of cubic Hamiltonian corrections on the density

The physics contained in the cubic term $`H_3`$ of the Hamiltonian is very rich. It includes interaction effects between the Bogoliubov modes of the previous section, allowing a generalization to quasi-condensates of the theory of energy shifts and Beliaev-Landau damping usually put forward for Bose-Einstein condensates Vincent ; Stringari\_Pitaevskii ; Shlyapnikov .

We are more modest here. Our motivation to include the cubic corrections is that the quadratic Hamiltonian $`H_2`$ brings actually no correction to the zeroth order approximation $`\rho _0`$ to the mean density, since the mean value of $`\delta \widehat{\rho }`$ vanishes at the level of the second order theory. This is highly non satisfactory as it brings some inconsistency in the calculation of an observable like $`g_1`$, the first order correlation function of the field: to get a non-trivial prediction for $`g_1`$ one has to include terms quadratic in the phase operator, which are second order in $`ϵ_2`$, which forces to also include second order corrections to the mean density, as will become very explicit in section III.

We therefore calculate the first order correction to the equations of motion of $`\delta \widehat{\rho }`$ and $`\widehat{\theta }`$ due to the cubic Hamiltonian term $`H_3`$, and we take the average over the quantum state corresponding to the density operator at thermal equilibrium for the Hamiltonian $`H_2`$. This gives source terms to add to the equations for the mean density and phase derived from $`H_2`$. We leave the details of the calculations for the appendix B and we give directly the result:

$`\mathrm{}_t\delta \widehat{\rho }_3`$ $`=`$ $`\rho _0^{1/2}\left[{\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{\Delta }+U+g_0\rho _0\mu \right]\left(2\rho _0^{1/2}\widehat{\theta }_3\right)`$ (59)
$`2\mathrm{}\sqrt{\rho _0}_t\widehat{\theta }_3`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{\Delta }+U+3g_0\rho _0\mu \right]\left({\displaystyle \frac{\delta \widehat{\rho }_3\widehat{B}^{}\widehat{B}_2}{\rho _0^{1/2}}}\right)`$ (60)
$`+`$ $`g_0\rho _0^{1/2}4\widehat{B}^{}\widehat{B}+\widehat{B}^2+B^2_22\widehat{P}^2_2\mu _0^{}_{N_0}\sqrt{\rho _0}`$
where the thermal average $`\mathrm{}_2`$ is taken with the unperturbed Hamiltonian $`H_2`$ and $`\mathrm{}_3`$ is taken with the perturbed Hamiltonian $`H_2+H_3`$ to first order in $`H_3`$. The expectation value of the ‘kinetic energy’ of the unbound phase variable in (58) is equal to $`k_BT/2`$ according to the equipartition theorem so that

$$\widehat{P}^2_2=\frac{k_BT}{\mu _0^{}}.$$
(61)
At equilibrium the expectation value of $`_t\delta \widehat{\rho }`$ and $`_t\widehat{\theta }`$ vanish. This fact is obvious for $`_t\delta \widehat{\rho }`$; it is less obvious for $`_t\widehat{\theta }`$ because of the presence of the unbound variable $`\widehat{Q}`$, we therefore produce a proof of that in the appendix C. We therefore have to solve (59,60) with the left-hand side set to zero. The first equation (59) imposes that the mean value of $`\widehat{\theta }`$ is position independent, a trivial result. In the second equation, the operator acting on $`\delta \widehat{\rho }_3`$ is strictly positive so that it is invertible and (60) determines the correction to the mean density in a unique way.

We now go through a sequence of transformations allowing us to get a physical understanding of the value of $`\delta \widehat{\rho }_3`$. The first step is to pull out the contribution of the ‘anomalous’ terms $`\widehat{P}`$, $`\widehat{Q}`$ in the modal expansion (54):

$$\widehat{B}(𝐫)i\sqrt{N_0}\widehat{Q}\varphi _0(𝐫)+\frac{1}{\sqrt{N_0}}\widehat{P}\varphi _a(𝐫)+\widehat{B}_n.$$
(62)
We calculate the expectation values of (60) involving the operator $`\widehat{B}`$, using the fact that all the crossed terms between the anomalous part and the operators $`b_s`$ have a vanishing expectation value:

$`\widehat{B}^{}\widehat{B}_2`$ $`=`$ $`{\displaystyle \frac{\varphi _a^2}{N_0}}\widehat{P}^2_2+N_0\varphi _0^2\widehat{Q}^2_2+\widehat{B}_n^{}\widehat{B}_n_2\varphi _a\varphi _0`$ (63)
$`\widehat{B}^2_2+\widehat{B}^2_2`$ $`=`$ $`2{\displaystyle \frac{\varphi _a^2}{N_0}}\widehat{P}^2_22N_0\varphi _0^2\widehat{Q}^2_2+\widehat{B}_n^2_2+\widehat{B}_n^2_2.`$ (64)
The term $`\varphi _0\varphi _a`$ in (63) comes from the non-commutation of $`\widehat{P}`$ and $`\widehat{Q}`$, see (55). The contributions of $`\widehat{Q}^2_2`$ in (63) and (64), when inserted into (60), are shown to compensate exactly when one uses the fact that $`\varphi _0`$ solves the Gross-Pitaevskii equation. This was expected from the $`U(1)`$ symmetry of the Hamiltonian: only differences of the phase operator in two points appear in the Hamiltonian, so that $`H`$ does not depend on $`\widehat{Q}`$ and the mean density does not depend on $`\widehat{Q}^2_2`$.

We therefore get an equation for $`\delta \widehat{\rho }_3`$ involving the expectation value of $`\widehat{P}^2`$ as a source term, and which looks rather involved:

$`0`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{\Delta }+U+3g_0\rho _0\mu \right]\left({\displaystyle \frac{\delta \widehat{\rho }_3\varphi _a^2N_0^1\widehat{P}^2_2(\widehat{B}_n^{}\widehat{B}_n_2\varphi _0\varphi _a)}{\rho _0^{1/2}}}\right)`$ (65)
$`+`$ $`g_0\rho _0^{1/2}4(\widehat{B}_n^{}\widehat{B}_n\varphi _0\varphi _a)+\widehat{B}_n^2+B_n^2_2+\widehat{P}^2_2\left[6g_0\rho _0^{1/2}_{N_0}\sqrt{\rho _0}2\mu _0^{}\right]_{N_0}\sqrt{\rho _0}`$
Fortunately the underlying physics is very simple and allows to predict the effect of this source term on the mean density. One first identifies the physical meaning of $`\widehat{P}`$ in (56). Using the well-known fact that the eigenmodes of $`_{\mathrm{GP}}`$ are orthogonal for the modified scalar product of signature $`(1,1)`$, one has yvan2 :

$$\varphi _0|u_s+\varphi _0|v_s\underset{𝐫}{}l^D\varphi _0(𝐫)\left[u_s(𝐫)+v_s(𝐫)\right]=0$$
(66)
so that the sum of $`\delta \rho _s`$ over all spatial nodes vanishes. As a consequence the operator $`\widehat{N}`$ giving the total number of particles in the gas is simply

$$\widehat{N}\underset{𝐫}{}l^D\widehat{\rho }(𝐫)=N_0+\widehat{P}$$
(67)
where we have used the identity

$$\underset{𝐫}{}l^D_{N_0}\rho _0(𝐫)=\frac{d}{dN_0}\underset{𝐫}{}l^D\rho _0(𝐫)=\frac{dN_0}{dN_0}=1.$$
(68)
The source terms involving $`\widehat{P}`$ therefore correspond to fluctuations in the total number of particles in the gas, due to the fact that we consider the grand canonical ensemble. The effect of these grand canonical fluctuations can be considered for the case of a pure quasi-condensate at the order of the present calculation so it is easy to calculate it directly. In the grand canonical ensemble the probability that the quasi-condensate has $`n`$ particles is

$$\mathrm{\Pi }_n\mathrm{exp}\left[\beta \left(E_0(n)\mu n\right)\right]$$
(69)
where $`E_0(n)`$ is the Gross-Pitaevskii energy for the density profile $`\rho _0(𝐫;n)`$:

$$E_0(n)=\underset{𝐫}{}l^D\left[\frac{\mathrm{}^2}{2m}\sqrt{\rho _0(𝐫;n)}\mathrm{\Delta }\sqrt{\rho _0(𝐫;n)}+U(𝐫)\rho _0(𝐫;n)+\frac{g_0}{2}\rho _0^2(𝐫;n)\right].$$
(70)
The corresponding mean grand canonical density is

$$\rho _{\mathrm{GC}}(𝐫)=𝑑n\mathrm{\Pi }_n\rho _0(𝐫;n)$$
(71)
where we treat $`n`$ as a continuous variable. The zeroth order approximation $`n=N_0`$ for the number of particles in the quasi-condensate is such that $`E_0(n)\mu n`$ has a minimum:

$$\frac{d}{dn}[E_0(n)\mu n]=\mu _0(n)\mu =0\text{for}n=N_0$$
(72)
as shown in (43). The corresponding density profile is $`\rho _0(𝐫;N_0)`$. The next order correction to that is obtained by expanding the $`n`$-dependent density profile to second order in $`nN_0`$ and by averaging over $`n`$:

$$\delta \rho _{\mathrm{GC}}(𝐫)=(nN_0)_{N_0}\rho _0(𝐫;N_0)+\frac{1}{2}(nN_0)^2_{N_0}^2\rho _0(𝐫;N_0).$$
(73)
The second moment of $`nN_0`$ is calculated to lowest non-vanishing order by a Gaussian approximation to $`\mathrm{\Pi }_n`$:

$$E_0(n)\mu nE_0(N_0)\mu N_0+\frac{1}{2}\frac{d^2E_0}{dN_0^2}(nN_0)^2=\text{const}+\frac{1}{2}\mu _0^{}(nN_0)^2.$$
(74)
This leads to

$$(nN_0)^2_{\mathrm{Gauss}}=\frac{k_BT}{\mu _0^{}}=\widehat{P}^2_2.$$
(75)
More care has to be taken in the calculation of the mean of $`nN_0`$: the Gaussian approximation to $`\mathrm{\Pi }_n`$ gives a vanishing contribution, so that the cubic distortion to it has to be included:

$`E_0(n)\mu n`$ $``$ $`E_0(N_0)\mu N_0+{\displaystyle \frac{1}{2}}\mu _0^{}(nN_0)^2+{\displaystyle \frac{1}{6}}\mu _0^{\prime \prime }(nN_0)^3`$ (76)
$`\mathrm{\Pi }_n`$ $``$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\beta \mu _0^{}(nN_0)^2\right]\left[1{\displaystyle \frac{1}{6}}\beta \mu _0^{\prime \prime }(nN_0)^3\right].`$ (77)
We then get a non-vanishing mean value for $`nN_0`$:

$$(nN_0)_{\mathrm{distor}}=\frac{1}{6}\beta \mu _0^{\prime \prime }(nN_0)^4_{\mathrm{Gauss}}=\frac{1}{2}\beta \mu _0^{\prime \prime }\left(\widehat{P}^2_2\right)^2.$$
(78)
We have therefore predicted in a very simple way the correction to the mean density due to grand canonical fluctuations:

$$\delta \rho _{\mathrm{GC}}(𝐫)=\frac{1}{2}\widehat{P}^2_2\left[_{N_0}^2\rho _0(𝐫;N_0)\frac{\mu _0^{\prime \prime }}{\mu _0^{}}_{N_0}\rho _0(𝐫;N_0)\right].$$
(79)
How does this compare to the general formalism (65)? We need to obtain a partial differential equation for $`\delta \rho _{\mathrm{GC}}`$. We just take the second order derivative of the Gross-Pitaevskii equation (41) with respect to $`N_0`$ and we replace $`\rho _0`$ by $`\sqrt{\rho _0}^2`$ in the resulting equation and in (79). This leads to the remarkable identity

$$\left[\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U+3g_0\rho _0\mu \right]\left(\frac{\delta \rho _{\mathrm{GC}}N_0^1\varphi _a^2\widehat{P}^2_2}{\rho _0^{1/2}}\right)=\widehat{P}^2_2\left[6g_0\rho _0^{1/2}_{N_0}\sqrt{\rho _0}2\mu _0^{}\right]_{N_0}\sqrt{\rho _0}.$$
(80)
The right-hand side of this identity coincides with the source term of (65) involving $`\widehat{P}`$! We have therefore successfully identified $`\delta \rho _{\mathrm{GC}}`$ as a piece of $`\delta \widehat{\rho }_3`$ and we are left with the simpler equation

$`0`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{\Delta }+U+3g_0\rho _0\mu \right]\left({\displaystyle \frac{\delta \widehat{\rho }_3\delta \rho _{\mathrm{GC}}(\widehat{B}_n^{}\widehat{B}_n_2\varphi _0\varphi _a)}{\rho _0^{1/2}}}\right)`$ (81)
$`+`$ $`g_0\rho _0^{1/2}4(\widehat{B}_n^{}\widehat{B}_n\varphi _0\varphi _a)+\widehat{B}_n^2+B_n^2_2.`$
We are not totally satisfied yet since the operator $`\widehat{B}_n`$ does not obey bosonic commutation relations when the system is not spatially homogeneous, in particular the field $`\widehat{B}_n`$ does not commute with itself when taken in two different points:

$`[\widehat{B}_n(𝐫),\widehat{B}_n(𝐫^{})]`$ $`=`$ $`\varphi _a(𝐫)\varphi _0(𝐫^{})\varphi _a(𝐫^{})\varphi _0(𝐫).`$ (82)
$`[\widehat{B}_n(𝐫),\widehat{B}_n^{}(𝐫^{})]`$ $`=`$ $`{\displaystyle \frac{1}{l^D}}\delta _{𝐫,𝐫^{}}\varphi _0(𝐫)\varphi _a(𝐫^{})\varphi _0(𝐫^{})\varphi _a(𝐫).`$ (83)
To circumvent this difficulty we split the field $`\widehat{B}_n`$ in its component along the quasi-condensate mode $`\varphi _0`$ and its orthogonal component:

$$\widehat{B}_n(𝐫)=\widehat{\alpha }\varphi _0(𝐫)+\widehat{\mathrm{\Lambda }}(𝐫).$$
(84)
The Bogoliubov functions $`u_s(𝐫)`$ and $`v_s(𝐫)`$ can be chosen here to be real. The operator $`\widehat{\alpha }`$ can be then written as

$$\widehat{\alpha }=\underset{s}{}\varphi _0|u_s\left(\widehat{b}_s\widehat{b}_s^{}\right)$$
(85)
where we have used the property (66). This clearly shows that the operator $`\widehat{\alpha }`$ is anti-Hermitian:

$$\widehat{\alpha }^{}=\widehat{\alpha }.$$
(86)
The field $`\widehat{\mathrm{\Lambda }}`$ has the following expansion on the $`\widehat{b}_s`$:

$$\widehat{\mathrm{\Lambda }}(𝐫)=\underset{s}{}u_s(𝐫)\widehat{b}_s+v_s(𝐫)\widehat{b}_s^{}$$
(87)
where the index $``$ indicates projection orthogonally to $`\varphi _0`$. This field has now the desired bosonic commutation relations

$`[\widehat{\mathrm{\Lambda }}(𝐫),\widehat{\mathrm{\Lambda }}(𝐫^{})]`$ $`=`$ $`0`$ (88)
$`[\widehat{\mathrm{\Lambda }}(𝐫),\widehat{\mathrm{\Lambda }}^{}(𝐫^{})]`$ $`=`$ $`{\displaystyle \frac{1}{l^D}}\delta _{𝐫,𝐫^{}}\varphi _0(𝐫)\varphi _0(𝐫^{}).`$ (89)
Note that $`\widehat{\alpha }`$ does not commute with $`\widehat{\mathrm{\Lambda }}`$:

$$[\widehat{\alpha },\widehat{\mathrm{\Lambda }}(𝐫)]=\frac{1}{2}\varphi _0(𝐫)\varphi _a(𝐫).$$
(90)
We insert the splitting of $`\widehat{B}_n`$ in (81). The terms quadratic in $`\widehat{\alpha }`$ cancel exactly, in the same way the terms in $`\widehat{Q}^2`$ canceled. The terms linear in $`\widehat{\alpha }`$ can all be expressed in terms of the expectation value of an anticommutator, $`\{\widehat{\alpha },\widehat{\mathrm{\Lambda }}\}_2`$ using the commutation relation (90) and the fact that $`\widehat{\alpha }\widehat{\mathrm{\Lambda }}_2`$ is a real quantity. Furthermore, using the techniques of Appendix E of cartago , as shown in the appendix D, one obtains a simple partial differential equation for the anticommutator:

$$\left[\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U+g_0\rho _0\mu \right]\{\widehat{\alpha },\widehat{\mathrm{\Lambda }}(𝐫)\}_2=\underset{𝐫^{}}{}l^Dg_0\rho _0(𝐫^{})\varphi _0(𝐫^{})\{\widehat{\mathrm{\Lambda }}(𝐫^{})+\widehat{\mathrm{\Lambda }}^{}(𝐫^{}),\widehat{\mathrm{\Lambda }}(𝐫)\}_2.$$
(91)
Remarkably this allows to eliminate completely the operator $`\widehat{\alpha }`$ in (81)! We finally get an equation for the correction to the mean density involving the operator $`\widehat{\mathrm{\Lambda }}`$ only:

$$0=\left[\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U+3g_0\rho _0\mu \right]\left(\frac{\delta \widehat{\rho }_3\delta \rho _{\mathrm{GC}}\widehat{\mathrm{\Lambda }}^{}\widehat{\mathrm{\Lambda }}_2}{\varphi _0}\right)+S(𝐫)$$
(92)
where we have introduced the source term

$$S(𝐫)g_0N_0\varphi _0(𝐫)4\widehat{\mathrm{\Lambda }}^{}(𝐫)\widehat{\mathrm{\Lambda }}(𝐫)+\widehat{\mathrm{\Lambda }}^2(𝐫)+\widehat{\mathrm{\Lambda }}^2(𝐫)\varphi _0^2(𝐫)_2\underset{𝐫^{}}{}l^Dg_0\rho _0(𝐫^{})\varphi _0(𝐫^{})\{\widehat{\mathrm{\Lambda }}(𝐫^{})+\widehat{\mathrm{\Lambda }}^{}(𝐫^{}),\widehat{\mathrm{\Lambda }}(𝐫)\}_2.$$
(93)
It will reveal convenient to introduce the function $`\chi (𝐫)`$ defined in a unique way by

$$0=\left[\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U+3g_0\rho _0\mu \right]\chi (𝐫)+\frac{1}{2}S(𝐫).$$
(94)
We then obtain the following final expression for the correction to the mean density due to the cubic Hamiltonian terms $`H_3`$:

$$\delta \widehat{\rho }_3(𝐫)=\delta \rho _{\mathrm{GC}}(𝐫)+2\varphi _0(𝐫)\chi (𝐫)+\widehat{\mathrm{\Lambda }}^+(𝐫)\widehat{\mathrm{\Lambda }}(𝐫)_2.$$
(95)
In the particular case where the gas is Bose condensed, our general theory for quasi-condensates also applies, of course. One then expects that the result (95) has already been obtained for the condensate and can be given a clear physical interpretation. This expectation is totally justified: as shown in the appendix E, the component of $`\chi (𝐫)/N_0`$ orthogonal to $`\varphi _0`$ is the correction given in yvan2  to the Gross-Pitaevskii condensate wavefunction $`\varphi _0`$ due to the interaction with the non-condensed particles; the component of $`\chi `$ along $`\varphi _0`$ describes the condensate depletion, and $`\widehat{\mathrm{\Lambda }}^+\widehat{\mathrm{\Lambda }}_2`$ is the mean density of non-condensed particles.

## III Applications of the formalism: general formulas

### III.1 Equation of state

What is referred to as the equation of state of the gas is the expression of the chemical potential as function of the mean total number of particles $`N`$ and the temperature $`T`$. It is useful in particular to predict properties of an inhomogeneous gas within the local density approximation.

We therefore have now to calculate $`\mu `$ for the quasi-condensate. This is equivalent to a calculation of $`N_0`$ as $`\mu `$ and $`N_0`$ are by definition related through (43). To lowest order of the theory one assumes a pure quasi-condensate with a density profile $`\rho (𝐫)=\rho _0(𝐫)`$, where $`\sqrt{\rho _0}`$ solves the Gross-Pitaevskii equation (41). One therefore gets $`N=N_0`$ so that $`\mu =\mu _0(N)`$.

The first non-vanishing correction to the density profile is given by (95). By integrating (95) over space we get the corresponding correction for the mean total number of particles:

$`N`$ $``$ $`N_0+\delta N`$ (96)
$`\delta N`$ $``$ $`\delta N_{\mathrm{GC}}+l^D{\displaystyle \underset{𝐫}{}}2\varphi _0(𝐫)\chi (𝐫)+l^D{\displaystyle \underset{𝐫}{}}\widehat{\mathrm{\Lambda }}^{}(𝐫)\widehat{\mathrm{\Lambda }}(𝐫)_2.`$ (97)
The contribution to $`\delta N`$ due to our use of the grand canonical ensemble can be calculated exactly from a spatial integration of (79), using the same technique as in (68):

$$\delta N_{\mathrm{GC}}=k_BT\frac{\mu _0^{\prime \prime }}{2\mu _{0}^{}{}_{}{}^{2}}.$$
(98)
The contribution of the term involving $`\chi `$ can also be made explicit by multiplication of (94) by $`\varphi _a(𝐫)`$ defined in (53) and by spatial integration. The function $`\varphi _a(𝐫)`$ is indeed known lewenstein  to solve the partial differential equation

$$\left[\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U+3g_0\rho _0\mu \right]\varphi _a=N_0\mu _0^{}\varphi _0(𝐫)$$
(99)
which can be checked easily, just by taking the derivative of (41) with respect to $`N_0`$. This leads to

$$l^D\underset{𝐫}{}2\varphi _0(𝐫)\chi (𝐫)=\frac{1}{N_0\mu _0^{}}\underset{𝐫}{}l^D\varphi _a(𝐫)S(𝐫)$$
(100)
where the source term $`S`$ is known explicitly, see (93). We just have now to replace $`N_0`$ by $`N\delta N`$ in (43) and expand to first order in $`\delta N`$:

$$\mu =\mu _0(N\delta N)\mu _0(N)\delta N\mu _0^{}(N_0).$$
(101)
We obtain the following expression for $`\mu `$:

$$\mu \mu _0(N)+k_BT\frac{\mu _0^{\prime \prime }}{2\mu _0^{}}\mu _0^{}(N_0)\underset{𝐫}{}l^D\widehat{\mathrm{\Lambda }}^{}(𝐫)\widehat{\mathrm{\Lambda }}(𝐫)_2+\frac{1}{N_0}\underset{𝐫}{}l^D\varphi _a(𝐫)S(𝐫).$$
(102)
Equivalently we can replace the source term by its explicit expression to get

$`\mu \mu _0(N)+k_BT{\displaystyle \frac{\mu _0^{\prime \prime }}{2\mu _0^{}}}\mu _0^{}(N_0)\left({\displaystyle \frac{1}{2}}+{\displaystyle \underset{𝐫}{}}l^D\widehat{\mathrm{\Lambda }}^{}(𝐫)\widehat{\mathrm{\Lambda }}(𝐫)_2\right)+{\displaystyle \underset{𝐫}{}}l^Dg_0\left(_{N_0}\rho _0(𝐫)\right)\left(2\widehat{\mathrm{\Lambda }}^{}\widehat{\mathrm{\Lambda }}_2+\text{Re}\widehat{\mathrm{\Lambda }}^2_2\right)`$
$`{\displaystyle \underset{𝐫}{}}l^Dg_0\varphi _0^3(𝐫)\{\widehat{\mathrm{\Lambda }}(𝐫)+\widehat{\mathrm{\Lambda }}^{}(𝐫),\widehat{\gamma }\}_2`$ (103)
where we have introduced the operator

$$\widehat{\gamma }=\underset{𝐫}{}l^D\varphi _a(𝐫)\widehat{\mathrm{\Lambda }}(𝐫)$$
(104)
and we have used the identity

$$\mu _0^{}=\underset{𝐫}{}l^Dg_0\varphi _0^2(𝐫)_{N_0}\rho _0(𝐫;N_0)$$
(105)
obtained by performing the scalar product of both sides of (99) with $`\varphi _0`$. Application to spatially homogeneous systems will be given in section IV; in this case both the operator $`\widehat{\gamma }`$ and $`\mu _0^{\prime \prime }`$ vanish.

### III.2 Ground state energy

We now show that the ground state energy of a quasi-condensate can be calculated with exactly the same Bogoliubov formula as for the ground state energy of a condensate.

We have to determine the ground state energy of $`H_2`$. We write it as the expectation value of (51) at zero temperature, that is here in the vacuum of the $`\widehat{b}_s`$ and of $`\widehat{P}`$:

$$\text{ground}(H_2)=l^D\underset{𝐫}{}\widehat{B}^{}\left(\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U+g_0\rho _0\mu \right)\widehat{B}_2+g_0\rho _0\left[\widehat{B}^{}\widehat{B}+\frac{1}{2}\left(\widehat{B}^2+\widehat{B}^2\right)\right]_2.$$
(106)
We reproduce the transformation of subsection II.4. We split $`\widehat{B}`$ in an anomalous part involving $`\widehat{P}`$, $`\widehat{Q}`$, plus the contribution of the antihermitian operator $`\widehat{\alpha }`$ and of $`\widehat{\mathrm{\Lambda }}`$, the orthogonal component of the normal part. In the first expectation value of the right-hand side of (106) the operators $`\widehat{Q}`$ and $`\widehat{\alpha }`$ disappear as they come with the factor $`\varphi _0(𝐫)`$ in $`\widehat{B}`$ and $`\varphi _0`$ solves the Gross-Pitaevskii equation (41). The expectation value of $`\widehat{P}^2`$ in the ground state of $`H_2`$ also vanishes, so that

$$\widehat{B}^{}\left(\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U+g_0\rho _0\mu \right)\widehat{B}_2=\widehat{\mathrm{\Lambda }}^{}\left(\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U+g_0\rho _0\mu \right)\widehat{\mathrm{\Lambda }}_2.$$
(107)
The same transformation is applied to the last expectation value in (106). Remarkably the terms involving $`\widehat{\alpha }`$ exactly cancel when one uses the relations (63), (64), (90) and (236). This leads to

$$\left[\widehat{B}^{}\widehat{B}+\frac{1}{2}\left(\widehat{B}^2+\widehat{B}^2\right)\right]_2=\frac{1}{2}\varphi _0^2+\left[\widehat{\mathrm{\Lambda }}^{}\widehat{\mathrm{\Lambda }}+\frac{1}{2}\left(\widehat{\mathrm{\Lambda }}^2+\widehat{\mathrm{\Lambda }}^2\right)\right]_2.$$
(108)
The expectation values involving $`\widehat{\mathrm{\Lambda }}`$ are readily calculated from the modal expansion (87):

$$\text{ground}(H_2)=\frac{1}{2}\underset{𝐫}{}l^Dg_0\rho _0\varphi _0^2+\underset{s}{}v_s|\left[\left(\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U+2g_0\rho _0\mu \right)|v_s+g_0\rho _0|u_s\right].$$
(109)
As $`(u_s,v_s)`$ is an eigenvector of $`_{\mathrm{GP}}`$, $`(u_s,v_s)`$ is an eigenvector of the operator $``$ defined in yvan2  and this expression can be further simplified to

$$\text{ground}(H_2)=\frac{1}{2}\underset{𝐫}{}l^Dg_0\rho _0\varphi _0^2\underset{s}{}ϵ_sv_s|v_s.$$
(110)
The last step is to include the contribution of $`H_0`$ and to remove the $`\mu \widehat{N}`$ term from the grand canonical Hamiltonian. The ground state energy of the canonical Hamiltonian for $`N`$ particles is therefore

$$E_{\mathrm{ground}}(N)\mu N+E_0(N_0)\mu N_0+\text{ground}(H_2)$$
(111)
where $`E_0`$ is the Gross-Pitaevskii energy (70). As we did in subsection III.1 we replace $`N_0`$ by $`N\delta N`$, where $`\delta N`$ is calculated from $`H_3`$, and we expand $`E_0(N\delta N)`$ to first order in $`\delta N`$:

$$\mu N+E_0(N_0)\mu N_0E_0(N)\delta N(\mu _0(N)\mu )E_0(N).$$
(112)
We recall that by definition $`\mu =\mu _0(N_0)`$. The first term in the right hand side of (110) amounts to performing a small change in the Gross-Pitaevskii energy functional, expressing the fact that a given particle interacts in the gas with $`N1`$ particles so that the mean field term should be proportional to $`N1`$ rather than to $`N`$. The final expression for the ground state energy is:

$$E_{\mathrm{ground}}(N)N\underset{𝐫}{}l^D\left[\frac{\mathrm{}^2}{2m}\varphi _0(𝐫;N)\mathrm{\Delta }\varphi _0(𝐫;N)+U(𝐫)\varphi _0^2(𝐫;N)+\frac{1}{2}g_0(N1)\varphi _0^2(𝐫;N)\right]\underset{s}{}ϵ_sv_s|v_s.$$
(113)
This exactly coincides with the Bogoliubov result, see e.g. equation (71) of yvan2 .

### III.3 Second order correlation function

The second order correlation function of the atomic field is defined as:

$$g_2(𝐫)\widehat{\psi }^{}(𝐫)\widehat{\psi }^{}(\mathrm{𝟎})\widehat{\psi }(\mathrm{𝟎})\widehat{\psi }(𝐫)$$
(114)
where we have taken for simplicity one of the two points as the origin of the coordinates. To calculate $`g_2`$ with the formalism of this paper we have to express $`g_2`$ in terms of the operator $`\widehat{\rho }`$ giving the density. This is achieved thanks to the commutation relation (4) of the bosonic field $`\widehat{\psi }`$:

$$g_2(𝐫)=\widehat{\rho }(𝐫)\widehat{\rho }(\mathrm{𝟎})\frac{\delta _{𝐫,\mathrm{𝟎}}}{l^D}\widehat{\rho }(\mathrm{𝟎}).$$
(115)
We then insert the splitting (32) of $`\widehat{\rho }`$ in terms of the quasi-condensate density $`\rho _0`$ and the fluctuations $`\delta \widehat{\rho }`$:

$$g_2(𝐫)=\rho _0(𝐫)\rho _0(\mathrm{𝟎})+\rho _0(\mathrm{𝟎})\delta \widehat{\rho }(𝐫)+\rho _0(𝐫)\delta \widehat{\rho }(\mathrm{𝟎})+\delta \widehat{\rho }(𝐫)\delta \widehat{\rho }(\mathrm{𝟎})\frac{\delta _{𝐫,\mathrm{𝟎}}}{l^D}\left[\rho _0(\mathrm{𝟎})+\delta \widehat{\rho }(\mathrm{𝟎})\right].$$
(116)
This expression of $`g_2`$ is still exact. We now perform approximations consistent with an expansion of $`g_2`$ up to second order in the small parameters $`ϵ_{1,2}`$. The expectation value of the term quadratic in $`\delta \widehat{\rho }`$ is calculated within the thermal equilibrium for the quadratic Hamiltonian $`H_2`$. The expectation value of $`\delta \widehat{\rho }`$ is evaluated in subsection II.4 by inclusion of the cubic perturbation $`H_3`$. The contribution of $`\delta \widehat{\rho }`$ in the last term of (116) is negligible as it is $`ϵ_1^4`$ times smaller than the leading term in $`g_2`$. We therefore obtain the explicit expression

$$g_2(𝐫)\rho _0(𝐫)\rho _0(\mathrm{𝟎})+\rho _0(\mathrm{𝟎})\delta \widehat{\rho }(𝐫)_3+\rho _0(𝐫)\delta \widehat{\rho }(\mathrm{𝟎})_3+\delta \widehat{\rho }(𝐫)\delta \widehat{\rho }(\mathrm{𝟎})_2\frac{\delta _{𝐫,\mathrm{𝟎}}}{l^D}\rho _0(\mathrm{𝟎}).$$
(117)
This writing is however not the optimal one as the last term in $`1/l^D`$ gives the wrong impression that $`g_2(0)`$ strongly depends on the discretisation length $`l`$ in the continuous limit $`l0`$. In fact this strong dependence exactly compensates a term in $`1/l^D`$ in the density fluctuations $`\delta \widehat{\rho }^2(\mathrm{𝟎})`$ coming from the fact that $`\delta \widehat{\rho }^2(\mathrm{𝟎})`$ is a product of field operators not in normal order. To reveal this fact we express $`\delta \widehat{\rho }`$ in terms of the operator $`\widehat{\mathrm{\Lambda }}`$ of (87)

$$\delta \widehat{\rho }(𝐫)=\sqrt{\rho _0(𝐫)}\left(\widehat{\mathrm{\Lambda }}(𝐫)+\widehat{\mathrm{\Lambda }}^{}(𝐫)\right)+\widehat{P}_{N_0}\rho _0(𝐫;N_0)$$
(118)
and we put the resulting expression in normal order with respect to the field $`\widehat{\mathrm{\Lambda }}`$ thanks to the bosonic commutation relation (89):

$$\delta \widehat{\rho }(𝐫)\delta \widehat{\rho }(\mathrm{𝟎})=:\delta \widehat{\rho }(𝐫)\delta \widehat{\rho }(\mathrm{𝟎}):+\frac{\delta _{𝐫,\mathrm{𝟎}}}{l^D}\rho _0(\mathrm{𝟎})N_0\varphi _0^2(𝐫)\varphi _0^2(\mathrm{𝟎})$$
(119)
where $`::`$ is the standard notation to represent normal order. The spurious term in $`1/l^D`$ is then exactly canceled:

$`g_2(𝐫)`$ $`=`$ $`N_0(N_01)\varphi _0^2(𝐫)\varphi _0^2(\mathrm{𝟎})`$ (120)
$`+`$ $`\rho _0(\mathrm{𝟎})\delta \widehat{\rho }(𝐫)_3+\rho _0(𝐫)\delta \widehat{\rho }(\mathrm{𝟎})_3+:\delta \widehat{\rho }(𝐫)\delta \widehat{\rho }(\mathrm{𝟎}):_2.`$
This expression allows to prove the equivalence with the prediction for $`g_2`$ in the Bogoliubov theory. We do not present the calculations here, as they are a straightforward application of appendix E. Finally we give a last alternative expression to $`g_2`$ equivalent to (120) at the present order:

$$g_2(𝐫)=(11/N)\rho (\mathrm{𝟎})\rho (𝐫)+:\delta \widehat{\rho }(𝐫)\delta \widehat{\rho }(\mathrm{𝟎}):_2$$
(121)
where $`N`$ is the mean total number of particles and $`\rho `$ is the mean total density:

$$\rho (𝐫)=\rho _0(𝐫)+\delta \widehat{\rho }(𝐫)_3.$$
(122)
### III.4 First order correlation function

The first order correlation function of the field is defined as:

$$g_1(𝐫)\widehat{\psi }^{}(𝐫)\widehat{\psi }(\mathrm{𝟎})=\sqrt{\widehat{\rho }(𝐫)}e^{i(\widehat{\theta }(\mathrm{𝟎})\widehat{\theta }(𝐫))}\sqrt{\widehat{\rho }(\mathrm{𝟎})}.$$
(123)
As previously done we perform the calculation up to second order in the small parameters $`ϵ_{1,2}`$. We therefore expand $`\sqrt{\widehat{\rho }}`$ up to second order in $`\delta \widehat{\rho }`$ using (34). Note that we do not expand the exponential in $`\widehat{\theta }(\mathrm{𝟎})\widehat{\theta }(𝐫)`$, contrarily to what we did in the Hamiltonian: as $`𝐫`$ and $`\mathrm{𝟎}`$ are not neighboring points of the lattice anymore, the phase difference of the field can be arbitrarily large. The expansion in $`\delta \widehat{\rho }`$ gives rise to six terms:

$`g_1(𝐫)`$ $`=`$ $`\rho _0^{1/2}(𝐫)\rho _0^{1/2}(\mathrm{𝟎})[e^{i\mathrm{\Delta }\theta }+{\displaystyle \frac{1}{2}}\delta \stackrel{~}{\rho }(𝐫)e^{i\mathrm{\Delta }\theta }+e^{i\mathrm{\Delta }\theta }\delta \stackrel{~}{\rho }(\mathrm{𝟎})`$ (124)
$``$ $`{\displaystyle \frac{1}{8}}\delta \stackrel{~}{\rho }^2(𝐫)e^{i\mathrm{\Delta }\theta }+e^{i\mathrm{\Delta }\theta }\delta \stackrel{~}{\rho }^2(\mathrm{𝟎})2\delta \stackrel{~}{\rho }(𝐫)e^{i\mathrm{\Delta }\theta }\delta \stackrel{~}{\rho }(\mathrm{𝟎})]`$
where we have introduced the following notations to simplify the writing:

$`\mathrm{\Delta }\theta `$ $``$ $`\widehat{\theta }(\mathrm{𝟎})\widehat{\theta }(𝐫)`$ (125)
$`\delta \stackrel{~}{\rho }(𝐫)`$ $``$ $`{\displaystyle \frac{\delta \widehat{\rho }(𝐫)}{\rho _0(𝐫)}}.`$ (126)
We calculate the expectation values in this expression in two steps, first using the thermal equilibrium distribution for $`H_2`$, and then including the corrections due to $`H_3`$.

The thermal expectation values corresponding to the quadratic Hamiltonian $`H_2`$ are evaluated using Wick’s theorem. One first expands the exponential in powers of $`\mathrm{\Delta }\theta `$, one calculates the expectation value of each term, and then one performs an exact resummation of the resulting series. This leads to the simple identities

$`e^{i\mathrm{\Delta }\theta }_2`$ $`=`$ $`e^{(\mathrm{\Delta }\theta )^2_2/2}`$ (127)
$`\delta \stackrel{~}{\rho }(𝐫)e^{i\mathrm{\Delta }\theta }_2`$ $`=`$ $`e^{(\mathrm{\Delta }\theta )^2_2/2}\delta \stackrel{~}{\rho }(𝐫)i\mathrm{\Delta }\theta _2`$ (128)
$`\delta \stackrel{~}{\rho }^2(𝐫)e^{i\mathrm{\Delta }\theta }_2`$ $`=`$ $`e^{(\mathrm{\Delta }\theta )^2_2/2}\left[\delta \stackrel{~}{\rho }^2(𝐫)_2+\left(\delta \stackrel{~}{\rho }(𝐫)i\mathrm{\Delta }\theta _2\right)^2\right]`$ (129)
$`\delta \stackrel{~}{\rho }(𝐫)e^{i\mathrm{\Delta }\theta }\delta \stackrel{~}{\rho }(\mathrm{𝟎})_2`$ $`=`$ $`e^{(\mathrm{\Delta }\theta )^2_2/2}\left[\delta \stackrel{~}{\rho }(𝐫)\delta \stackrel{~}{\rho }(\mathrm{𝟎})_2+\delta \stackrel{~}{\rho }(𝐫)i\mathrm{\Delta }\theta _2i\mathrm{\Delta }\theta \delta \stackrel{~}{\rho }(\mathrm{𝟎})_2\right].`$ (130)
The expectation value of a product of the density fluctuation $`\delta \widehat{\rho }`$ and of the phase variation $`i\mathrm{\Delta }\theta `$ is particularly simple. In a classical field theory this expectation value would obviously vanish, as there is no crossed term in $`H_2`$ between the density fluctuations and the phase. In the present quantum field theory this is not exactly the case as $`\widehat{\rho }`$ and $`\widehat{\theta }`$ do not commute. To show that, we use the fact that the Bogoliubov mode functions $`u_s(𝐫)`$, $`v_s(𝐫)`$ can be chosen to be real, so that $`\delta \widehat{\rho }`$ and $`i\mathrm{\Delta }\theta `$ are linear combinations of $`\widehat{b}_s,\widehat{b}_s^{},\widehat{P}`$ with real coefficients. As a consequence

$$\delta \stackrel{~}{\rho }(𝐫)i\mathrm{\Delta }\theta _2=\delta \stackrel{~}{\rho }(𝐫)i\mathrm{\Delta }\theta _2^{}=i\mathrm{\Delta }\theta \delta \stackrel{~}{\rho }(𝐫).$$
(131)
This leads to

$$\delta \stackrel{~}{\rho }(𝐫)i\mathrm{\Delta }\theta _2=\frac{i}{2}[\delta \stackrel{~}{\rho }(𝐫),\mathrm{\Delta }\theta ]_2=\frac{1\delta _{𝐫,\mathrm{𝟎}}}{2\rho _0(𝐫)l^D}.$$
(132)
The same reasoning can be made for the other expectation value

$$i\mathrm{\Delta }\theta \delta \stackrel{~}{\rho }(\mathrm{𝟎})_2=\frac{1\delta _{𝐫,\mathrm{𝟎}}}{2\rho _0(\mathrm{𝟎})l^D}.$$
(133)
These expressions are second order in $`ϵ_{1,2}`$. An important consequence is that the product of such crossed phase-density expectation values in (129,130) are actually negligible at the present order of the calculation. The resulting form for $`g_1`$, at the level of $`H_2`$, is quite simple:

$$g_1(𝐫)|_{H_2}=\rho _0^{1/2}(𝐫)\rho _0^{1/2}(\mathrm{𝟎})e^{(\mathrm{\Delta }\theta )^2_2/2}\left[1\frac{1}{8}(\mathrm{\Delta }\delta \stackrel{~}{\rho })^2_2+\frac{1}{4}\text{scoria}(𝐫)\right]$$
(134)
The notation $`\mathrm{\Delta }\delta \stackrel{~}{\rho }`$ is similar to the one for the phase:

$$\mathrm{\Delta }\delta \stackrel{~}{\rho }\delta \stackrel{~}{\rho }(\mathrm{𝟎})\delta \stackrel{~}{\rho }(𝐫)$$
(135)
and the ‘scoria’ comes from the crossed expectation value of $`\mathrm{\Delta }\theta `$ and $`\delta \stackrel{~}{\rho }`$:

$$\text{scoria}(𝐫)\left(1\delta _{𝐫,\mathrm{𝟎}}\right)\left(\frac{1}{\rho _0(𝐫)l^D}+\frac{1}{\rho _0(\mathrm{𝟎})l^D}\right).$$
(136)
At this point we face the same apparent problem as in the calculation of $`g_2`$: the ‘scoria’ scales as $`1/l^D`$ and gives the wrong impression that our expression for $`g_1`$ will depend dramatically on $`l`$ in the continuous limit $`l0`$. As in the case of $`g_2`$, we solve this problem by expressing $`\delta \widehat{\rho }`$ and $`\widehat{\theta }`$ in terms of the field $`\widehat{\mathrm{\Lambda }}`$ and putting the operators $`\widehat{\mathrm{\Lambda }}`$, $`\widehat{\mathrm{\Lambda }}^{}`$ in normal order. We use (118) for the expression of $`\delta \widehat{\rho }`$. For the difference of two phase operators, $`\widehat{Q}`$ and the antihermitian operator $`\widehat{\alpha }`$ cancel so that

$$\mathrm{\Delta }\theta =\frac{1}{2i}\left(\mathrm{\Delta }\stackrel{~}{\mathrm{\Lambda }}\mathrm{\Delta }\stackrel{~}{\mathrm{\Lambda }}^{}\right)$$
(137)
where we have introduced the notations

$`\stackrel{~}{\mathrm{\Lambda }}(𝐫)`$ $``$ $`{\displaystyle \frac{\widehat{\mathrm{\Lambda }}(𝐫)}{\rho _0^{1/2}(𝐫)}}`$ (138)
$`\mathrm{\Delta }\stackrel{~}{\mathrm{\Lambda }}`$ $``$ $`\stackrel{~}{\mathrm{\Lambda }}(\mathrm{𝟎})\stackrel{~}{\mathrm{\Lambda }}(𝐫).`$ (139)
After some calculations we arrive at

$`(\mathrm{\Delta }\delta \stackrel{~}{\rho })^2_2`$ $`=`$ $`:(\mathrm{\Delta }\delta \stackrel{~}{\rho })^2:_2+\text{scoria}(𝐫)`$ (140)
$`(\mathrm{\Delta }\theta )^2_2`$ $`=`$ $`:(\mathrm{\Delta }\theta )^2:_2+{\displaystyle \frac{1}{4}}\text{scoria}(𝐫).`$ (141)
As the ‘scoria’ is second order in $`ϵ_{1,2}`$, the exponential function of it can be expanded to first order. We then find as expected that all the $`1/l^D`$ terms exactly cancel:

$$g_1(𝐫)|_{H_2}=\rho _0^{1/2}(𝐫)\rho _0^{1/2}(\mathrm{𝟎})e^{:(\mathrm{\Delta }\theta )^2:_2/2}[1\frac{1}{8}:(\mathrm{\Delta }\delta \stackrel{~}{\rho })^2:_2].$$
(142)
The last step is to include the first order correction to $`g_1`$ coming from the cubic Hamiltonian $`H_3`$. One then has to calculate expectation values with the thermal equilibrium density operator $`\mathrm{exp}[\beta (H_2+H_3)]`$ to first order in $`H_3`$. This thermal density operator can be viewed as the evolution operator during the imaginary time $`i\mathrm{}\beta `$ so that one can use first order time dependent perturbation theory to get

$$\widehat{O}_3\widehat{O}_2_0^\beta 𝑑\tau e^{\tau H_2}H_3e^{\tau H_2}\widehat{O}_2$$
(143)
where $`\widehat{O}`$ is an arbitrary operator of the gas and where we have used the fact that $`H_3`$ has a vanishing expectation value in the thermal equilibrium state for $`H_2`$. One is back to the calculation of expectation values of some operators in the thermal state corresponding to $`H_2`$. Wick’s theorem can be applied. The resulting calculations are very similar to the ones leading to (142), but more involved and are detailed in the appendix F. The same phenomena takes place, that terms of a higher order than the present calculation can be neglected. One then gets

$`e^{i\mathrm{\Delta }\theta }_3`$ $``$ $`e^{(\mathrm{\Delta }\theta )^2_2/2}\left[1+i\mathrm{\Delta }\theta _3\right]`$ (144)
$`\delta \stackrel{~}{\rho }(𝐫)e^{i\mathrm{\Delta }\theta }_3`$ $``$ $`e^{(\mathrm{\Delta }\theta )^2_2/2}\left[\delta \stackrel{~}{\rho }(𝐫)i\mathrm{\Delta }\theta _2+\delta \stackrel{~}{\rho }(𝐫)_3\right].`$ (145)
The first terms in the right-hand side of (144) and (145) already appeared at the level of $`H_2`$, the second terms are corrections due to $`H_3`$ that we now take into account. There is no need to include $`H_3`$ corrections to the other terms of (124) since they are quadratic in $`\delta \widehat{\rho }`$ and are therefore already of second order. The expectation values of the phase $`\widehat{\theta }`$ and of the density fluctuations $`\delta \widehat{\rho }`$ have been calculated in subsection II.4. It was found that the expectation value of the phase operator is space independent so that $`\mathrm{\Delta }\theta _3`$ vanishes. The expectation value of the density fluctuations including the effect of $`H_3`$ was given in (95) and is in general different from zero. Remarkably the whole effect on the correlation function $`g_1`$ of the first order correction in $`H_3`$ is to replace $`\rho _0(𝐫)`$ by the total mean density $`\rho (𝐫)`$ defined in (122)!

We write our final expression for the first order correlation function of the field, calculated consistently up to $`ϵ_{1,2}^2`$:

$$g_1(𝐫)=\sqrt{\rho (𝐫)\rho (\mathrm{𝟎})}\mathrm{exp}[\frac{1}{2}:(\mathrm{\Delta }\theta )^2:_2\frac{1}{8}:(\mathrm{\Delta }\delta \stackrel{~}{\rho })^2:_2].$$
(146)
Note that we have inserted the contribution of the density fluctuations inside the exponential factor, which is allowed at the order of the present calculation since this contribution is of order $`ϵ_{1,2}^2`$.

What happens in the regime where a true condensate is present ? Both phase and density fluctuations are small, so that the exponential function in (146) can be expanded to first order. We then express $`\mathrm{\Delta }\theta `$ and $`\mathrm{\Delta }\stackrel{~}{\rho }`$ in terms of the operator $`\widehat{\mathrm{\Lambda }}`$ and of the operator $`\widehat{P}`$. Since the Bogoliubov theory is usually considered in the canonical ensemble we remove the terms corresponding to the grand canonical fluctuations of the particle number. We then recover exactly the Bogoliubov prediction:

$$g_1^{\mathrm{Bog}}(𝐫)=\mathrm{\Psi }_c(𝐫)\mathrm{\Psi }_c(\mathrm{𝟎})+\widehat{\mathrm{\Lambda }}^{}(𝐫)\widehat{\mathrm{\Lambda }}(\mathrm{𝟎})$$
(147)
where $`\mathrm{\Psi }_c(𝐫)=\sqrt{N_0}\varphi _0(𝐫)+\chi (𝐫)/\sqrt{N_0}`$ is the condensate field. Amazingly the general formula (146) for quasi-condensates can be related to the Bogoliubov formula in the following very simple way:

$$g_1(𝐫)=\sqrt{\rho (𝐫)\rho (\mathrm{𝟎})}\mathrm{exp}\left[\frac{g_1^{\mathrm{Bog}}(𝐫)}{\sqrt{\rho (𝐫)\rho (\mathrm{𝟎})}}1\right].$$
(148)
## IV Explicit results for the spatially homogeneous case

In this section we apply our approach to a spatially homogeneous Bose gas. The quasi-condensate density is then uniform:

$$\rho _0(𝐫)=\frac{N_0}{L^D}=\frac{\mu }{g_0}.$$
(149)
The Bogoliubov equations (50) can then be exactly solved for any dimension of space and lead to $`u_k(r)=\overline{u_k}e^{ikr}/L^{D/2}`$ and $`v_k(r)=\overline{v_k}e^{ikr}/L^{D/2}`$ with:

$$\overline{u_k}\overline{v_k}=\left[\frac{\frac{\mathrm{}^2k^2}{2m}+2\mu }{\frac{\mathrm{}^2k^2}{2m}}\right]^{1/4}\text{and}\overline{u_k}+\overline{v_k}=\left[\frac{\frac{\mathrm{}^2k^2}{2m}}{\frac{\mathrm{}^2k^2}{2m}+2\mu }\right]^{1/4}.$$
(150)
The corresponding eigenenergies are given by:

$$ϵ_k=\left[\frac{\mathrm{}^2k^2}{2m}(\frac{\mathrm{}^2k^2}{2m}+2\mu )\right]^{1/2}.$$
(151)
### IV.1 Equation of state

From the general expression (103) for the chemical potential of the gas we arrive in the thermodynamical limit at

$$\mu =\rho g_0+g_0_𝒟\frac{d𝐤}{(2\pi )^D}[(\overline{u}_k+\overline{v}_k)^2n_k+\overline{v}_k(\overline{u}_k+\overline{v}_k)]$$
(152)
where $`n_k=1/(\mathrm{exp}(\beta ϵ_k)1)`$ is the mean occupation number of the Bogoliubov mode $`𝐤`$. $`𝒟=[\pi /l,\pi /l]^D`$ is the square domain of integration in the $`𝐤`$ space. The integral over the wavevector $`𝐤`$ does not contain any infrared divergence for any dimension of space. However the long wavevector behavior given by

$$\overline{v}_k(\overline{u}_k+\overline{v}_k)\frac{m\mu }{\mathrm{}^2k^2}$$
(153)
gives rise to an integral convergent in $`1D`$ and divergent in $`2D`$ and $`3D`$ in the $`l0`$ limit. This gives the impression that the result depends strongly on $`l`$. The solution of this paradox comes from the link between the bare coupling constant $`g_0`$ of the model potential in the discretized space and the low energy two-body scattering properties of the exact potential in the continuous space. This gives to $`g_0`$ in two and three dimensions a dependence in $`l`$ so that our expression for $`\mu `$ does not depend on $`l`$ anymore in the $`l0`$ limit. In one dimension, the bare coupling $`g_0`$ is simply equal to the actual coupling strength $`g`$ for $`l0`$ and there is no divergence. A $`T=0`$, the equation (152) leads to:

$$\mu =g\rho \left(1\frac{1}{\pi \rho \xi }\right)$$
(154)
where $`\xi `$ is the healing length defined in (23). This agrees with the result of Lieb-Liniger in the weak interaction limit lieb . In three dimensions, we refer to the appendix of cartago  where the calculation has been done. One finds:

$$g_0=\frac{g}{1g{\displaystyle _𝒟}{\displaystyle \frac{d𝐤}{(2\pi )^3}}{\displaystyle \frac{1}{\mathrm{}^2k^2/m}}}.$$
(155)
$`g`$ is the usual $`3D`$ coupling strength given by :

$$g=\frac{4\pi \mathrm{}^2a}{m}$$
(156)
where $`a`$ is the exact potential scattering length. A more explicit form of (155) is

$$g_0=\frac{g}{1Ka/l}$$
(157)
where $`K=2.442\mathrm{}`$. It has to be noted that the difference between $`g_0`$ and $`g`$ is still small in the validity domain of our approach since it is a second order correction in $`ϵ_{1,2}`$: taking $`l\xi `$ one finds $`a/l1/\rho l^3`$. Replacing the first factor $`g_0`$ in (152) with the formula (155) expanded up to second order in $`ϵ_{1,2}`$ gives:

$$\mu =\rho g+g_0_𝒟\frac{d𝐤}{(2\pi )^3}\left((\overline{u}_k+\overline{v}_k)^2n_k+\overline{v}_k(\overline{u}_k+\overline{v}_k)+\frac{m\mu }{\mathrm{}^2k^2}\right)$$
(158)
One can then safely take the $`l0`$ limit. At $`T=0`$, the integration gives:

$$\mu =g\rho \left(1+\frac{32\sqrt{\pi }}{3}\sqrt{\rho a^3}\right)$$
(159)
which is the same result as Lee and Yang’s Leeyang . In two dimensions, the low energy two-body scattering of a general short range potential is described by a single length $`a`$ also named the scattering length. In a continuous space, the $`T`$ matrix can be calculated in the low energy limit:

$$𝐤|T(E+i\eta )|𝐤^{}\frac{2\pi \mathrm{}^2}{m\left[\mathrm{ln}\left(\frac{ak_0}{2}\right)+Ci\pi /2\right]}$$
(160)
where $`C=.57721\mathrm{}`$ is the Euler constant, $`a`$ is the scattering length, $`E=\mathrm{}^2k_0^2/m`$ and $`\eta 0^+`$. We can also calculate the $`T`$ matrix for the discrete delta potential defined by Eq. (27) which can also be expressed as:

$$V=\frac{g_0}{l^2}|𝐫=0𝐫=0|.$$
(161)
The general scattering theory gives the relations between the $`T`$ matrix, the propagator $`G`$ and the free propagator $`G_0`$:

$`T`$ $`=`$ $`V+VGV`$ (162)
$`G`$ $`=`$ $`G_0+G_0VG`$ (163)
Using these relations and Eq. (161) for the potential, we find:

$$𝐤|T_{\mathrm{grid}}(E+i\eta )|𝐤^{}=\frac{g_0}{1g_0𝐫=0|G_0(E+i\eta )|𝐫=0}.$$
(164)
The only term we need to calculate is the free propagator taken at the origin which is conveniently performed with a Fourier transform:

$$𝐫=0|G_0(E+i\eta )|𝐫=0=_𝒟\frac{d𝐤}{(2\pi )^2}\frac{1}{E+i\eta \mathrm{}^2k^2/m}.$$
(165)
We split the square $`𝒟`$ into a disk of radius $`\pi /l`$ and the complementary domain. Integration over the complementary domain gives simply a constant term in the low energy limit $`E\mathrm{}^2/ml^2`$:

$$J\frac{2\pi \mathrm{}^2}{m}_{𝒟\mathrm{disk}}\frac{d𝐤}{(2\pi )^2}\frac{1}{E+i\eta \mathrm{}^2k^2/m}\frac{1}{2\pi }_{𝒟\mathrm{disk}}\frac{d𝐤}{k^2}=\frac{2G}{\pi }\mathrm{ln}(2)$$
(166)
where $`G=.91596`$ is the Catalan constant. The disk integration is straightforward and leads to the following expression for the $`T`$ matrix:

$$𝐤|T_{\mathrm{grid}}(E+i\eta )|𝐤^{}=\frac{1}{{\displaystyle \frac{1}{g_0}}{\displaystyle \frac{m}{2\pi \mathrm{}^2}}\mathrm{ln}\left({\displaystyle \frac{lk_0}{\pi }}\right)+{\displaystyle \frac{im}{4\mathrm{}^2}}{\displaystyle \frac{m}{2\pi \mathrm{}^2}}J}.$$
(167)
We now take $`T_{\mathrm{grid}}=T`$, where $`T`$ is approximated by (160), in order to reproduce the low-energy scattering properties of the exact potential. This leads to:

$$\frac{1}{g_0}=\frac{m}{2\pi \mathrm{}^2}\left(\mathrm{ln}\left(\frac{l}{\pi a}\right)C+\frac{2G}{\pi }\right).$$
(168)
Note that the condition (26) has to be satisfied in our approach. In two dimensions, this gives $`\mathrm{}^2/mg_01`$ or using Eq. (168):

$$\frac{1}{2\pi }\mathrm{ln}\left(\frac{l}{a}\right)1.$$
(169)
We now show that the logarithmic dependence in $`l`$ appearing in $`g_0`$, Eq. (168), exactly cancels the one appearing in the equation of state. Eq. (152) can be rewritten as:

$$\rho =\frac{\mu }{g_0}_𝒟\frac{d𝐤}{(2\pi )^2}\left[(\overline{u_k}+\overline{v_k})^2n_k+\overline{v_k}(\overline{u_k}+\overline{v_k})\right].$$
(170)
In the thermal part, one can immediately take the $`l0`$ limit. In order to calculate the integral corresponding to the $`T=0`$ case, we use the same technique as for the calculation of $`g_0`$: the integration is made on a disk domain of radius $`\pi /l`$ and we keep as a correction the integration over the complementary domain. The complementary domain integration is made by using the high wavevector behavior of $`\overline{v_k}(\overline{u_k}+\overline{v_k})`$, Eq. (153). This leads to:

$$_𝒟\frac{d𝐤}{(2\pi )^2}\overline{v_k}(\overline{u_k}+\overline{v_k})=\frac{m\mu }{4\pi \mathrm{}^2}\left[\mathrm{ln}\left(\frac{\pi ^2\mathrm{}^2}{ml^2\mu }\right)12J\right].$$
(171)
Using Eqs.(168) and (171) in Eq. (170), we arrive at an implicit equation of state:

$$\rho =\frac{m\mu }{4\pi \mathrm{}^2}\mathrm{ln}\left(\frac{4\mathrm{}^2}{a^2m\mu e^{2C+1}}\right)\frac{d𝐤}{(2\pi )^2}(\overline{u_k}+\overline{v_k})^2n_k.$$
(172)
Remarkably this is identical to the result (20.45) obtained by the functional integral method in popov . A $`T=0`$, one can show from the condition $`\rho \xi ^21`$, see Eq.(26), that the validity condition of our approach is $`\mathrm{ln}(1/\rho a^2)4\pi `$. If one approximately inverts Eq.(172) neglecting constant terms and $`\mathrm{ln}(\mathrm{ln}(1/\rho a^2))`$ with respect to $`\mathrm{ln}(1/\rho a^2)`$, one recovers Schick’s formula schick .

### IV.2 Are density and gradient-of-phase fluctuations small ?

As mentioned in section II.2, our approach relies in particular on two assumptions: the assumption that the relative density fluctuation $`ϵ_1`$ is small, and the assumption that the phase variation $`ϵ_2`$ between two neighboring points of the grid is small.

Let us consider first the relative density fluctuations. Thanks to (119) their mean square value can be separated in two parts:

$$ϵ_1^2=\frac{\delta \widehat{\rho }^2(\mathrm{𝟎})_2}{\rho _0^2}=\frac{1}{\rho _0l^D}+\frac{:\delta \widehat{\rho }^2(\mathrm{𝟎}):_2}{\rho _0^2}$$
(173)
where we have neglected $`1/L^D`$ with respect to $`1/l^D`$ in the thermodynamical limit. The second term in (173), involving the normal order, is expressed in terms of the $`\overline{u}_k,\overline{v}_k`$ in the thermodynamical limit as

$$\frac{:\delta \widehat{\rho }^2(\mathrm{𝟎}):_2}{\rho _0^2}=\frac{2}{\rho _0}_𝒟\frac{d𝐤}{(2\pi )^D}\left[(\overline{u}_k+\overline{v}_k)^2n_k+\overline{v}_k(\overline{u}_k+\overline{v}_k)\right]$$
(174)
where the integration domain is $`𝒟=[\pi /l,\pi /l]^D`$. At zero temperature one introduces the change of variable $`𝐪=𝐤\xi `$ in the integral: one finds that (174) is of the order of $`1/\rho _0\xi `$ in 1D, of the order of $`\mathrm{ln}(\xi /l)/\rho _0\xi ^2`$ in 2D and of the order of $`1/\rho _0\xi ^2l`$ in 3D. Since $`l<\xi `$ the second term in (173) is dominated by the first term and one has indeed

$$ϵ_1^2\frac{1}{\rho _0l^D}.$$
(175)
At finite temperature we have to calculate the thermal contribution to (174) involving the occupation number $`n_k`$.

At a temperature $`k_BT<\mu `$ we use the low momentum expansion of the $`\overline{u}_k+\overline{v}_k`$ and $`ϵ_k`$ and we find that the thermal contribution is $`(k_BT/\mu )^{D+1}(l/\xi )^D`$ times smaller than $`1/\rho _0l^D`$.

At a temperature $`k_BT>\mu `$, that is $`\lambda <\xi `$, the treatment depends on the dimension of space. In 1D the main contribution to the integral comes from the domain $`ϵ_k\mu `$, over which one can approximate the Bose formula by its low energy limit $`k_BT/ϵ_k`$. This leads to a normal ordered fluctuation (174) of the order of $`k_BT/(\mu \rho _0\xi )`$. This is larger than $`1/\rho _0\lambda `$ so that the condition $`l<\lambda `$ then no longer implies that the first term $`1/\rho _0l`$ in (173) is the dominant one. For convenience one can however adjust $`l`$ to a value such that

$$\frac{1}{\rho _0l}\frac{k_BT}{\mu }\frac{1}{\rho _0\xi }.$$
(176)
The condition for weak density fluctuations then becomes

$$ϵ_1^2\frac{k_BT}{\mu }\frac{1}{\rho _0\xi }1.$$
(177)
Using $`\rho _0\rho `$ and $`\mu g\rho `$ we recover a condition already obtained in yvan1  with a pure classical field approach. Note that this condition can be rewritten as $`\xi l_c`$ where the coherence length of the field will be defined in (187). In 2D both the low energy domain $`ϵ_k<k_BT`$ and the high energy domain $`ϵ_k>k_BT`$ have an important contribution. In the low energy domain we approximate the Bose law by its low energy limit. In the high energy domain we keep the full Bose law but $`ϵ_k`$ being then larger than $`\mu `$, we approximate $`\overline{u}_k+\overline{v}_k`$ by unity and $`ϵ_k`$ by $`\mathrm{}^2k^2/2m`$. This leads to a normal ordered fluctuation (174) of the order of $`\mathrm{ln}(k_BT/\mu )k_BT/(\mu \rho _0\xi ^2)`$, a quantity which is larger than $`1/\rho _0\lambda ^2`$. As in 1D we therefore adjust $`l`$ so that

$$ϵ_1^2\frac{1}{\rho _0l^2}\frac{k_BT}{\mu }\mathrm{ln}\left(\frac{k_BT}{\mu }\right)\frac{1}{\rho _0\xi ^2}.$$
(178)
In 3D the high energy domain $`ϵ_k>k_BT`$ gives the dominant contribution so that the normal ordered expectation value (174) scales as $`1/\rho _0\lambda ^3`$. This is dominated by the first term in (173) so that the estimate (175) applies as soon as $`l<\lambda ,\xi `$.

Let us consider now the condition that the mean squared phase change over a grid cell:

$$ϵ_2^2=(l\widehat{\theta })^2_2=\frac{l^2}{2\rho _0}_𝒟\frac{d𝐤}{(2\pi )^D}k^2(\overline{u}_k\overline{v}_k)^2(n_k+1/2)$$
(179)
is much smaller than unity. The presence of the factor $`k^2`$ inside the integral, coming from the action of $``$, has the consequence that the contribution to the integral is dominated by the high energy domain. At zero temperature one can replace $`\overline{u}_k\overline{v}_k`$ by unity since the integral is dominated by wavevectors $`k1/l>1/\xi `$. This leads to

$$ϵ_2^2\frac{1}{\rho _0l^D}$$
(180)
whatever the dimension $`D`$.

At a temperature $`k_BT<\mu `$ we estimate the thermal contribution by replacing $`\overline{u}_k\overline{v}_k`$ and $`ϵ_k`$ by their low momentum approximations: the thermal contribution is then $`(l/\lambda )^{2+D}(\xi /\lambda )^D`$ times smaller than the zero temperature result (180) and is therefore negligible since $`l<\xi <\lambda `$.

At a temperature $`k_BT>\mu `$ we use the high energy approximation replacing $`\overline{u}_k\overline{v}_k`$ by unity and $`ϵ_k`$ by $`\mathrm{}^2k^2/2m`$. Note that this works even in 1D thanks to the presence of the $`k^2`$ factor in the integral (179). This leads to a thermal contribution which is $`(l/\lambda )^{2+D}`$ times smaller than the zero temperature contribution (180), and which is negligible since $`l<\lambda `$.

We conclude that the small parameter $`ϵ_2`$ of the theory, ensuring that there is a weak phase variation over a grid cell, is always given by (180) provided that the conditions (22), (25) and (26) are satisfied.

One may wonder if the condition that the corrections of the mean density due to the interaction $`H_3`$ between the Bogoliubov modes lead to an extra validity condition of our treatment. For the considered case of a spatially homogeneous gas it turns out that the answer to this question is no. One has indeed the remarkable identity in the thermodynamical limit:

$$\frac{1}{\rho _0}\delta \widehat{\rho }(𝐫)_3=\frac{1}{2\rho _0^2}:\delta \widehat{\rho }(𝐫)^2:_2.$$
(181)
If the relative density fluctuations are weak, the relative correction to the density will be also weak.

To end this subsection we discuss briefly the second order correlation function of the field $`g_2(𝐫)`$. Restricting the general formula (121) to the spatially homogeneous case in the thermodynamical limit we obtain

$$g_2(𝐫)=\rho ^2+2\rho _𝒟\frac{d𝐤}{(2\pi )^D}\left[(\overline{u}_k+\overline{v}_k)^2n_k+\overline{v}_k(\overline{u}_k+\overline{v}_k)\right]\mathrm{cos}(𝐤𝐫).$$
(182)
Limiting cases of this general formula can be compared to existing results in the literature. At zero temperature for a 1D Bose gas one gets for $`r=0`$:

$$g_2(0)=\rho ^2\left(1\frac{2}{\pi \rho \xi }\right).$$
(183)
This formula can be checked from lieb : the mean interaction energy per particle $`v`$ is equal to $`g_2(0)`$ multiplied by $`g/2\rho `$, and $`v`$ can be calculated in the weakly interacting regime by combining (3.29) of lieb  (relating $`v`$ to the derivative of the ground state energy with respect to $`g`$) and (4.2) of lieb  (giving the ground state energy in the Bogoliubov approximation). This exactly leads to (183). This prediction for $`g_2(0)`$ also appears in shlyap\_gangardt .

### IV.3 First order correlation function

Thanks to the general formula (148) the first order correlation function of the field for the quasi-condensate is immediately related to the one of the Bogoliubov theory, here in the thermodynamical limit:

$$\mathrm{ln}[g_1(𝐫)/\rho ]=\frac{g_1^{\mathrm{Bog}}(𝐫)}{\rho }1=\frac{1}{\rho }\frac{d^Dk}{(2\pi )^D}\left[\left(\overline{u}_k^2+\overline{v}_k^2\right)n_k+\overline{v}_k^2\right](1\mathrm{cos}𝐤𝐫).$$
(184)
We have also taken here the continuous limit $`l0`$, which does not lead to any divergence.

We concentrate our analysis to the 1D case and we make the link with existing results in the literature. These existing results deal with the asymptotic behavior of $`g_1`$ for large $`r`$, where $`r`$ is the absolute value of the spatial coordinate. At zero temperature, we find for $`r\xi `$:

$$g_1(r)\rho \left(\frac{r_1}{r}\right)^{1/2\pi \rho \xi }$$
(185)
with $`r_1=e^{2C}\xi /41.037\xi `$ where $`C=0.57721\mathrm{}`$ is Euler’s constant note\_ci . This reproduces a result obtained in a non explicit way by popovR . At a finite temperature, $`g_1/\rho `$ is the exponential of an integral of the form

$$_0^+\mathrm{}\frac{A(k)}{k^2}[1\mathrm{cos}(kr)]$$
(186)
where the function $`A(k)`$ is a regular and even function of $`k`$, therefore behaving quadratically with $`k`$ around $`k=0`$ note\_A . Writing $`A(k)`$ as $`[A(k)A(0)]+A(0)`$ and splitting accordingly the integral one obtains for $`r`$ much larger than both $`\xi `$ and $`\lambda `$:

$$\mathrm{ln}[g_1(r)/\rho ]=\frac{r}{l_c}+K+o(1/r^n)$$
(187)
where the coherence length $`l_c=\rho \lambda ^2/\pi `$ coincides with the one of schwartz2  and the constant $`K`$ is given by

$$K=_0^+\mathrm{}\frac{A(k)A(0)}{k^2}.$$
(188)
Since $`A(k)`$ is even one can show by repeated integration by parts that the remainder in (187) tends to 0 faster than any power law, contrarily to what is stated in schwartz2 .

Of course our formula gives access to $`g_1`$ for any value of the distance. This is illustrated in figure 1 where we have plotted the logarithm of $`g_1`$ as function of $`r/\xi `$ for various temperatures.

As a consequence we can for example calculate the momentum distribution of the atoms:

$$\mathrm{\Pi }(p)=2_0^+\mathrm{}𝑑rg_1(r)\mathrm{cos}(pr/\mathrm{})$$
(189)
normalized here as $`𝑑p\mathrm{\Pi }(p)=2\pi \mathrm{}\rho `$ so that $`\mathrm{\Pi }(p)`$ is dimensionless. This is illustrated in figure 2 where we have plotted the momentum distribution for various temperatures and for $`\rho \xi =10`$. Using integration by part we can show that the behavior of $`\mathrm{\Pi }`$ for large $`p`$ is related to the fact that the third order derivative of $`g_1`$ in $`r=0^+`$ does not vanish:

$$\mathrm{\Pi }(p)\frac{2\mathrm{}^4g_1^{(3)}(0^+)}{p^4}\text{with}g_1^{(3)}(0^+)=\mu ^2m^2/(2\mathrm{}^4).$$
(190)
This prediction, valid at zero or finite temperature, agrees with the weak interaction limit of a recently obtained exact result based on the Bethe ansatz olshanii . At zero temperature we find that the momentum distribution diverges in $`p=0`$ as

$$\mathrm{\Pi }(p)\frac{\mathrm{}\rho \pi \nu }{p}(r_1p/\mathrm{})^\nu $$
(191)
where $`\nu =1/(2\pi \rho \xi )1`$.

## V Conclusion

We have studied the thermal equilibrium of weakly interacting degenerate Bose gases in the regime of weak density fluctuations, the so-called quasi-condensate regime. The method can be considered as a Bogoliubov method in the density-phase representation of the field operator.

In a first step one discretizes the real space in cells of size $`l`$: $`l`$ is small enough that the macroscopic properties of the gas are not affected by the discretisation, and $`l`$ is large enough that each cell contains on the average a large number of particles. The macroscopic occupation of each cell allows to give a precise definition of the phase operator, as pioneered by Girardeau girardeau .

In a second step one performs a systematic expansion of the full Hamiltonian in terms of two small parameters, the relative density fluctuations inside a cell and the phase change over a grid cell. This procedure leads to an exact expansion of the observables of the gas in the regime of weak interactions and low density fluctuations, in 1D, 2D and 3D. In particular it is free of any ultraviolet or infrared divergences and exactly matches the usual Bogoliubov predictions when the gas contains a true Bose-Einstein condensate.

As a first application of the general formalism, we have given in this paper formulas for the equation of state of the gas, the ground state energy, the first order and second order correlation functions of the field. We have applied these formulas to the spatially homogeneous case in 1D, 2D and 3D, recovering in this way known results, but obtaining also to our knowledge new results, like the full position dependence of the first order correlation function of the field.

###### Acknowledgements.

We acknowledge useful discussions with Iacopo Carusotto, Gora Shlyapnikov, Dimitri Gangardt, Gordon Baym. Laboratoire Kastler Brossel is a research unit of Ecole normale supérieure and Université Paris 6, associated to CNRS.

## Appendix A Expansion of the Hamiltonian

As explained in II.2 we expand the Hamiltonian (28) up to third order in powers of the small parameters $`ϵ_1`$ and $`ϵ_2`$ defined in (33), (35). This will produce terms $`H^{(n_1,n_2)}`$ of order $`ϵ_1^{n_1}ϵ_2^{n_2}`$ with $`n_1+n_23`$. The expansion of the potential energy part $`H_{\mathrm{pot}}`$ defined in (30) is very simple as it only involves the operator giving the density. The only point is to realize that the term $`1/l^D`$ is $`ϵ_1^2`$ times smaller than the zeroth order density $`\rho _0`$. This leads to

$`H_{\mathrm{pot}}^{(0,0)}`$ $`=`$ $`{\displaystyle \underset{𝐫}{}}l^D\rho _0\left[U(𝐫)\mu +{\displaystyle \frac{g_0}{2}}\rho _0\right]`$ (192)
$`H_{\mathrm{pot}}^{(1,0)}`$ $`=`$ $`{\displaystyle \underset{𝐫}{}}l^D\delta \widehat{\rho }\left[U(𝐫)\mu +g_0\rho _0\right]`$ (193)
$`H_{\mathrm{pot}}^{(2,0)}`$ $`=`$ $`{\displaystyle \underset{𝐫}{}}l^D{\displaystyle \frac{g_0}{2}}\left[\delta \widehat{\rho }^2{\displaystyle \frac{\rho _0}{l^D}}\right].`$ (194)
$`H_{\mathrm{pot}}^{(3,0)}`$ $`=`$ $`{\displaystyle \underset{𝐫}{}}{\displaystyle \frac{g_0}{2}}\delta \widehat{\rho }.`$ (195)
The expansion of the kinetic energy part (31) is more complicated as it involves also the phase operator $`\widehat{\theta }`$, which furthermore does not commute with $`\delta \widehat{\rho }`$. An expression slightly more convenient than (31) can be given for the kinetic energy. Thanks to the periodic boundary conditions one can freely shift the summation variable in the term of (31) involving $`\widehat{\rho }_j`$, so that

$$H_{\mathrm{kin}}=\frac{\mathrm{}^2}{2ml^2}\underset{𝐫,j}{}l^D\left\{\left[\sqrt{\widehat{\rho }}e^{i(\widehat{\theta }_{+j}\widehat{\theta })}\sqrt{\widehat{\rho }_{+j}}+\text{h.c.}\right]2\sqrt{\widehat{\rho }}\right\}.$$
(196)
The calculation to zeroth order in $`ϵ_2`$ can be done first easily: using the expansion (37) to zeroth order, we get from (31) to all orders in $`ϵ_1`$:

$$H_{\mathrm{kin}}^{(+\mathrm{},0)}=\frac{\mathrm{}^2}{2m}\underset{𝐫}{}l^D\sqrt{\widehat{\rho }}\mathrm{\Delta }\sqrt{\widehat{\rho }}.$$
(197)
It involves a function of $`\widehat{\rho }`$ only that it is easily expanded in powers of $`ϵ_1`$ using (34). A simplification occurs after summation over the lattice, as the matrix $`\mathrm{\Delta }`$ is symmetric for the considered periodic boundary conditions:

$$\underset{𝐫}{}u\mathrm{\Delta }v=\underset{𝐫}{}(\mathrm{\Delta }u)v$$
(198)
where $`u`$ and $`v`$ are arbitrary functions on the lattice. This leads to

$`H_{\mathrm{kin}}^{(0,0)}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}l^D{\displaystyle \underset{𝐫}{}}\sqrt{\rho _0}\mathrm{\Delta }\sqrt{\rho _0}`$ (199)
$`H_{\mathrm{kin}}^{(1,0)}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}l^D{\displaystyle \underset{𝐫}{}}{\displaystyle \frac{\delta \widehat{\rho }}{\sqrt{\rho _0}}}\mathrm{\Delta }\sqrt{\rho _0}`$ (200)
$`H_{\mathrm{kin}}^{(2,0)}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}l^D{\displaystyle \underset{𝐫}{}}\left[{\displaystyle \frac{\delta \widehat{\rho }}{4\sqrt{\rho _0}}}\mathrm{\Delta }{\displaystyle \frac{\delta \widehat{\rho }}{\sqrt{\rho _0}}}{\displaystyle \frac{\delta \widehat{\rho }^2}{4\rho _0^{3/2}}}\mathrm{\Delta }\sqrt{\rho _0}\right]`$ (201)
$`H_{\mathrm{kin}}^{(3,0)}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}l^D{\displaystyle \underset{𝐫}{}}\left[{\displaystyle \frac{1}{8}}{\displaystyle \frac{\delta \widehat{\rho }^3}{\rho _0^{5/2}}}\mathrm{\Delta }\sqrt{\rho _0}{\displaystyle \frac{1}{8}}{\displaystyle \frac{\delta \widehat{\rho }^2}{\rho _0^{3/2}}}\mathrm{\Delta }{\displaystyle \frac{\delta \widehat{\rho }}{\rho _0^{1/2}}}\right].`$ (202)
The second order term of vanishing order in $`ϵ_1`$ is also immediately obtained:

$$H_{\mathrm{kin}}^{(0,2)}=\frac{\mathrm{}^2}{2ml^2}l^D\underset{𝐫,j}{}\sqrt{\rho _0\rho _{0,+j}}(\widehat{\theta }_{+j}\widehat{\theta })^2.$$
(203)
The last second order quantity to calculate is $`H_{\mathrm{kin}}^{(1,1)}`$, which is first order in $`ϵ_1`$ and first order in $`ϵ_2`$. There are four terms, two involving $`\widehat{\theta }_{+j}`$ and two being their hermitian conjugate. One can then collect the terms to form commutators:

$`H_{\mathrm{kin}}^{(1,1)}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2ml^2}}l^D{\displaystyle \underset{𝐫,j}{}}{\displaystyle \frac{i}{2}}\left({\displaystyle \frac{\rho _0}{\rho _{0,+j}}}\right)^{1/2}[\widehat{\theta }_{+j}\widehat{\theta },\delta \widehat{\rho }_{+j}]{\displaystyle \frac{i}{2}}\left({\displaystyle \frac{\rho _{0,+j}}{\rho _0}}\right)^{1/2}[\widehat{\theta }_{+j}\widehat{\theta },\delta \widehat{\rho }].`$ (204)
$`=`$ $`{\displaystyle \frac{\mathrm{}^2}{4ml^2}}{\displaystyle \underset{𝐫,j}{}}\left[\left({\displaystyle \frac{\rho _{0,+j}}{\rho _0}}\right)^{1/2}+\left({\displaystyle \frac{\rho _0}{\rho _{0,+j}}}\right)^{1/2}\right]`$ (205)
where we have used the commutation relation of $`\widehat{\rho }`$ and $`\widehat{\theta }`$, see (20).

We collect all the second order c-number contributions to the Hamiltonian $`\overline{H}`$ in a single energy functional of the density profile of the quasi-condensate,

$$E_2[\rho _0]=\frac{g_0}{2}\underset{𝐫}{}\rho _0\frac{\mathrm{}^2}{4ml^2}\underset{𝐫,j}{}\left[\left(\frac{\rho _{0,+j}}{\rho _0}\right)^{1/2}+\left(\frac{\rho _0}{\rho _{0,+j}}\right)^{1/2}\right].$$
(206)
The technique used to calculate $`H_{\mathrm{kin}}^{(1,1)}`$ can be extended to the calculation of $`H_{\mathrm{kin}}^{(2,1)}`$. There are now three terms and their hermitian conjugates. Two of these terms, when combined with their hermitian conjugate, form a commutator that is calculated according to (20). The third term and its hermitian conjugate involve the expression

$$\delta \widehat{\rho }(\widehat{\theta }_+\widehat{\theta })\delta \widehat{\rho }_+\delta \widehat{\rho }_+(\widehat{\theta }_+\widehat{\theta })\delta \widehat{\rho }=\delta \widehat{\rho }[\widehat{\theta }_+,\delta \widehat{\rho }_+][\delta \widehat{\rho },\widehat{\theta }]\delta \widehat{\rho }_+$$
(207)
which is a sum of two commutators, easy to evaluate. This leads to

$$H_{\mathrm{kin}}^{(2,1)}=\frac{\mathrm{}^2}{8m}\underset{𝐫}{}\frac{\delta \widehat{\rho }}{\rho _0}\left(\rho _0^{1/2}\mathrm{\Delta }\rho _0^{1/2}\rho _0^{1/2}\mathrm{\Delta }\rho _0^{1/2}\right).$$
(208)
To calculate $`H_{\mathrm{kin}}^{(1,2)}`$ we first evaluate

$$H_{\mathrm{kin}}^{(+\mathrm{},2)}=\frac{\mathrm{}^2}{4ml^2}\underset{𝐫,j}{}l^D\left[\sqrt{\widehat{\rho }}(\widehat{\theta }_{+j}\widehat{\theta })^2\sqrt{\widehat{\rho }_{+j}}+\text{h.c.}\right]$$
(209)
and we expand to first order in $`\delta \widehat{\rho }`$, which leads to a sum of terms that are not individually hermitian. We then use the commutation relation (20) to produce hermitian terms, e.g.

$$\delta \widehat{\rho }(\widehat{\theta }_{+j}\widehat{\theta })^2=(\widehat{\theta }_{+j}\widehat{\theta })\delta \widehat{\rho }(\widehat{\theta }_{+j}\widehat{\theta })\frac{i}{l^D}(\widehat{\theta }_{+j}\widehat{\theta }).$$
(210)
The last term of the right hand side of this expression is antihermitian and does not contribute to the final result

$$H_{\mathrm{kin}}^{(1,2)}=\frac{\mathrm{}^2}{4ml^2}\underset{𝐫,j}{}l^D(\widehat{\theta }_{+j}\widehat{\theta })\left(\frac{\rho _{0,+j}^{1/2}}{\rho _0^{1/2}}\delta \widehat{\rho }+\frac{\rho _0^{1/2}}{\rho _{0,+j}^{1/2}}\delta \widehat{\rho }_{+j}\right)(\widehat{\theta }_{+j}\widehat{\theta }).$$
(211)
Finally $`H_{\mathrm{kin}}^{(0,1)}`$ and $`H_{\mathrm{kin}}^{(0,3)}`$ vanish as the odd order expansion of $`\mathrm{exp}[i(\widehat{\theta }_{+j}\widehat{\theta })]`$ is antihermitian.

## Appendix B Corrections to the equations of motion due to $`H_3`$

The Hamiltonian $`H_3`$ gives rise to quadratic corrections to the equations of motion for $`\delta \widehat{\rho }`$ and $`\widehat{\theta }`$. In this appendix, these corrections are calculated explicitly and the thermal average is taken over the equations of motion with the Hamiltonian $`H_2+H_3`$ for the linear part and the Hamiltonian $`H_2`$ for the quadratic corrections. This allows to calculate the first correction to the mean density due to $`H_3`$.

The corrections to the equation of motion for the density fluctuations are given by:

$$\mathrm{}_t\delta \widehat{\rho }|_{H_3}=\frac{\mathrm{}^2}{4ml^2}\underset{j}{}[\{\widehat{\theta }\widehat{\theta }_{+j},\left(\frac{\rho _{0,+j}}{\rho _0}\right)^{1/2}\delta \widehat{\rho }+\left(\frac{\rho _0}{\rho _{0,+j}}\right)^{1/2}\delta \widehat{\rho }_{+j}\}+(+jj)]$$
(212)
where $`\{A,B\}`$ stands for the anticommutator $`AB+BA`$ of two operators. When we take the average with the Hamiltonian $`H_2`$, we use the explicit modal expansion of $`\delta \widehat{\rho }`$ and $`\widehat{\theta }`$ given by Eq.(56). The operator $`\widehat{Q}`$ disappears since Eq. (212) involves only differences of $`\widehat{\theta }`$. Terms with $`\widehat{P}`$ also disappear since $`\widehat{P}_2=0`$. The expectation value of the product $`\widehat{\theta }(𝐫)\delta \widehat{\rho }(𝐫^{})`$, where $`\widehat{\theta }`$ is written without the $`\widehat{Q}`$ and $`\delta \widehat{\rho }`$ is written without the $`\widehat{P}`$, is actually purely imaginary: as $`u_s`$ and $`v_s`$ can be chosen to be real, $`\theta _s^{}=\theta _s`$, see Eq.(II.3). Since $`_t\delta \widehat{\rho }`$ is real, all imaginary contributions to it have to cancel so that the correction to the motion of $`\delta \widehat{\rho }`$ due to $`H_3`$ finally vanish when we take the thermal average:

$$\mathrm{}_t\delta \widehat{\rho }|_{H_3}=0$$
(213)
The corrections to the equation of motion for $`\widehat{\theta }`$ are more involved:

$`\mathrm{}_t\widehat{\theta }|_{H_3}={\displaystyle \frac{1}{2\sqrt{\rho _0}}}({\displaystyle \frac{\mathrm{}^2}{4m\rho _0}}\delta \widehat{\rho }\mathrm{\Delta }\left({\displaystyle \frac{\delta \widehat{\rho }}{\sqrt{\rho _0}}}\right)+{\displaystyle \frac{3\mathrm{}^2}{8m\rho _0^2}}\delta \widehat{\rho }^2\mathrm{\Delta }\left(\sqrt{\rho _0}\right){\displaystyle \frac{\mathrm{}^2}{8m}}\mathrm{\Delta }\left({\displaystyle \frac{\delta \widehat{\rho }^2}{\rho _0^{3/2}}}\right)`$
$`{\displaystyle \frac{\mathrm{}^2}{2ml^2}}{\displaystyle \underset{j}{}}(\sqrt{\rho _{0,+j}}(\widehat{\theta }_{+j}\widehat{\theta })^2+\sqrt{\rho _{0,j}}(\widehat{\theta }_j\widehat{\theta })^2)+{\displaystyle \frac{g_0\sqrt{\rho _0}}{l^D}}{\displaystyle \frac{\mathrm{}^2}{4ml^D\sqrt{\rho _0}}}(\rho _0^{1/2}\mathrm{\Delta }\rho _0^{1/2}\rho _0^{1/2}\mathrm{\Delta }\rho _0^{1/2}))`$ (214)
Fortunately we can use the linear equations of motion (46,47) to significantly simplify the above equation of motion. We rewrite the first term in parentheses of Eq.(B) as:

$$\frac{\mathrm{}^2}{4m\rho _0}\delta \widehat{\rho }\mathrm{\Delta }\left(\frac{\delta \widehat{\rho }}{\sqrt{\rho _0}}\right)=\frac{\delta \widehat{\rho }}{2\rho _0}\left[(U\mu +3g_0\rho _0)\left(\frac{\delta \widehat{\rho }}{\sqrt{\rho _0}}\right)+2\sqrt{\rho _0}\mathrm{}_t\widehat{\theta }\right]$$
(215)
The second term in parentheses of Eq.(B) gives, as $`\sqrt{\rho _0}`$ solves the Gross-Pitaevskii equation:

$$\frac{3\mathrm{}^2}{8m\rho _0^2}\delta \widehat{\rho }^2\mathrm{\Delta }\left(\sqrt{\rho _0}\right)=\frac{3\delta \widehat{\rho }^2}{4\rho _0^{3/2}}(U\mu +g_0\rho _0).$$
(216)
The sum of Eqs.(215),(216) and the third term in parentheses of Eq.(B) lead to:

$$\frac{1}{4}\left(\frac{\mathrm{}^2\mathrm{\Delta }}{2m}+U\mu +g_0\rho _0\right)\left(\frac{\delta \widehat{\rho }^2}{\rho _0^{3/2}}\right)g_0\frac{\delta \widehat{\rho }^2}{\sqrt{\rho _0}}\mathrm{}(_t\widehat{\theta })\frac{\delta \widehat{\rho }}{\sqrt{\rho _0}}.$$
(217)
To rewrite the fourth term in parentheses of Eq.(B), it is convenient to use the following identity:

$`\sqrt{\rho _{0,+j}}(\widehat{\theta }_{+j}\widehat{\theta })^2+\sqrt{\rho _{0,j}}(\widehat{\theta }_j\widehat{\theta })^2+2\widehat{\theta }\left(\sqrt{\rho _{0,+j}}(\widehat{\theta }_{+j}\widehat{\theta })+\sqrt{\rho _{0,j}}(\widehat{\theta }_j\widehat{\theta })\right)`$ $`=`$ $`\sqrt{\rho _{0,+j}}(\widehat{\theta }_{+j}^2\widehat{\theta }^2)+\sqrt{\rho _{0,j}}(\widehat{\theta }_j^2\widehat{\theta }^2)`$ (218)
$`=`$ $`l^2\left(\mathrm{\Delta }\left(\sqrt{\rho _0}\widehat{\theta }^2\right)\widehat{\theta }^2\mathrm{\Delta }\left(\sqrt{\rho _0}\right)\right)`$
leading to:

$`\sqrt{\rho _{0,+j}}(\widehat{\theta }_{+j}\widehat{\theta })^2+\sqrt{\rho _{0,j}}(\widehat{\theta }_j\widehat{\theta })^2=l^2\left(\widehat{\theta }^2\mathrm{\Delta }\sqrt{\rho _0}2\widehat{\theta }\mathrm{\Delta }\left(\sqrt{\rho _0}\widehat{\theta }\right)+\mathrm{\Delta }\left(\sqrt{\rho _0}\widehat{\theta }^2\right)\right).`$ (219)
Using this equality, the Gross-Pitaevskii equation (41) and the equation of motion (47), the fourth term in parentheses of Eq.(B) can be written as:

$$\frac{\mathrm{}^2}{2ml^2}\underset{j}{}\left(\sqrt{\rho _{0,+j}}(\widehat{\theta }_{+j}\widehat{\theta })^2+\sqrt{\rho _{0,j}}(\widehat{\theta }_j\widehat{\theta })^2\right)=\left(\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U\mu +g_0\rho _0\right)\left(\sqrt{\rho _0}\widehat{\theta }^2\right)\mathrm{}\widehat{\theta }\frac{_t\delta \widehat{\rho }}{\sqrt{\rho _0}}.$$
(220)
The sixth (and last) term in parentheses of Eq.(B) can also be transformed using the Gross-Pitaevskii equation (41):

$$\frac{\mathrm{}^2}{4ml^D\sqrt{\rho _0}}\left(\rho _0^{1/2}\mathrm{\Delta }\rho _0^{1/2}\rho _0^{1/2}\mathrm{\Delta }\rho _0^{1/2}\right)=\left(\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U\mu +g_0\rho _0\right)\left(\frac{1}{2l^D\sqrt{\rho _0}}\right).$$
(221)
This leads finally to a rewriting of the thermal average of Eq. (B) as:

$$2\sqrt{\rho _0}\mathrm{}_t\delta \widehat{\theta }|_{H_3}=\left(\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U\mu +g_0\rho _0\right)\left(\frac{\delta \widehat{\rho }^2_2}{4\rho _0^{3/2}}+\sqrt{\rho _0}\widehat{\theta }^2_2\frac{1}{2l^D\sqrt{\rho _0}}\right)g_0\frac{\delta \widehat{\rho }^2_2\rho _0/l^D}{\sqrt{\rho _0}}\frac{\mathrm{}_t\widehat{\theta }\delta \widehat{\rho }_2}{\sqrt{\rho _0}}.$$
(222)
The last term of this expression can be calculated using Eq.(56). The harmonic modes do not contribute since the expectation value of products of $`\widehat{b_s}`$ and $`\widehat{b_s}^{}`$ is time independent. We are left with:

$$_t\widehat{Q}\widehat{P}_2=_t\widehat{Q(0)}\widehat{P}+t\frac{\mu _0^{}}{\mathrm{}}\widehat{P}^2_2=\frac{\mu _0^{}}{\mathrm{}}\widehat{P}^2_2$$
(223)
which gives:

$$\frac{\mathrm{}_t\widehat{\theta }\delta \widehat{\rho }_2}{\sqrt{\rho _0}}=2\mu _0^{}\widehat{P}^2_2_{N_0}\sqrt{\rho _0}.$$
(224)
As a conclusion, the quadratic correction to the first equation of motion can be written as in (60) if one uses the identities:

$`\widehat{B}^{}\widehat{B}`$ $`=`$ $`{\displaystyle \frac{\delta \widehat{\rho }^2}{4\rho _0}}+\rho _0\widehat{\theta }^2{\displaystyle \frac{1}{2l^D}}`$
$`{\displaystyle \frac{\delta \widehat{\rho }^2}{\sqrt{\rho _0}}}`$ $`=`$ $`\sqrt{\rho _0}\left(\widehat{B}+\widehat{B}^{}\right)^2=\sqrt{\rho _0}\left(2\widehat{B}^{}\widehat{B}+\widehat{B}^2+\widehat{B}^2+{\displaystyle \frac{1}{l^D}}\right).`$ (225)
## Appendix C The mean value of $`_t\widehat{\theta }`$ vanishes at equilibrium

As the field degree of freedom $`\widehat{Q}`$, that is the global phase of the field, is not subject to a restoring force in $`H_2`$, it is not totally obvious that the perturbation $`H_3`$ cannot set it into permanent motion. We therefore check this point explicitly here.

The first step is to calculate the mean value of $`\widehat{P}`$ to first order in $`H_3`$. We approximate the unnormalized density operator of the gas at thermal equilibrium to first order in $`H_3`$ using perturbation theory:

$$\sigma =e^{\beta (H_2+H_3)}=e^{\beta H_2}_0^\beta 𝑑\tau e^{(\beta \tau )H_2}H_3e^{\tau H_2}+\mathrm{}.$$
(226)
$`\widehat{P}`$ commutes with $`H_2`$ and has a vanishing mean value in the thermal state corresponding to $`H_2`$ so that, to first order in $`H_3`$:

$$\widehat{P}_3=\beta \widehat{P}H_3_2.$$
(227)
The Hamiltonian $`H_3`$ is a polynomial of degree 3 in $`\widehat{P}`$:

$$H_3=A_0+A_1\widehat{P}+A_2\widehat{P}^2+A_3\widehat{P}^3$$
(228)
where the $`A_i`$ are still operators with respect to the harmonic oscillator variables $`b_s`$. This leads to

$$\widehat{P}_3=\beta \left[A_1_2\widehat{P}^2_2+A_3_2\widehat{P}^4_2\right].$$
(229)
From Wick’s theorem, $`\widehat{P}^4_2=3\widehat{P}^2_2^2`$.

In a second step we calculate $`d\widehat{Q}/dt`$ to first order in $`H_3`$:

$`d\widehat{Q}/dt`$ $``$ $`_{\widehat{P}}(H_2+H_3)_3`$ (230)
$``$ $`\mu _0^{}\widehat{P}_3+A_1_2+3A_3_2\widehat{P}^2_2`$
where the terms coming from $`_{\widehat{P}}H_3`$ are calculated in the thermal state for $`H_2`$ since they are already first order in the perturbation. From the value of $`\widehat{P}_3`$ obtained from (229) and from (61) we obtain the desired result:

$$d\widehat{Q}/dt_3=0$$
(231)
to first order in $`H_3`$.

## Appendix D An equation for $`\{\widehat{\alpha },\widehat{\mathrm{\Lambda }}\}`$

In this appendix, we derive the partial differential equation (91). We first note that $`\widehat{B}_n`$, being a sum of eigenmodes of the operator $`_{\mathrm{GP}}`$, obeys the differential equation for the evolution governed by $`H_2`$:

$$i\mathrm{}_t\left(\begin{array}{c}\widehat{B}_n\\ \widehat{B}_n^{}\end{array}\right)=_{\mathrm{GP}}\left(\begin{array}{c}\widehat{B}_n\\ \widehat{B}_n^{}.\end{array}\right)$$
(232)
We project this equation orthogonally to $`\varphi _0`$ and along $`\varphi _0`$, so that we get the quantum analog of equations (E9), (E10) of cartago , with the simplification that $`\varphi _0(𝐫)`$ is real:

$`i\mathrm{}_t\left(\begin{array}{c}\widehat{\mathrm{\Lambda }}\\ \widehat{\mathrm{\Lambda }}^{}\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}𝒬& 0\\ 0& 𝒬\end{array}\right)_{\mathrm{GP}}\left(\begin{array}{c}\widehat{\mathrm{\Lambda }}\\ \widehat{\mathrm{\Lambda }}^{}\end{array}\right)+(\widehat{\alpha }+\widehat{\alpha }^{})\left(\begin{array}{c}𝒬g_0\rho _0\varphi _0\\ 𝒬g_0\rho _0\varphi _0\end{array}\right)`$ (233)
$`i\mathrm{}{\displaystyle \frac{d\widehat{\alpha }}{dt}}`$ $`=`$ $`l^D{\displaystyle \underset{𝐫}{}}g_0\rho _0\varphi _0\left(\widehat{B}_n+\widehat{B}_n^{}\right)=l^D{\displaystyle \underset{𝐫}{}}g_0\rho _0\varphi _0\left(\widehat{\mathrm{\Lambda }}+\widehat{\mathrm{\Lambda }}^{}\right).`$ (234)
We have introduced the projection matrix

$$𝐫|𝒬|𝐫^{}=\delta _{𝐫,𝐫^{}}l^D\varphi _0(𝐫)\varphi _0(𝐫^{}).$$
(235)
As $`\widehat{\alpha }`$ is antihermitian, the source term in (233) vanishes and one can replace $`\widehat{B}_n`$ by $`\widehat{\mathrm{\Lambda }}`$ in (234).

We use these two equations of motion to calculate the first order time derivative of $`A(𝐫)\{\widehat{\alpha },\widehat{\mathrm{\Lambda }}(𝐫)\}_2`$. We do not give the intermediate result. As $`A`$ is real here, we have the property

$$\{\widehat{\alpha },\widehat{\mathrm{\Lambda }}^{}(𝐫)\}_2=\{\widehat{\alpha },\widehat{\mathrm{\Lambda }}(𝐫)\}_2.$$
(236)
As $`\widehat{\mathrm{\Lambda }}`$ is orthogonal to $`\varphi _0`$ one has

$$𝒬A=A.$$
(237)
All this leads to (91).

## Appendix E Interpretation of $`\chi `$ in the number conserving Bogoliubov approach

We assume here that the gas is a quasi-pure condensate so that $`\varphi _0`$ is now the condensate wavefunction in the Gross-Pitaevskii approximation. We then show that $`\chi (𝐫)/N_0`$, where $`\chi `$ is defined in (94), essentially coincides with the lowest order deviation of the exact condensate wavefunction from the Gross-Pitaevskii prediction $`\varphi _0`$. This deviation was calculated in yvan2 .

We split $`\chi `$ in a component orthogonal to $`\varphi _0`$ and a component colinear to $`\varphi _0`$:

$$\chi (𝐫)=\gamma \varphi _0(𝐫)+\chi _{}(𝐫).$$
(238)
The component $`\gamma `$ has a simple physical interpretation: we sum (95) over $`𝐫`$ after multiplication by $`l^D`$. If we omit the grand canonical term (absent in the canonical treatment of yvan2 ) we obtain

$$N=N_0+2\gamma +\delta N$$
(239)
where

$$\delta Nl^D\underset{𝐫}{}\widehat{\mathrm{\Lambda }}^{}\widehat{\mathrm{\Lambda }}_2$$
(240)
exactly coincides with the mean number of non-condensed particles predicted in yvan2 . The physical interpretation of $`2\gamma `$ is then simple:

$$\delta N_02\gamma $$
(241)
is the correction to apply to the pure condensate prediction for the number of condensate particles in order to recover the correct Bogoliubov prediction. Applying to (94) the matrix $`𝒬`$ (235) projecting orthogonally to $`\varphi _0`$ we obtain

$$\left[\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+U+g_0\rho _0\mu \right]\chi _{}+2𝒬g_0\rho _0\chi _{}+𝒬\left(2g_0\rho _0\gamma \varphi _0+\frac{1}{2}S\right)=0.$$
(242)
We modify slightly the form of the source term $`S`$, eliminating the anticommutator:

$$\{\widehat{\mathrm{\Lambda }}^{}(𝐫^{}),\widehat{\mathrm{\Lambda }}(𝐫)\}=2\widehat{\mathrm{\Lambda }}^{}(𝐫^{})\widehat{\mathrm{\Lambda }}(𝐫)+\frac{1}{l^D}𝐫|𝒬|𝐫^{}.$$
(243)
This leads to the system

$$\left(\begin{array}{cc}𝒬& 0\\ 0& 𝒬\end{array}\right)_{\mathrm{GP}}\left(\begin{array}{c}\chi _{}\\ \chi _{}\end{array}\right)+\left(\begin{array}{c}𝒬S_{\mathrm{eff}}\\ 𝒬S_{\mathrm{eff}}\end{array}\right)=0$$
(244)
with the effective source term

$$S_{\mathrm{eff}}(𝐫)=g_0\rho _0(𝐫)\varphi _0(𝐫)(\delta N_01)+g_0N_0\varphi _0(𝐫)\left(2\widehat{\mathrm{\Lambda }}^{}(𝐫)\widehat{\mathrm{\Lambda }}(𝐫)_2+\widehat{\mathrm{\Lambda }}^2(𝐫)_2\right)l^D\underset{𝐫^{}}{}g_0\rho _0(𝐫^{})\varphi _0(𝐫^{})\left(\widehat{\mathrm{\Lambda }}(𝐫^{})+\widehat{\mathrm{\Lambda }}^{}(𝐫^{})\right)\widehat{\mathrm{\Lambda }}(𝐫)_2$$
(245)
where we used the fact that here $`\widehat{\mathrm{\Lambda }}^2_2=\widehat{\mathrm{\Lambda }}^2_2`$ since the condensate wavefunction is real. Equation (244) is the steady version of Eq.(95) of yvan2 , which gives $`N`$ times the correction to the condensate wavefunction, and the source term (245) exactly coincides with the one of Eq.(96) of yvan2  if one realizes that $`N=N_0`$, so that $`\delta N_0=\delta N`$, in the systematic expansion used in yvan2 .

## Appendix F Corrections to $`g_1`$ due to the cubic Hamiltonian

We calculate the corrections to the first order correlation function due to $`H_3`$ using the perturbative formula (143). A first remark is that

$$H_3(\tau )e^{\tau H_2}H_3e^{\tau H_2}$$
(246)
is still cubic in the operators $`\widehat{b}_s`$ since one has

$`e^{\tau H_2}\widehat{b}_se^{\tau H_2}`$ $`=`$ $`e^{\tau ϵ_s}\widehat{b}_s`$ (247)
$`e^{\tau H_2}\widehat{b}_s^{}e^{\tau H_2}`$ $`=`$ $`e^{\tau ϵ_s}\widehat{b}_s^{}`$ (248)
where $`ϵ_s`$ is the energy of the Bogoliubov mode $`s`$. The second step is to use Wick’s theorem to calculate the expectation values in the thermal state corresponding to Hamiltonian $`H_2`$. One can derive the general formulas

$`A_1A_2A_3e^{i\mathrm{\Delta }\theta }_2`$ $`=`$ $`\left[A_1A_2A_3i\mathrm{\Delta }\theta _2+A_1i\mathrm{\Delta }\theta _2A_2i\mathrm{\Delta }\theta _2A_3i\mathrm{\Delta }\theta _2\right]e^{(\mathrm{\Delta }\theta )^2_2/2}`$ (249)
$`A_1A_2A_3A_4e^{i\mathrm{\Delta }\theta }_2`$ $`=`$ $`\{A_1A_2A_3A_4_2[1+{\displaystyle \frac{1}{2}}(\mathrm{\Delta }\theta )^2_2]{\displaystyle \frac{1}{2}}A_1A_2A_3A_4(\mathrm{\Delta }\theta )^2_2`$ (250)
$`+`$ $`A_1i\mathrm{\Delta }\theta _2A_2i\mathrm{\Delta }\theta _2A_3i\mathrm{\Delta }\theta _2A_4i\mathrm{\Delta }\theta _2\}e^{(\mathrm{\Delta }\theta )^2_2/2}`$
where the $`A_i`$ are linear in $`\delta \widehat{\rho }`$ and $`\widehat{\theta }`$ and have a vanishing mean value. A last point is to realize that some of the obtained terms contain a larger number of factors equal to $`\mathrm{\Delta }\theta `$ than other ones. Since $`\mathrm{\Delta }\theta `$ scales as $`1/\sqrt{\rho _0l^D}`$, see e.g. the expression of $`\widehat{\theta }`$ in terms of the mode functions $`u_s`$, $`v_s`$ in (II.3), the terms with an excess of $`\mathrm{\Delta }\theta `$ factors are higher order in the expansion and are therefore negligible. Note that strictly speaking this argument is correct provided that each factor $`A_i\mathrm{\Delta }\theta _2`$ remains bounded whatever the distance from 0 to $`𝐫`$. This can be checked indeed to be the case: from the form of $`H_3`$ one sees that $`A_i`$ is either $`\delta \widehat{\rho }`$ or the phase difference between two neighboring points of the grid. One can therefore use the approximate identities:

$`A_1A_2A_3e^{i\mathrm{\Delta }\theta }_2`$ $``$ $`A_1A_2A_3i\mathrm{\Delta }\theta _2e^{(\mathrm{\Delta }\theta )^2_2/2}`$ (251)
$`A_1A_2A_3A_4e^{i\mathrm{\Delta }\theta }_2`$ $``$ $`A_1A_2A_3A_4_2e^{(\mathrm{\Delta }\theta )^2_2/2}.`$ (252)
This immediately leads to the identities (144,145).