text:
\section{Introduction and Analysis of Problem}

The appearance of Bose-Einstein condensate in a Bose system is
usually associated with the breaking of the global gauge symmetry,
which is commonly realized by means of the Bogolubov shift in the
field operators [1--4]. The idea of symmetry breaking is in line with
the general understanding that phase transitions of different nature
are accompanied by some symmetry changes. In theoretical description,
there exist several ways of symmetry breaking (see review article [5]).
The most known is the Bogolubov method of quasiaverages, when the
symmetry of the Hamiltonian is disturbed by infinitesimal sources [4].
For breaking the global gauge symmetry, the latter technique is not
always convenient, while the Bogolubov shift is a simple and sufficient
condition for the symmetry breaking [6].

When the global gauge symmetry is broken, then the description of Bose
systems, based on the standard grand canonical ensemble, encounters a
dilemma, first discussed by Hohenberg and Martin [7], who emphasized that
this description results in either nonconserving theories or yields an
unphysical gap in the spectrum of particles. In nonconserving theories,
local conservation laws are not valid, which, at the same time, is
connected with inconsistencies in thermodynamics.

The Hohenberg-Martin dilemma of conserving versus gapless theories is
usually illustrated by considering some approximate schemes [7], for
instance, a kind of mean-field approximations. The origin of this
dilemma can be explained as follows.

Bose-Einstein condensation implies the macroscopic occupation of
the ground-state energy level, when the number of condensed particles
$N_0$ is such that the condensate fraction $N_0/N$ is nonzero in
the thermodynamic limit. The total number of particles in the system,
$N=N_0+N_1$, becomes a sum of the condensed-particle number $N_0$ and
the number $N_1$ of uncondensed particles. According to Bogolubov [1--4],
the number of condensed particles $N_0$ is such that it provides
thermodynamic stability for the system, minimizing the thermodynamic
potential. At the same time, $N_0=N-N_1$, where $N_1=N_1(\mu,T,\rho)$
is a function of the chemical potential $\mu$, temperature $T$, and
density $\rho\equiv N/V$, with $V$ being the system volume. Hence,
$N_0=N_0(\mu,T,\rho)$ is a function of the same variables. Conversely,
the chemical potential $\mu=\mu(T,\rho)$ is a function of $T$ and $\rho$.
Thus, the condition of the stability, realized as the minimization of
the thermodynamic potential, defines the chemical potential $\mu$. It
is important to emphasize that this procedure of minimization does not
depend on approximations, but is generally valid, as was strictly proved
by Ginibre [8].

From another side, the particle spectrum, under the broken gauge symmetry,
has to be gapless. This is, actually, the condition for the existence of
stable Bose-Einstein condensate. Since, if there would be a gap in the
spectrum, there could be no macroscopic occupation of a single ground-state
level. Hugenholtz and Pines [9] showed that the chemical potential
$\mu=\Sigma_{11}(0,0)-\Sigma_{12}(0,0)$ is expressed through the self-energies
$\Sigma_{\alpha\beta}({\bf k},\omega)$, and this relation makes the spectrum gapless for
any Bose system. As far as the self-energies $\Sigma_{\alpha\beta}(0,0)$ are the
functions of $T$ and $\rho$, the Hugenholtz-Pines relation defines the
chemical potential $\mu=\mu(T,\rho)$.

In this way, the condition of thermodynamic stability, in the frame of the
Bogolubov-Ginibre minimization procedure, defines a chemical potential that
we shall denote as $\mu_{BG}$. At the same time, the existence of a gapless
spectrum, connected to the Hugenholtz-Pines relation, also defines a chemical
potential, which may be denoted as $\mu_{HP}$. These chemical potentials,
found in two different ways do not necessarily coincide, and, generally,
they are different: $\mu_{BG}\neq\mu_{HP}$. Therefore, if one accepts as the
system chemical potential the Bogolubov-Ginibre value $\mu_{BG}$, providing
thermodynamic stability, then one acquires an unphysical gap in the spectrum,
proportional to the difference $|\mu_{HP}-\mu_{BG}|^{1/2}$. Conversely,
accepting as the chemical potential the Hugenholtz-Pines form $\mu_{HP}$,
one gets a nonconserving theory with incorrect thermodynamics, since the
stability condition does not hold. This is the origin of the Hohenberg-Martin
dilemma of conserving versus gapless theories [7].

In any case, the fact that $\mu_{BG}\neq\mu_{HP}$ makes the system unstable
and its description not self-consistent. This problem does not arise only
in the limit of asymptotically weak interactions, when the Bogolubov
weakly-nonideal gas approximation is applicable [1,2]. In this limit,
$\mu_{BG}$ and $\mu_{HP}$ asymptotically coincide. However, in any higher
approximation one has $\mu_{BG}\neq\mu_{HP}$.

For practical applications, one usually does the following. When one
is interested solely in the system dynamics, but not in its spectrum,
one derives the equations from a variational principle, which
guarantees the validity of conservation laws [10,11]. For an
equilibrium system, this is equivalent to the Bogolubov-Ginibre
variational procedure. And, when one studies only the system spectrum,
one accepts the Hugenholtz-Pines relation, forgetting about inconsistent
thermodynamics and instability. Clearly, such palliative ways are not
satisfactory. The principal problem remains how to make the general
theory both conserving as well as gapless.

There have been several attempts to cure the problem, which could be
classified into three groups:

The most often used trick is the omission of anomalous averages. As is
clear, this is a rather unjustified way, since, as soon as the global
gauge symmetry is broken, the anomalous averages do exist and are not
zero. One, sometimes, ascribes this trick to Popov, calling it the "Popov
approximation". However, it is sufficient to look attentively at the original
works by Popov [12--14], which are usually cited in this respect, in order
to realize that he has never suggested anything like that. He considered the
properties of a Bose gas in the vicinity of the critical temperature $T_c$,
honestly calculating all terms, normal and anomalous. When temperature tends
to $T_c$, the condensate density tends to zero together with the anomalous
averages. Actually, both these quantities, the condensate density and the
anomalous averages are the order parameters, appearing together in the
broken-symmetry phase, and disappearing also together, when the symmetry is
getting restored at the critical temperature. Contrary to this, the normal
average, that is, the density of uncondensed particles, increases when
approaching the critical temperature, reaching at $T_c$ the total system
density. This is why at the close vicinity of $T_c$, the anomalous averages
become, without any special assumptions, smaller than the normal ones. However
at low temperatures $T\ll T_c$, the anomalous averages not merely can be of the
order of the normal averages, but can even be much larger than the latter, as
direct calculations show [15]. Moreover, omitting the anomalous averages makes
the system principally unstable [15,16]. Thus, this trick neither has anything
to do with Popov nor can be accepted as a reasonable approximation at low
temperatures. Additionally, if one would wish to ascribe a name to this trick,
one should know that the first person, who really suggested and used it, was
Shohno [17]. And this was known in literature as the Shohno model. The word
"model" is appropriate here, since this is, actually, just a model, but not
a justified approximation. For example, Reatto and Straley [18] used the term
"Shohno model" and studied its properties.

Another way of removing the gap in the spectrum is to calculate the chemical
potential and self energy in different approximations. Then one defines the
chemical potential from the thermodynamic stability condition in one
approximation, but for calculating the self-energy, one invokes a higher-order
approximation, so that to cancel the gap. The additional higher-order terms
can be motivated by Bethe-Salpeter or scattering-matrix approximations [19,20].
This way is what was called by Bogolubov [4] the mismatch of approximations.
Bogolubov already mentioned [4] that such a mismatch can really influence
the appearance or disappearance of the spectrum gap, but renders the theory
not self-consistent and the system unstable.

One can also kill the spectrum gap by adding to the self-energy phenomenological
terms, such that to cancel the gap [21--23]. This way is, clearly, equivalent to
the previous one, since changing the self-energy by invoking some higher-order
approximations is not unique and, hence, is also phenomenological.

As is evident, all these ways, attempting to cure the problem, are based on
a kind of the mismatch of approximations and, as a result, they have the same
common defects:

(i) They are not uniquely defined since there exists an infinite number of
particular tricks for removing the spectrum gap. So, the ambiguity remains
[23].

(ii) They are not self-consistent, involving in this or that way the mismatch
of approximations [4].

(iii) They, as a rule, correspond to an unstable system, either with a not
minimal thermodynamic potential or with a divergent susceptibility [15,16].

(iv) The order of the condensation transitions changes, resulting in a
first-order phase transition, instead of the correct second-order one. This
happens because of the internal inconsistencies in the description. The
disruption of the phase-transition order is a common feature of inconsistent
approximations, which either do not satisfy the stability conditions or
possess a spectrum gap [18,21--23]. This was already noticed in the early
works analysing the Hartree-Fock-Bogolubov approximation having a gap in the
spectrum [18,24--27]. The generality of such a change of the phase-transition
order from second to first in different mean-field approximations was
emphasized in a detailed discussion by Baym and Grinstein [28] and recently
by Kita [23]. As is clear, the thermodynamics of a system with a first-order
phase transition is rather different from that of a system displaying a
second-order transition. And in the vicinity of the critical point, the
thermodynamics in these two cases differs drastically [23,28].

(v) Finally, the mismatch of approximations, by its own, is not a regular
procedure. For each given approximation, it is necessary to invent special
tricks, which, as is stressed above, are not uniquely defined, hence, ambiguous.
There is no a general rule how to do this in a unique way for approximations
of different order.

Thus, we have to conclude that the Hohenberg-Martin dilemma [7] remains
unsolved. The methods, based on the mismatch of approximations, are not
self-consistent, having several internal defects discussed above. In order
that an effective theory be self-consistent, all dynamic and thermodynamic
equations must be derived from the same Hamiltonian, or Lagrangian, and
treated in one chosen approximation, without involving the approximation
mismatch [29,30]. The most recent thorough discussion of the Hohenberg-Martin
dilemma, with many citations, can be found in the review article by
Andersen [31].

In the paper [6], the idea was advanced that the problem can be solved by
employing a representative statistical ensemble. In Refs. [32,33], it was
shown how to make the Hartree-Fock-Bogolubov (HFB) approximation for a
dilute gas both conserving and gapless. The aim of the present paper is
to develop a general approach, independent of particular approximations,
for the self-consistent treatment of arbitrary Bose systems with broken
gauge symmetry and to demonstrate on the most general footing that the
resulting theory is really completely self-consistent. The approach is
based on the notion of representative statistical ensembles, whose
general formulation is given in Sec. 2 for both equilibrium as well
as nonequilibrium systems. This notion is specified in Sec. 3 for
Bose systems with broken gauge symmetry. Thermodynamic self-consistency
of the theory is emphasized in Sec. 4. Operator equations of motion are
obtained in Sec. 5 for a Hamiltonian with an arbitrary interaction
potential. Sec. 6 demonstrates that the local conservation laws are
valid on the operator level, hence, being automatically satisfied for
the related average quantities. The equation for the condensate wave
function is derived and analyzed in Sec. 7 for an arbitrary Bose system.
In Sec. 8 a uniform Bose system is considered and illustrated for the HFB
approximation. The behaviour of the condensate and superfluid fractions
is studied in Sec. 9. The equations for Green functions are presented in
Sec. 10, where it is shown, by using the Bogolubov theorem [4], that the
Hugenholtz-Pines relation follows, thus, proving the complete
self-consistency of the developed approach. Sec. 11 is the conclusion.

Throughout the paper, the system of units is used, where the Planck and
Boltzmann constants are set to unity, $\hbar=1$, $k_B=1$.

\section{Representative Statistical Ensembles}

The idea of representative statistical ensembles goes back to Gibbs [34],
who mentioned that to prescribe a canonical distribution for a system
may be not sufficient, but this distribution has to be complimented by
those constraints and conditions that provide a correct representation
of the considered statistical system. The term "representative ensembles"
was used by ter Haar [35,36], who investigated the problem of a proper
representation of equilibrium statistical systems. Equilibrium and
quasiequilibrium representative ensembles were described in the
review article [5]. In this section, we shall formulate the notion
of representative statistical ensembles for both equilibrium and
nonequilibrium systems.

Representative statistical ensembles for equilibrium systems are sometimes
also termed as generalized Gibbs ensembles, subjective or conditional ensembles.
Their mathematical construction is based on the conditional maximization of the
Gibbs entropy, as was done by Janes [37,38]. This concept was also employed
by Girardeau [39].

Let us start with equilibrium systems. An {\it equilibrium statistical
ensemble}, by definition, is a pair $\{{\cal F},\hat\rho\}$ of the space ${\cal F}$
of microstates and a statistical operator $\hat\rho$. In order to correctly
define the latter, it is necessary, according to Gibbs [34], to take into
account all conditions and constraints, imposed on the system. Suppose, we
have a set $\{\hat C_i\}$ of self-adjoint operators, defined on the space
${\cal F}$, which will be called condition operators. This is because these
operators enter the {\it statistical conditions}
\begin{equation}
\label{1}
C_i = \; <\hat C_i>\; = \; {\rm Tr}\hat\rho\; \hat C_i \; ,
\end{equation}
which are necessary to take into account for correctly representing the
considered system. The trace operation in Eq. (1) is over the given space
${\cal F}$, which can be defined as a Fock space. The first evident constraint
is the normalization condition for the statistical operator,
\begin{equation}
\label{2}
1 \; = \; < \hat 1_{\cal F}> \; =  \; {\rm Tr}\hat\rho \; ,
\end{equation}
where $\hat 1_{\cal F}$ is the unity operator in ${\cal F}$. The definition of the
internal energy
\begin{equation}
\label{3}
E\; = \; <\hat H> \; = \; {\rm Tr}\hat\rho\; \hat H \; ,
\end{equation}
as the average of a Hamiltonian $\hat H$, is another standard statistical
condition. But, in addition to constraints (2) and (3), there can exist any
number of other statistical conditions (1). The conditional maximization of
the Gibbs entropy
\begin{equation}
\label{4}
S \equiv - {\rm Tr}\hat\rho \; \ln\hat\rho
\end{equation}
is equivalent to the unconditional minimization of the {\it information
functional} [40], defined as
\begin{equation}
\label{5}
I[\hat\rho] = - S + \lambda_0\left ({\rm Tr}\hat\rho - 1\right ) +
\beta \left ({\rm Tr}\hat\rho\; \hat H - E \right ) + \beta \sum_i \nu_i
\left ({\rm Tr}\hat\rho\; \hat C_i - C_i \right ) \; ,
\end{equation}
in which the standard conditions (2) and (3) are included explicitly.
The quantities $\lambda_0$, $\beta$, and $\beta\nu_i$ are the appropriate
Lagrange multipliers, with $\beta\equiv 1/T$ being inverse temperature.
Minimizing the information functional (5) with respect to the statistical
operator $\hat\rho$ gives
\begin{equation}
\label{6}
\hat\rho = \frac{1}{Z} \; e^{-\beta H} \; ,
\end{equation}
where $Z\equiv\exp(1+\lambda_0)$ is the partition function and the {\it
grand Hamiltonian}
\begin{equation}
\label{7}
H \equiv \hat H + \sum_i \nu_i \hat C_i
\end{equation}
is introduced. The representative statistical ensemble, under
constraints (1), (2), and (3), is then the pair $\{{\cal F},\hat\rho\}$
of a space ${\cal F}$ of microstates and the statistical operator (6),
with the grand Hamiltonian (7). When one of the condition operators
$\hat C_i$ is the number-of-particle operator $\hat N$ and the related
Lagrange multiplier $\nu_i=-\mu$, one gets a particular form of the
grand Hamiltonian $H=\hat H-\mu\hat N$. However, any other necessary
constraints can be included, resulting in the general expression (7)
for the grand Hamiltonian.

The condition operators $\hat C_i$ are to be self-adjoint,
$\hat C_i^+=\hat C_i$, so that the grand Hamiltonian (7) be also
self-adjoint. In many cases, the condition operators are taken as
integrals of motion, such that $[\hat C_i,\hat H]=0$. But this is not
compulsory. For instance, the number-of-particle operator $\hat N$
does not commute with the Hamiltonian energy operator $\hat H$, when
the global gauge symmetry is broken.

The representative ensemble $\{{\cal F},\hat\rho\}$, with the statistical
operator (6) and the grand Hamiltonian (7), define all thermodynamics
of an equilibrium, or stationary, system. The construction of this
ensemble has been more or less straightforward, following the ideas of
Gibbs [34], ter Haar [35,36], and Janes [37,38], as is reviewed in Refs.
[5,40]. But the definition of representative ensembles for arbitrary
nonequilibrium systems is not evident. Below, we give the generalization
of the notion of representative ensembles for nonequilibrium statistical
systems.

To describe a nonequilibrium system, one needs, in addition to the space
of microstates ${\cal F}$ and the initial value of the statistical operator
$\hat\rho=\hat\rho(0)$, to define the temporal evolution of the system.
This evolution can be described by the time-dependent statistical
operators $\hat\rho(t)$, satisfying the Liouville equation, or by the
time dependence of physical operators, satisfying the Heisenberg equation.
Equivalently, the time evolution can be associated with the evolution
operator, satisfying the Schr\"odinger equation. Keeping in mind any of
these ways, we may denote the prescribed temporal evolution by the symbol
$\partial/\partial t$. Then a {\it nonequilibrium statistical ensemble} is a
triplet $\{{\cal F},\hat\rho,\partial/\partial t\}$. Clearly, when the time evolution
is absent, or trivial, this definition reduces to that for the equilibrium
case.

To be more specific, we may remember that each system is characterized by
some dynamical variables, such as field operators. Let us keep in mind a
set of field operators $\psi(x,t)$, whose particular representation is not
important at this stage. For example, the variable $x$ can represent
real-space coordinates or momenta. It may also include other continuous
or discrete variables, such as the spin indices or component labels. All
physical operators, such as the Hamiltonian energy operator $\hat H[\psi]$,
are functionals of the field operators.

The most general way for describing the system dynamics, as is known [41],
is the extremization of the action functional. Implementing this for our
case, we need the Lagrangian
\begin{equation}
\label{8}
\hat L[\psi] \equiv \hat E[\psi] - \hat H[\psi] \; ,
\end{equation}
in which the notation
\begin{equation}
\label{9}
\hat E[\psi] \equiv \int \psi^\dagger(x,t)\; i \;
\frac{\partial}{\partial t}\; \psi(x,t)\; dx
\end{equation}
for the temporal energy operator is used.

To make the ensemble representative, we have to take account of all
additional conditions uniquely characterizing the system. This implies
that, similarly to Eq. (1), we have to take care of the {\it statistical
conditions}
\begin{equation}
\label{10}
C_i \; = \; <\hat C_i[\psi]> \; ,
\end{equation}
where $\hat C_i[\psi]$ are the appropriate condition operators and
$$
<\hat C_i[\psi]> \; \equiv \; {\rm Tr}\hat\rho(0)\;
\hat C_i[\psi(x,t)] = {\rm Tr}\hat\rho(t)\; \hat C_i [\psi(x,0)] \; .
$$

The principle of action extremization, under the given statistical
conditions (10), is equivalent to the unconditional extremization of
the effective action
\begin{equation}
\label{11}
A[\psi] = \int \left\{ \hat L[\psi] -\nu_i \hat C_i[\psi]
\right \}\; dt
\end{equation}
with the Lagrange multipliers $\nu_i$ guaranteeing the validity
of conditions (10). Combining Eqs. (8) and (11), we can rewrite the
effective action (11) as
\begin{equation}
\label{12}
A[\psi] = \int \left\{ \hat E[\psi] - H[\psi]\right \} \; dt \; ,
\end{equation}
with the {\it grand Hamiltonian}
\begin{equation}
\label{13}
H[\psi] \equiv \hat H[\psi] + \sum_i \nu_i \; \hat C_i[\psi] \; ,
\end{equation}
having the same form as in Eq. (7). The extremization of the action
functional with respect to field operators means the variational
equation
\begin{equation}
\label{14}
\frac{\delta A[\psi]}{\delta\psi^\dagger(x,t)} = 0
\end{equation}
plus its Hermitian conjugate. Equation (14), in view of action (12),
is identical to the equation
\begin{equation}
\label{15}
i\; \frac{\partial}{\partial t}\; \psi(x,t) =
\frac{\delta H[\psi]}{\delta\psi^\dagger(x,t)} \; .
\end{equation}

The evolution equation (15) is what one needs for a complete
definition of the nonequilibrium representative ensemble. It is
important to stress that the system dynamics is governed by the same
grand Hamiltonian as its thermodynamics.

We may also note that in the Heisenberg representation, as is well
known (see, e.g., Refs. [40,41]), the variational equation (15) is
the same as the Heisenberg equation
$$
i\; \frac{\partial}{\partial t}\; \psi(x,t) =
\left [ \psi(x,t), \; H[\psi] \right ] \; .
$$
The time evolution of the field operators is given by the form
$$
\psi(x,t) = \hat U^+(t) \psi(x,0) \hat U(t) \; ,
$$
expressed through the evolution operator $\hat U(t)$ satisfying the
Schr\"odinger equation
$$
i\; \frac{d}{dt} \; \hat U(t) =  H[\psi(x,0)] \hat U(t) \; .
$$
Respectively, the time dependence of the statistical operator
$\hat\rho(t)$ stems from the Liouville equation, yielding
$$
\hat\rho(t) = \hat U(t)\; \hat\rho(0)\; \hat U^+(t) \; .
$$
In any case, it is the grand Hamiltonian (13), which governs the temporal
evolution of a nonequilibrium system.

\section{Broken Gauge Symmetry}

Now, let us specify the general notion of representative statistical
ensembles, formulated above, for Bose systems, in which there exists the
critical temperature $T_c$ below which the global $U(1)$ gauge symmetry
becomes broken. For concreteness, we keep in mind a one-component system
of spinless particles, characterized by the field operators satisfying
the Bose commutation relations.

Above the critical temperature $T_c$, the system is described by field
operators $\psi=\psi({\bf r},t)$ and the conjugate $\psi^\dagger$, being
functions of spatial, ${\bf r}$, and temporal, $t$, variables. These
operators generate the Fock space ${\cal F}(\psi)$ on which all physical
operators are defined. The related mathematical details can be found in
the books [40--42]. The Hamiltonian energy operator $\hat H[\psi]$ is a
gauge-invariant functional of the field operators. The number-of-particle
operator $\hat N[\psi]$ is normalized to the total number of particles
$N=<\hat N[\psi]>$. Therefore the grand Hamiltonian is
$$
H[\psi] = \hat H[\psi] - \mu \hat N[\psi] \qquad (T > T_c) \; .
$$
This grand Hamiltonian, entering the statistical operator (6),
characterizes the representative statistical ensemble for the normal Bose
system, above the critical temperature, where the global gauge symmetry is
preserved.

Below the critical temperature $T_c$, the global gauge symmetry becomes
broken. The symmetry breaking is realized by the Bogolubov shift
\begin{equation}
\label{16}
\psi({\bf r},t) \longrightarrow \; \hat \psi({\bf r},t) \; \equiv
\; \eta({\bf r},t) + \psi_1({\bf r},t) \; ,
\end{equation}
where $\eta({\bf r},t)$ is the condensate wave function, while $\psi_1({\bf r},t)$
is the field operator of uncondensed particles. The field operators $\psi_1$
and $\psi_1^\dagger$ generate the Fock space ${\cal F}(\psi_1)$, which all physical
operators are to be defined on. In Eq. (16), the condensate wave function
$\eta({\bf r},t)$, strictly speaking, is assumed to be factored with the unity
operator $\hat 1_{\cal F}$ in ${\cal F}(\psi_1)$. However here and in what follows,
we shall use the common way of omitting the explicit appearance of the unity
operator, in order not to make formulas too cumbersome. The condensate
function $\eta({\bf r},t)$ can also be termed as the coherent field, since the
related coherent state $|\eta>$ is the vacuum state in the space ${\cal F}(\psi_1)$.
This and other mathematical details can be found in Refs. [40,42,43].

Thus, instead of one field variable $\psi({\bf r},t)$ above $T_c$, for the
Bose system below $T_c$, where the gauge symmetry is broken, there arise
two field variables, the condensate function (coherent field) $\eta({\bf r},t)$
and the field operator $\psi_1({\bf r},t)$ of uncondensed particles. These two
dynamical variables are, of course, assumed to be linearly independent,
being orthogonal to each other,
\begin{equation}
\label{17}
\int \eta^*({\bf r},t) \psi_1({\bf r},t)\; d{\bf r} = 0 \; .
\end{equation}

It is of principal importance to emphasize that the spaces ${\cal F}(\psi)$
for $T>T_c$ and ${\cal F}(\psi_1)$ for $T<T_c$ are mutually orthogonal [43,44].
The field operators $\psi$ and $\psi_1$ are defined on different spaces,
${\cal F}(\psi)$ and ${\cal F}(\psi_1)$, respectively, realizing two different
unitary nonequivalent operator representations, with the Bose commutation
relations [43,44]. As soon as the gauge symmetry is broken, one has to deal
with the space ${\cal F}(\psi_1)$. It would be principally incorrect to work,
first, in the space ${\cal F}(\psi)$, accomplishing there some transformations,
and then to pass to the space ${\cal F}(\psi_1)$ by breaking the symmetry with
the Bogolubov shift (16). As is shown in Ref. [6], such a procedure leads to
wrong results. From the mathematical point of view, it is absolutely obvious
that any manipulations must be accomplished in one given space, where all
operators are defined.

Now, in the space ${\cal F}(\psi_1)$, we have two normalization conditions for
two linearly independent field variables, $\eta({\bf r},t)$ and $\psi_1({\bf r},t)$.
The condensate function is normalized to the number of condensed particles
\begin{equation}
\label{18}
N_0 = \int |\eta({\bf r},t)|^2 \; d{\bf r} \; ,
\end{equation}
which is assumed to be macroscopic, such that the condensate fraction
$N_0/N$ be nonzero in the thermodynamic limit. This normalization can be
rewritten in the standard form of the statistical conditions (10) by using
the operator
\begin{equation}
\label{19}
\hat N_0 \equiv N_0 \hat 1_{\cal F} \; ,
\end{equation}
where $\hat 1_{\cal F}$ is the unity operator in ${\cal F}(\psi_1)$. Then Eq. (18)
transforms to the statistical condition
\begin{equation}
\label{20}
N_0 \; =  \; < \hat N_0> \; ,
\end{equation}
in which, as in what follows, the averaging is over the space
${\cal F}(\psi_1)$. The second normalization is, clearly, for the number of
uncondensed particles
\begin{equation}
\label{21}
N_1 \; = \; < \hat N_1> \; ,
\end{equation}
with the corresponding operator
\begin{equation}
\label{22}
\hat N_1 \; \equiv \; \int \psi_1^\dagger({\bf r},t)
\psi_1({\bf r},t)\; d{\bf r} \; .
\end{equation}
The total number of particles
\begin{equation}
\label{23}
N \; = \; <\hat N> \; = \; N_0 + N_1
\end{equation}
is the average of the operator
\begin{equation}
\label{24}
\hat N \; \equiv \; \int \hat\psi^\dagger({\bf r},t)
\hat\psi({\bf r},t)\; d{\bf r} \; = \; \hat N_0 + \hat N_1 \; .
\end{equation}
Normalization (23) follows from Eqs. (20) and (21). Therefore, among three
normalization conditions, (20), (21), and (23), only two can be treated as
independent. Generally, any combination of the pairs, for $N_0$ and $N_1$,
or for $N_0$ and $N$, or for $N_1$ and $N$, could be chosen. For the sake
of symmetry, we prefer to choose the normalization conditions (20) and (21).

When the gauge symmetry is broken, the average $<\psi_1>$ may become nonzero.
This, however, would mean that quantum numbers, as spin or momentum, are not
conserved. In order to avoid this unpleasant situation, one has to impose an
additional constraint
\begin{equation}
\label{25}
<\psi_1({\bf r},t)> \; = \; 0 \; .
\end{equation}
The latter can be reduced to the standard form of the statistical conditions
(10) by defining a self-adjoint operator
\begin{equation}
\label{26}
\hat\Lambda[\psi_1] \; \equiv \;
\int \left [ \lambda({\bf r},t) \psi_1^\dagger({\bf r},t) + \lambda^*({\bf r},t)\psi_1({\bf r},t)
\right ]\; d{\bf r} \; ,
\end{equation}
in which $\lambda({\bf r},t)$ is a complex function. Then constraint (25) can be
rewritten as the {\it quantum-number conservation condition}
\begin{equation}
\label{27}
<\hat\Lambda[\psi_1]> \; = \; 0 \; .
\end{equation}

In this way, for the correct representation of a Bose system with
broken gauge symmetry, we must work in the space ${\cal F}(\psi_1)$ and take
into account three statistical conditions, (20), (21), and (27). The
corresponding representative ensemble is constructed following the general
procedure, formulated in Sec. 2.

For an equilibrium system, according to Eq. (5), we have the information
functional
$$
I[\hat\rho] ={\rm Tr}\hat\rho\; \ln\hat\rho + \lambda_0 \left (
{\rm Tr}\hat\rho -1 \right ) + \beta \left ( {\rm Tr}\hat\rho\;
\hat H[\eta,\psi_1] - E \right )
$$
\begin{equation}
\label{28}
 - \; \beta\mu_0 \left ( {\rm Tr}\hat\rho \; \hat N_0 - N_0 \right ) -
\beta\mu_1 \left ( {\rm Tr}\hat\rho \; \hat N_1 - N_1 \right ) -
\beta {\rm Tr}\hat\rho \; \hat \Lambda[\psi_1] \; .
\end{equation}
Minimizing the latter yields the statistical operator
\begin{equation}
\label{29}
\hat\rho = \frac{1}{Z} \; \exp\left\{ -\beta H[\eta,\psi_1]
\right \} \; ,
\end{equation}
with the partition function
$$
Z \equiv {\rm Tr} \exp\left\{ -\beta H[\eta,\psi_1]
\right \} \; ,
$$
and the grand Hamiltonian
\begin{equation}
\label{30}
H[\eta,\psi_1] \equiv \hat H [\eta,\psi_1] -\mu_0 \hat N_0 -
\mu_1\hat N_1 -\hat\Lambda[\psi_1] \; ,
\end{equation}
in agreement with Eq. (7).

For the general case of a nonequilibrium system, we have the Lagrangian
$$
\hat L [\eta,\psi_1] = \hat E[\eta,\psi_1] - \hat H[\eta,\psi_1] \; ,
$$
in which
$$
\hat E[\eta,\psi_1] =\int \left\{ \eta^*({\bf r},t)\; i \;
\frac{\partial}{\partial t}\; \eta({\bf r},t) + \psi_1^\dagger({\bf r},t)\; i \;
\frac{\partial}{\partial t} \; \psi_1({\bf r},t)\right \} \; d{\bf r} \; .
$$
Then, similarly to Eq. (12), the effective action is
\begin{equation}
\label{31}
A[\eta,\psi_1] = \int \left\{ \hat E[\eta,\psi_1] - H][\eta,\psi_1]
\right \} \; dt \; ,
\end{equation}
with the same grand Hamiltonian (30). Since we have now two linearly
independent field variables, the extremization of action (31) gives
two variational equations, for the condensate function,
\begin{equation}
\label{32}
\frac{\delta A[\eta,\psi_1]}{\delta\eta^*({\bf r},t)} = 0 \; ,
\end{equation}
and for the field operators of uncondensed particles,
\begin{equation}
\label{33}
\frac{\delta A[\eta,\psi_1]}{\delta\psi^\dagger_1({\bf r},t)} = 0 \; .
\end{equation}
Equations (32) and (33), in view of the effective action (31), are
equivalent to the evolution equations
\begin{equation}
\label{34}
i\; \frac{\partial}{\partial t}\; \eta({\bf r},t) =
\frac{\delta H[\eta,\psi_1]}{\delta\eta^*({\bf r},t)}
\end{equation}
and, respectively,
\begin{equation}
\label{35}
i\; \frac{\partial}{\partial t}\; \psi_1({\bf r},t) =
\frac{\delta H[\eta,\psi_1]}{\delta\psi_1^\dagger({\bf r},t)} \; .
\end{equation}
And, as usual, these equations are to be complimented by their
Hermitian conjugate.

Thus, the representative statistical ensemble
$\{{\cal F}(\psi_1),\hat\rho,\partial/\partial t\}$ for a Bose system with broken
global gauge symmetry is defined in complete agreement with the general
theory of Sec. 2. The dynamics and thermodynamics of such a system are
governed by the grand Hamiltonian (30). It is only by accurately taking
into account all imposed constraints that it is possible to correctly
define the representative ensemble and to develop a self-consistent
theory, avoiding any internal defects and paradoxes. The imposed
statistical conditions (20), (21), and (27) lead to the grand Hamiltonian
(30), with the Lagrange multipliers $\mu_0$, $\mu_1$, and $\lambda({\bf r},t)$.
The form of this Hamiltonian is more general than that of the trivial
expression $\hat H-\mu\hat N$. For a system with broken gauge symmetry,
the number of the Lagrange multipliers in the form $\hat H -\mu\hat N$ is
smaller than the number of imposed constraints. As a result, the problem
becomes mathematically overdefined, which leads to the inconsistencies
described in Sec. 1.

\section{Self-Consistent Thermodynamic Relations}

For an equilibrium system, the condensate function does not depend on
time, $\eta({\bf r},t)=\eta({\bf r})$. The grand thermodynamic potential
\begin{equation}
\label{36}
\Omega = - T\ln{\rm Tr}\exp\left ( -\beta H[\eta,\psi_1] \right )
\end{equation}
is defined through the grand Hamiltonian (30). As is evident by this
definition,
$$
\frac{\partial\Omega}{\partial\mu_0} = -N_0 \; , \qquad
\frac{\partial\Omega}{\partial\mu_1} = -N_1 \; .
$$
Varying potential (36) with respect to $\eta({\bf r})$ gives the equation
\begin{equation}
\label{37}
\frac{\delta\Omega}{\delta\eta^*({\bf r})} = \;
< \frac{\delta H[\eta,\psi_1]}{\delta\eta^*({\bf r})} > \; = \; 0
\end{equation}
for the condensate function, in agreement with Eq. (34), when in
equilibrium $\partial\eta({\bf r})/\partial t=0$. In order that the system with
Bose-Einstein condensate be stable, the thermodynamic potential (36)
is to be minimal with respect to the number of condensed particles $N_0$,
so that
\begin{equation}
\label{38}
\frac{\delta\Omega}{\delta N_0} = \;
< \frac{\partial H[\eta,\psi_1]}{\partial N_0} > \; = \; 0 \; ,
\end{equation}
which is the Bogolubov-Ginibre stability condition. Under normalization
(18), the dependence on $N_0$ comes through the condensate function
$\eta({\bf r})$, because of which
$$
\frac{\partial\Omega}{\partial N_0} = \int \left [
\frac{\delta\Omega}{\delta\eta({\bf r})}\; \frac{\partial\eta({\bf r})}{\partial N_0}
+ \frac{\delta\Omega}{\delta\eta^*({\bf r})}\; \frac{\partial\eta^*({\bf r})}{\partial N_0}
\right ]\; d{\bf r} \; .
$$
Therefore Eq. (38) is a direct consequence of Eq. (37). Of course, for
a uniform system, Eqs. (37) and (38) identically coincide.

The average densities of condensed, uncondensed, and all particles are
\begin{equation}
\label{39}
\rho_0 \equiv \frac{N_0}{V} \; , \qquad \rho_1 \equiv \frac{N_1}{V} \; ,
\qquad \rho \equiv \frac{N}{V} \; ,
\end{equation}
respectively. The related fractions of condensed and uncondensed particles
are
\begin{equation}
\label{40}
n_0 \equiv \frac{N_0}{N} \equiv \frac{\rho_0}{\rho} \; , \qquad
n_1 \equiv \frac{N_1}{N} \equiv \frac{\rho_1}{\rho} \; .
\end{equation}
From these definitions, one has
\begin{equation}
\label{41}
\rho = \rho_0 + \rho_1 \; , \qquad n_0 + n_1 =  1 \; .
\end{equation}
What one can usually fix in experiment is temperature $T$ and the total
density $\rho$.

The system free energy  writes as
\begin{equation}
\label{42}
F = \Omega + \mu_0 N_0 + \mu_1 N_1 \; .
\end{equation}
Recall that $N_0 = N_0(T,\rho)$ is such that to guarantee the system
thermodynamic stability. From the equations of motion, together with the
relation $N=N_0+N_1$, one finds $N_1=N_1(T,\rho)$. As is mentioned above,
in experiment, only the total number of particles can be fixed, except
temperature and volume. Hence, the free energy can be written as
\begin{equation}
\label{43}
F = \Omega + \mu N \; ,
\end{equation}
which can be considered as the definition of the system chemical
potential $\mu$. Comparing Eqs. (42) and (43) gives
\begin{equation}
\label{44}
\mu=\mu_0 n_0 +  \mu_1 n_1 \; .
\end{equation}
Then one can define the free energy $F=F(T,V,N)$ as a function of temperature,
volume, and the total number of particles, with the differential
\begin{equation}
\label{45}
dF = -S\; dT - P\; dV + \mu\; dN \; ,
\end{equation}
in which $S$ is entropy and $P$ is pressure. From here, all measurable
thermodynamic quantities are calculated in the standard way. For example,
the system chemical potential is
\begin{equation}
\label{46}
\mu = \left ( \frac{\partial F}{\partial N} \right )_{TV} \; .
\end{equation}
As is evident, Eqs. (44) and (46) are identical by the above definitions.

At zero temperature, the free energy $F=E-TS$ reduces to the internal
energy
\begin{equation}
\label{47}
E \equiv \; < \hat H[\eta,\psi_1]> \; = \; < H[\eta,\psi_1]> + \mu N \; .
\end{equation}
Therefore the chemical potential (46) becomes
$$
\mu = \left ( \frac{\partial E}{\partial N} \right )_V \qquad
(T=0) \; .
$$

The grand potential (36) is a function $\Omega=\Omega(T,V,\mu)$ of temperature,
volume, and chemical potential, with the differential
\begin{equation}
\label{48}
d\Omega = - S\; dT - P\; dV - N\; d\mu \; .
\end{equation}
It is easy to check that the thermodynamic and Gibbs entropies coincide,
so that
\begin{equation}
\label{49}
S = - \; \frac{\partial\Omega}{\partial T} = - {\rm Tr}\hat\rho \;
\ln\hat\rho \; .
\end{equation}
This immediately follows from the direct differentiation of the grand
potential (36) and the form of the statistical operator (29).

It is also possible to prove (see the proof in the Appendix A), that the
standard relation
\begin{equation}
\label{50}
\frac{\partial N}{\partial\mu} = \beta\Delta^2(\hat N)
\end{equation}
holds, in which $\Delta^2(\hat N)$ is the dispersion of the number-of-particle
operator $\hat N=\hat N_0+\hat N_1$. The dispersion of a self-adjoint operator
$\hat A$ is
$$
\Delta^2(\hat A) \equiv \; < \hat A^2> - <\hat A>^2 \; .
$$

In this way, all thermodynamics quantities and relations are defined
self-consistently.

\section{Operator Equations of Motion}

The equations of motion for the field variables $\eta({\bf r},t)$ and
$\psi_1({\bf r},t)$ are given by Eqs. (34) and (35). To specify them, let
us take the Hamiltonian energy operator in the standard form
$$
\hat H[\eta,\psi_1] = \int \hat\psi^\dagger({\bf r}) \left ( - \;
\frac{\nabla^2}{2m} + U \right ) \hat\psi({\bf r}) \; d{\bf r}
$$
\begin{equation}
\label{51}
+ \; \frac{1}{2}\; \int \hat\psi^\dagger({\bf r}) \hat\psi^\dagger({\bf r}')
\Phi({\bf r}-{\bf r}') \hat\psi({\bf r}')\hat\psi({\bf r}) \; d{\bf r} d{\bf r}' \; ,
\end{equation}
in which $\hat\psi({\bf r})=\hat\psi({\bf r},t)$ is the shifted field operator
defined in Eq. (16). In order to avoid cumbersome notations, we omit
the explicit dependence on time of the variables $\eta({\bf r})=\eta({\bf r},t)$
and $\psi_1({\bf r})=\psi_1({\bf r},t)$, though this dependence is assumed. In
what follows, such an abreviated notation will be often used, where it
does not lead to confusion. The external potential $U=U({\bf r},t)$ can,
in general, depend on time. The interaction potential $\Phi({\bf r})$ is
arbitrary, provided it is symmetric, such that $\Phi({\bf r})=\Phi(-{\bf r})$,
and enjoys the Fourier transformation.

The grand Hamiltonian (30) can be written as a sum of five terms,
\begin{equation}
\label{52}
H[\eta,\psi_1] = \sum_{n=0}^4 H^{(n)} \; ,
\end{equation}
whose order is labelled by the number of operators $\psi_1$ in their
products. The zero-order term contains no $\psi_1$, being
$$
H^{(0)} = \int \eta^*({\bf r}) \left ( -\; \frac{\nabla^2}{2m} + U -
\mu_0 \right ) \eta({\bf r}) \; d{\bf r}
$$
\begin{equation}
\label{53}
+ \; \frac{1}{2} \int \Phi({\bf r}-{\bf r}') |\eta({\bf r}')|^2 |\eta({\bf r})|^2 \;
d{\bf r} d{\bf r}' \; .
\end{equation}
To satisfy the quantum-number conservation condition (27), the Hamiltonian
must not contain the terms linear in $\psi_1$. This prescribes to take the
Lagrange multipliers in Eq. (26) as
\begin{equation}
\label{54}
\lambda({\bf r},t) = \left [ - \; \frac{\nabla^2}{2m} + U +
\int \Phi({\bf r}-{\bf r}') |\eta({\bf r}',t)|^2\; d{\bf r}' \right ] \eta({\bf r},t) \; ,
\end{equation}
so that $H^{(1)}=0$. The necessity of removing the terms linear in $\psi_1$
or $\psi_1^\dagger$, by means of condition (54), in order to satisfy constraint
(25), can be easily proved for quadratic Hamiltonians, involving linear
terms, by employing the method of canonical transformations [42,45]. In the
general case of arbitrary Hamiltonians, the proof of Eq. (54), yielding the
cancellation of linear terms, is presented in the Appendix B.

The second-order term is
$$
H^{(2)} = \int \psi_1^\dagger({\bf r}) \left ( - \; \frac{\nabla^2}{2m} +
U - \mu_1 \right ) \psi_1({\bf r})\; d{\bf r}
$$
$$
+ \; \int \Phi({\bf r}-{\bf r}') \left [
|\eta({\bf r})|^2 \psi_1^\dagger({\bf r}') \psi_1({\bf r}') +
\eta^*({\bf r})\eta({\bf r}')\psi_1^\dagger({\bf r}')\psi_1({\bf r}) \right.
$$
\begin{equation}
\label{55}
+ \; \left. \frac{1}{2}\; \eta^*({\bf r}) \eta^*({\bf r}')
\psi_1({\bf r}')\psi_1({\bf r}) + \frac{1}{2}\; \eta({\bf r})\eta({\bf r}')
\psi_1^\dagger({\bf r}') \psi_1^\dagger({\bf r}) \right ] \; d{\bf r} d{\bf r}' \; .
\end{equation}
Respectively, we have the third-order term
\begin{equation}
\label{56}
H^{(3)} = \int \Phi({\bf r}-{\bf r}') \left [
\eta^*({\bf r})\psi_1^\dagger({\bf r}')\psi_1({\bf r}')\psi_1({\bf r}) +
\psi_1^\dagger({\bf r}) \psi_1^\dagger({\bf r}') \psi_1({\bf r}') \eta({\bf r}) \right ]\;
d{\bf r} d{\bf r}'
\end{equation}
and the fourth-order term
\begin{equation}
\label{57}
H^{(4)} = \frac{1}{2} \; \int \psi_1^\dagger({\bf r}) \psi_1^\dagger({\bf r}')
\Phi({\bf r}-{\bf r}') \psi_1({\bf r}')\psi_1({\bf r})\; d{\bf r} d{\bf r}' \; .
\end{equation}

From Eq. (34), we obtain
\begin{equation}
\label{58}
i\; \frac{\partial}{\partial t} \; \eta({\bf r},t) = \left ( - \;
\frac{\nabla^2}{2m} + U - \mu_0 \right ) \eta({\bf r},t) +
\int \Phi({\bf r}-{\bf r}') \left [ |\eta({\bf r}')|^2 \eta({\bf r}) +
\hat X({\bf r},{\bf r}') \right ] \; d{\bf r}' \; ,
\end{equation}
where the last term in the integrand is the correlation operator
$$
\hat X({\bf r},{\bf r}') \equiv \psi_1^\dagger({\bf r}')\psi_1({\bf r}')\eta({\bf r})
$$
\begin{equation}
\label{59}
+ \; \psi_1^\dagger({\bf r}')\eta({\bf r}')\psi_1({\bf r}) +
\eta^*({\bf r}')\psi_1({\bf r}')\psi_1({\bf r}) +
\psi_1^\dagger({\bf r}')\psi_1({\bf r}')\psi_1({\bf r}) \; .
\end{equation}
And Eq. (35) yields
$$
i\; \frac{\partial}{\partial t}\; \psi_1({\bf r},t) = \left ( - \;
\frac{\nabla^2}{2m} + U - \mu_1 \right )\psi_1({\bf r},t)
$$
\begin{equation}
\label{60}
+\; \int \Phi({\bf r}-{\bf r}') \left [ |\eta({\bf r}')|^2 \psi_1({\bf r}) +
\eta^*({\bf r}')\eta({\bf r})\psi_1({\bf r}') +
\eta({\bf r}')\eta({\bf r})\psi_1^\dagger({\bf r}') + \hat X({\bf r},{\bf r}')\right ]\;
d{\bf r}' \; .
\end{equation}
Again, for brevity, the time dependence is not explicitly shown in the
integrals of Eqs. (59) and (60). Also, recall that in the operator equation
(58), the condensate function is assumed to be factored with $\hat 1_{\cal F}$,
the unity operator in ${\cal F}(\psi_1)$.

\section{Local Conservation Laws}

The equations of motion (58) and (60) are derived from the variational
principle of action extremization. Therefore they guarantee the validity
of all local conservation laws on the operator level. As an illustration,
let us consider the time variation of the local densities.

The local condensate density is
\begin{equation}
\label{61}
\rho_0({\bf r},t) \equiv |\eta({\bf r},t)|^2 \; ,
\end{equation}
and the local condensate current is
\begin{equation}
\label{62}
{\bf j}_0({\bf r},t) \equiv -\; \frac{i}{2m} \left [ \eta^*({\bf r},t)\;
\nabla\eta({\bf r},t) - \eta({\bf r},t)\; \nabla\eta^*({\bf r},t)\right ] \; .
\end{equation}
Respectively, we define the operator density of uncondensed particles
\begin{equation}
\label{63}
\hat\rho_1({\bf r},t) \equiv \psi_1^\dagger({\bf r},t)\psi_1({\bf r},t)
\end{equation}
and the related current operator
\begin{equation}
\label{64}
{\bf j}_1({\bf r},t) \equiv -\; \frac{i}{2m} \left\{ \psi_1^\dagger({\bf r},t)\;
\nabla\psi_1({\bf r},t) - \left [\nabla\psi_1^\dagger({\bf r},t)\right ]
\psi_1({\bf r},t)\right \} \; .
\end{equation}

Differentiating the condensate density (61), we find
\begin{equation}
\label{65}
\frac{\partial}{\partial t}\; \rho_0({\bf r},t) + \nabla\cdot{\bf j}_0({\bf r},t) =
- \hat \Gamma({\bf r},t) \; ,
\end{equation}
where the source operator is
\begin{equation}
\label{66}
\hat\Gamma({\bf r},t) \equiv \int \Phi({\bf r}-{\bf r}') \left [
\hat R({\bf r},{\bf r}') + \hat R^+({\bf r},{\bf r}')\right ]\; d{\bf r}' \; ,
\end{equation}
in which
$$
\hat R({\bf r},{\bf r}') \equiv i\eta^*({\bf r}) \left [
\psi_1^\dagger({\bf r}')\eta({\bf r}') + \eta^*({\bf r}')\psi_1({\bf r}') +
\psi_1^\dagger({\bf r}')\psi_1({\bf r}')\right ] \psi_1({\bf r}) \; .
$$
And the time variation of density (63) yields
\begin{equation}
\label{67}
\frac{\partial}{\partial t} \; \hat\rho_1({\bf r},t) +
\nabla\cdot\hat{\bf j}_1({\bf r},t) = \hat\Gamma({\bf r},t) \; .
\end{equation}
For the total operator density
\begin{equation}
\label{68}
\hat\rho({\bf r},t) = \rho_0({\bf r},t) + \hat\rho_1({\bf r},t)
\end{equation}
and the total operator of current
\begin{equation}
\label{69}
\hat{\bf j}({\bf r},t) = {\bf j}_0({\bf r},t) + \hat{\bf j}_1({\bf r},t) \; ,
\end{equation}
we obtain the {\it continuity equation}
\begin{equation}
\label{70}
\frac{\partial}{\partial t}\; \rho({\bf r},t) +
\nabla\cdot \hat{\bf j}({\bf r},t) = 0 \; .
\end{equation}
In the same way, one can derive any other local conservation laws, following
the standard procedure [40,45] and employing the equations of motion (58)
and (60).

\section{Condensate Wave Function}

The equation for the condensate wave function $\eta({\bf r},t)$ follows from
averaging Eq. (58) over ${\cal F}(\psi_1)$. To this end, we need to introduce
several notations. The normal density matrix is
\begin{equation}
\label{71}
\rho_1({\bf r},{\bf r}') \; \equiv \; < \psi_1^\dagger({\bf r}')\psi_1({\bf r})> \; .
\end{equation}
Under the broken gauge symmetry, there appears the so-called anomalous
density matrix
\begin{equation}
\label{72}
\sigma_1({\bf r},{\bf r}') \; \equiv \; <\psi_1({\bf r}') \psi_1({\bf r})> \; .
\end{equation}
The density of condensed particles is $\rho_0({\bf r})$, defined in Eq. (61),
and the density of uncondensed particles is
\begin{equation}
\label{73}
\rho_1({\bf r}) \equiv \rho_1({\bf r},{\bf r}) \; = \;
<\psi_1^\dagger({\bf r})\psi_1({\bf r})> \; .
\end{equation}
We shall also need the anomalous average
\begin{equation}
\label{74}
\sigma_1({\bf r}) \equiv \sigma_1({\bf r},{\bf r}) \; = \;
< \psi_1({\bf r}) \psi_1({\bf r})> \; .
\end{equation}
The absolute value $|\sigma_1({\bf r})|$ has the meaning of the density of
pair-correlated particles [43]. The total density
\begin{equation}
\label{75}
\rho({\bf r}) = \rho_0({\bf r}) + \rho_1({\bf r})
\end{equation}
is the average of the operator density (68). According to the normalization
conditions (20) and (21), the partial densities are normalized to the number
of corresponding particles,
\begin{equation}
\label{76}
N_0 = \int \rho_0({\bf r})\; d{\bf r} \; , \qquad
N_1 = \int \rho_1({\bf r})\; d{\bf r} \; .
\end{equation}
One more notation is the triple correlator
\begin{equation}
\label{77}
\xi({\bf r},{\bf r}') \; \equiv \;
<\psi_1^\dagger({\bf r}')\psi_1({\bf r}')\psi_1({\bf r})> \; .
\end{equation}
Using these definitions, the average of the correlation operator (59)
becomes
$$
<\hat X({\bf r},{\bf r}')> \; = \; \rho_1({\bf r}')\eta({\bf r}) +
\rho_1({\bf r},{\bf r}')\eta({\bf r}') + \sigma_1({\bf r},{\bf r}') \eta^*({\bf r}') +
\xi({\bf r},{\bf r}') \; .
$$

Finally, averaging Eq. (58) yields the evolution equation for the
condensate wave function
$$
i\; \frac{\partial}{\partial t}\; \eta({\bf r},t) = \left ( -\;
\frac{\nabla^2}{2m} + U - \mu_0 \right ) \eta({\bf r},t)
$$
\begin{equation}
\label{78}
+\; \int \Phi({\bf r}-{\bf r}') \left [ \rho({\bf r}')\eta({\bf r}) +
\rho_1({\bf r},{\bf r}')\eta({\bf r}') + \sigma_1({\bf r},{\bf r}')\eta^*({\bf r}') +
\xi({\bf r},{\bf r}')\right ]\; d{\bf r}' \; .
\end{equation}
This is the general equation for an arbitrary Bose-condensed system. This
equation can also be represented in a shorter, though not explicit, form
$$
i\; \frac{\partial}{\partial t}\; \eta({\bf r},t) = \left ( -\;
\frac{\nabla^2}{2m} + U - \mu_0 \right ) \eta({\bf r},t) +
\int \Phi({\bf r}-{\bf r}')<\hat\psi^\dagger({\bf r}')\hat\psi({\bf r}')\hat\psi({\bf r})>
d{\bf r}' \; ,
$$
in which one has to substitute the shifted operator $\hat\psi=\eta+\psi_1$.

We may notice that Eq. (78) is invariant under the transformation
\begin{equation}
\label{79}
\eta({\bf r},t) \; \longrightarrow \eta({\bf r},t)e^{i\alpha t} \; , \qquad
\psi_1({\bf r},t) \; \longrightarrow \psi_1({\bf r},t)e^{i\alpha t} \; , \qquad
\mu_0 \; \longrightarrow \mu_0 + \alpha \; .
\end{equation}
This means that there exists a freedom in choosing the phase factor
$\exp(i\alpha t)$ of the condensate function. But the phase factor becomes
fixed by defining the Lagrange multiplier $\mu_0$. The arbitrariness in
the condensate phase implies that the stationary solutions for $\eta({\bf r},t)$
would be proportional to an undefined factor $e^{i\alpha t}$. By fixing the
multiplier $\mu_0$, we require that in equilibrium the condensate function
would not be dependent on time, that is,
\begin{equation}
\label{80}
\frac{\partial}{\partial t}\; \eta({\bf r}) = 0 \qquad ({\rm equilibrium}) \; .
\end{equation}

In equilibrium, Eq. (78), according to condition (80), reduces to the
eigenvalue problem
$$
\mu_0\eta({\bf r}) = \left [ -\; \frac{\nabla^2}{2m} + U({\bf r})\right ]
\eta({\bf r})
$$
\begin{equation}
\label{81}
+ \; \int \Phi({\bf r}-{\bf r}') \left [ \rho({\bf r}')\eta({\bf r}) +
\rho_1({\bf r},{\bf r}')\eta({\bf r}') + \sigma_1({\bf r},{\bf r}')\eta^*({\bf r}') +
\xi({\bf r},{\bf r}') \right ] \; d{\bf r}' \; ,
\end{equation}
defining, together with the normalization condition (76), the eigenfunction
$\eta({\bf r})$ and the eigenvalue $\mu_0$. Let us emphasize that without the
normalization condition (76) the condensate function cannot be uniquely
defined from Eq. (81). For example, if the system is uniform, with no
external potential $U\rightarrow 0$, when $\eta({\bf r})=\eta$, $\rho({\bf r})=\rho$, and
$\eta=\sqrt{\rho_0}$, then Eq. (81) gives
\begin{equation}
\label{82}
\mu_0 = \rho \Phi_0 + \int \Phi({\bf r}) \left [ \rho_1({\bf r},0) +
\sigma_1({\bf r},0) + \frac{\xi({\bf r},0)}{\sqrt{\rho_0}} \right ] \; d{\bf r} \; ,
\end{equation}
where $\Phi_0\equiv\int\Phi({\bf r})d{\bf r}$, or in the compact form
$$
\mu_0 = \frac{1}{\sqrt{\rho_0}} \; \int \Phi({\bf r})
<\hat\psi^\dagger(0)\hat\psi(0)\hat\psi({\bf r})>d{\bf r} \; .
$$
Expression (82) is valid for an arbitrary uniform equilibrium system with
an interaction potential $\Phi({\bf r})$. No approximations have been used in
obtaining Eq. (82).

In recent days, the physics of dilute Bose gases is intensively explored
(see the book [46] and review articles [31,47--49]). The interaction
potential for these gases is taken in the contact form
\begin{equation}
\label{83}
\Phi({\bf r}) = \Phi_0\delta({\bf r}) \; , \qquad
\Phi_0 \equiv 4\pi\; \frac{a_s}{m} \; ,
\end{equation}
where $a_s$ is the scattering length. With this interaction potential,
the eigenvalue problem (81) for a nonuniform system becomes
$$
\mu_0 \eta({\bf r}) = \left [ -\; \frac{\nabla^2}{2m} + U({\bf r}) \right ]
\eta({\bf r})
$$
\begin{equation}
\label{84}
+ \; \Phi_0 \left\{ \left [ \rho({\bf r}) + \rho_1({\bf r})\right ]\eta +
\sigma_1({\bf r})\eta^*({\bf r}) + \xi({\bf r},{\bf r})\right \} \; .
\end{equation}
For a uniform system, $\mu_0$ is given by Eq. (82), which, in the case
of the contact potential (83), yields
\begin{equation}
\label{85}
\mu_0 = \left ( \rho +\rho_1 +\sigma_1 +
\frac{\xi}{\sqrt{\rho_0}} \right ) \Phi_0 \; ,
\end{equation}
where $\sigma_1\equiv\sigma_1({\bf r},{\bf r})$ and $\xi\equiv\xi({\bf r},{\bf r})$.

When one assumes that the Bose gas is dilute, being characterized by
the contact potential (83), the temperature is zero, and the interaction
is so weak that the condensate depletion can be neglected, so that all
particles are condensed, then $N_0=N$, $N_1=0$, and $\rho_1=\sigma_1=\xi=0$.
In this approximation, because of relation (44), the multiplier $\mu_0=\mu$
coincides with the system chemical potential. The eigenvalue problem (84)
reduces then to the Gross-Pitaevskii equation
$$
\left [ -\; \frac{\nabla^2}{2m} + U({\bf r}) + \Phi_0|\eta({\bf r})|^2
\right ] \eta({\bf r}) = \mu \eta({\bf r}) \; .
$$
But when the condensate depletion is not negligible, then $\mu_0$,
according to Eq. (44), is not the same as $\mu$.

\section{Uniform Bose system}

Let us illustrate in more detail the application of the representative
ensemble with the grand Hamiltonian (52), to an equilibrium uniform
system. Then the field operators of uncondensed particles can be expanded
in plane waves,
\begin{equation}
\label{86}
\psi_1({\bf r}) =\sum_{k\neq 0} a_k \varphi_k({\bf r}) \; ,
\end{equation}
where $\varphi_k({\bf r})=e^{i{\bf k}\cdot{\bf r}}/\sqrt{V}$. The condensate function
reduces to the constant $\eta=\sqrt{\rho_0}$, while the condensate
multiplier $\mu_0$ is given by Eq. (82).

In the momentum representation, the momentum distribution of particles
\begin{equation}
\label{87}
n_k \equiv \; <a_k^\dagger a_k >
\end{equation}
is usually termed the normal average, as compared to the anomalous average
\begin{equation}
\label{88}
\sigma_k \equiv \; < a_k a_{-k} > \; .
\end{equation}
The normal and anomalous density matrices (71) and (72) take the form of
the expansions
\begin{equation}
\label{89}
\rho_1({\bf r},{\bf r}') =\sum_{k\neq 0} n_k \varphi_k({\bf r}) \varphi_k^*({\bf r}') \; ,
\qquad
\sigma_1({\bf r},{\bf r}') =\sum_{k\neq 0} \sigma_k \varphi_k({\bf r}) \varphi_k^*({\bf r}') \; ,
\end{equation}
in which the properties
$$
<a_k^\dagger a_p> \; = \; \delta_{kp} n_k \; , \qquad
<a_k a_p> \; = \; \delta_{-kp}\sigma_k
$$
are taken into account. The diagonal elements of the matrices in Eq. (89)
give the densities
\begin{equation}
\label{90}
\rho_1 = \rho_1({\bf r},{\bf r}) = \frac{1}{V} \; \sum_{k\neq 0} n_k \; ,
\qquad
\sigma_1 = \sigma_1({\bf r},{\bf r}) = \frac{1}{V} \; \sum_{k\neq 0} \sigma_k \; .
\end{equation}
The interaction potential is assumed to allow the Fourier transformation
\begin{equation}
\label{91}
\Phi({\bf r}) = \frac{1}{V} \; \sum_k \Phi_k e^{i{\bf k}\cdot{\bf r}} \; \qquad
\Phi_k = \int \Phi({\bf r}) e^{-i{\bf k}\cdot{\bf r}} \; d{\bf r} \; .
\end{equation}
Then we find the following terms of the grand Hamiltonian (52). The
zero-order term (53) becomes
\begin{equation}
\label{92}
H^{(0)} = \left ( \frac{1}{2}\; \rho_0 \Phi_0 - \mu_0 \right ) N_0 \; .
\end{equation}
The first-order term $H^{1)}=0$ is automatically zero. The second-order
term (55) is
\begin{equation}
\label{93}
H^{(2)} = \sum_{k\neq 0} \left [ \frac{k^2}{2m} +
\rho_0 (\Phi_0+\Phi_k) - \mu_1 \right ]  a_k^\dagger a_k +
\frac{1}{2} \; \sum_{k\neq 0} \rho_0 \Phi_k \left (
a_k^\dagger a_{-k}^\dagger + a_{-k} a_k \right ) \; .
\end{equation}
For the third-order term (56), we have
\begin{equation}
\label{94}
H^{(3)} = \sqrt{\frac{\rho_0}{V}} \; {\sum_{kp}}' \; \Phi_p
\left ( a_k^\dagger a_{k+p} a_{-p} + a_{-p}^\dagger a_{k+p}^\dagger a_k
\right ) \; ,
\end{equation}
where the prime on the summation symbol implies that ${\bf k}\neq 0$,
${\bf p}\neq 0$, ${\bf k}+{\bf p}\neq 0$. The fourth-order term (57) is
\begin{equation}
\label{95}
H^{(4)} = \frac{1}{2V} \;
\sum_q {\sum_{kp}}' \; \Phi_q a_k^\dagger a_p^\dagger a_{p+q} a_{k-q} \; ,
\end{equation}
where the prime shows that ${\bf k}\neq 0$, ${\bf p}\neq 0$, ${\bf p}+{\bf q}\neq 0$,
${\bf k}-{\bf q}\neq 0$.

To be more specific, it is necessary to resort to some approximation. The
natural mean-field approximation for a system with broken gauge symmetry
is the HFB approximation. The latter is usually blamed to display an
unphysical gap in the spectrum, because of which it is qualified as gapful
(see detailed discussion in Refs. [7,31,32]). However, as is explained in
the Introduction, this defect comes into play only owing to the usage of
a nonrepresentative statistical ensemble. But for the representative
ensemble, with the grand Hamiltonian (52), there appear no such problems.
Below, we show this for the case of arbitrary temperature and interaction
potential $\Phi({\bf r})$.

We apply, in the standard way [40,45], the HFB approximation to the third-
and fourth-order products of the operators $a_k$. Then the third-order
term (94) becomes identically zero, because of the condition $<a_k>=0$,
\begin{equation}
\label{96}
H^{(3)} = 0 \; .
\end{equation}
And for the fourth-order term (95), we get
$$
H^{(4)} = \sum_{k\neq 0} \rho_1 \Phi_0 \left ( a_k^\dagger a_k \; - \;
\frac{1}{2}\; n_k \right )
$$
\begin{equation}
\label{97}
+\; \frac{1}{V} \; \sum_{k,p\neq 0} \Phi_k \left [
n_{k+p} a_p^\dagger a_p + \frac{1}{2} \left (
\sigma_{k+p} a_p^\dagger a_{-p}^\dagger + \sigma_{k+p}^* a_{-p} a_p \right ) \;
- \; \frac{1}{2}\left ( n_{k+p} n_p + \sigma_{k+p} \sigma_p \right )
\right ] \; .
\end{equation}

It is convenient to introduce the notations
\begin{equation}
\label{98}
\omega_k \equiv \frac{k^2}{2m} + \rho \Phi_0 + \rho_0 \Phi_k +
\frac{1}{V} \; \sum_{p\neq 0} n_p \Phi_{k+p} - \mu_1
\end{equation}
and
\begin{equation}
\label{99}
\Delta_k \equiv \rho_0 \Phi_k + \frac{1}{V} \sum_{p\neq 0}
\sigma_p \Phi_{k+p} \; .
\end{equation}
Since the interaction potential $\Phi({\bf r})=\Phi(-{\bf r})$ is symmetric
and real, we have the properties $\sigma_k=\sigma_k^*=\sigma_{-k}$ and
$\Delta_k=\Delta_k^*=\Delta_{-k}$.

Summing up all terms in the grand Hamiltonian (52), we obtain in the HFB
approximation
\begin{equation}
\label{100}
H_{HFB} =  E_{HFB} + \sum_{k\neq 0} \omega_k a_k^\dagger a_k +
\frac{1}{2} \sum_{k\neq 0} \Delta_k \left ( a_k^\dagger a_{-k}^\dagger +
a_{-k} a_k \right ) \; ,
\end{equation}
where the nonoperator term is
\begin{equation}
\label{101}
E_{HFB} =  H^{(0)} \; - \; \frac{1}{2} \; \rho_1^2 \Phi_0 V \; - \;
\frac{1}{2V} \; \sum_{k,p\neq 0} \Phi_{k+p} \left ( n_k n_p +
\sigma_k \sigma_p \right ) \; .
\end{equation}
The quadratic form (100) can be diagonalized by the Bogolubov canonical
transformation [3,4,45]
$$
a_k = u_k b_k + v_{-k}^* b_{-k}^\dagger \; .
$$
As a result, the grand Hamiltonian (100) is transformed to the Bogolubov
representation
\begin{equation}
\label{102}
H_B = E_B + \sum_{k\neq 0} \varepsilon_k b_k^\dagger b_k \; ,
\end{equation}
in which
\begin{equation}
\label{103}
E_B = E_{HFB}  + \frac{1}{2}\; \sum_{k\neq 0} (\varepsilon_k - \omega_k) \; .
\end{equation}
For the operators $b_k$, one has the properties
\begin{equation}
\label{104}
<b_k>\; = \; <b_k b_p> \; = \; 0 \; , \qquad
<b_k^\dagger b_p> \; = \; \delta_{kp} \pi_k \; ,
\end{equation}
with the momentum distribution of quasiparticles
\begin{equation}
\label{105}
\pi_k \equiv \; < b_k^\dagger b_k > \; = \; \left ( e^{\beta\varepsilon_k} -
1 \right )^{-1} = \frac{1}{2} \left [ {\rm coth}\left (
\frac{\varepsilon_k}{2T}\right ) -1 \right ] \; .
\end{equation}
The coefficient functions $u_k$ and $v_k$, and the spectrum $\varepsilon_k$ are
defined by the Bogolubov - de Gennes equations
$$
(\omega_k -\varepsilon_k ) u_k + \Delta_k v_k = 0 \; , \qquad
\Delta_k u_k + (\omega_k + \varepsilon_k ) v_k = 0 \; ,
$$
with the normalization condition $u_k^2-v_k^2=1$. This gives
$$
u_k^2  = \frac{\omega_k+\varepsilon_k}{2\varepsilon_k} \; , \qquad
v_k^2 = \frac{\omega-\varepsilon_k}{2\varepsilon_k} \; .
$$
And the Bogolubov spectrum is
\begin{equation}
\label{106}
\varepsilon_k = \sqrt{\omega_k^2 - \Delta_k^2} \; .
\end{equation}

The existence of the Bose-Einstein condensate, as is known, requires
that the spectrum (106) be gapless, such that
\begin{equation}
\label{107}
\lim_{k\rightarrow 0} \varepsilon_k = 0 \; ,
\end{equation}
with the stability condition $\varepsilon_k\geq 0$. Then Eq. (106) yields
\begin{equation}
\label{108}
\mu_1 = \rho \Phi_0 + \frac{1}{V} \;
\sum_{p\neq 0} ( n_p - \sigma_p) \Phi_p \; .
\end{equation}
As will be shown in Sec. 10, this value is in exact agreement with the
Hugenholtz-Pines relation. This should be compared with the condensate
potential (82), which in the HFB approximation, when $\xi({\bf r},{\bf r}')=0$,
becomes
\begin{equation}
\label{109}
\mu_0 = \rho \Phi_0 + \frac{1}{V} \;
\sum_{p\neq 0} ( n_p + \sigma_p) \Phi_p \; .
\end{equation}
In the particular case of the contact potential (83), we get
\begin{equation}
\label{110}
\mu_1 = (\rho + \rho_1 -\sigma_1) \Phi_0
\end{equation}
and, respectively,
\begin{equation}
\label{111}
\mu_0 = (\rho + \rho_1 + \sigma_1) \Phi_0 \; .
\end{equation}
In any case, generally, $\mu_0\neq\mu_1$. The multipliers $\mu_0$
and $\mu_1$ become equal only in the limit of the asymptotically small
condensate depletion, which corresponds to the Bogolubov approximation
[1,2], when both $\rho_1$ and $\sigma_1$ in Eqs. (108) to (111) are
neglected. Then, clearly, $\mu_0=\mu_1=\rho_0\Phi_0$, and in view of
relation (44), $\mu_0=\mu_1=\mu$. However, as soon as the condensate
depletion is not negligible, $\mu_0\neq\mu_1$.

With the multiplier (108), Eq. (98) can be represented as
\begin{equation}
\label{112}
\omega_k = \frac{k^2}{2m} + \rho_0 \Phi_k + \frac{1}{V} \;
\sum_{p\neq 0} \left ( n_p \Phi_{k+p} - n_p \Phi_p +
\sigma_p \Phi_p \right ) \; .
\end{equation}
For the long-wave limit of $\Delta_k$, we introduce the notation
\begin{equation}
\label{113}
\Delta \equiv \lim_{k\rightarrow 0} \Delta_k \equiv mc^2 \; .
\end{equation}
From Eq. (99) it follows
\begin{equation}
\label{114}
\Delta = \rho_0 \Phi_0 + \frac{1}{V} \;
\sum_{p\neq 0} \sigma_p \Phi_p \; .
\end{equation}
Then the long-wave limit of spectrum (106) is explicitly of the phonon
type, $\varepsilon\simeq ck$, as $k\rightarrow 0$. According to the notation (113), the
sound velocity is
\begin{equation}
\label{115}
c =\sqrt{\frac{\Delta}{m} } \; .
\end{equation}
Strictly speaking, the sound velocity is defined as $c=\sqrt{\Delta/m^*}$,
with an effective mass to be given in Sec. 10. For short-range
interactions, $m^*\cong m$. Using Eqs. (113) and (114), we may rewrite
Eq. (112) in the form
\begin{equation}
\label{116}
\omega_k = \frac{k^2}{2m} + mc^2 + \rho_0 (\Phi_k -\Phi_0) +
\int n_p (\Phi_{k+p} - \Phi_p) \; \frac{d{\bf k}}{(2\pi)^3} \; .
\end{equation}
For the particle momentum distribution (87), we find
\begin{equation}
\label{117}
n_k = \frac{\omega_k}{2\varepsilon_k} \; {\rm coth}\left (
\frac{\varepsilon_k}{2T} \right ) \; - \; \frac{1}{2} \; ,
\end{equation}
while for the anomalous average (88), we obtain
\begin{equation}
\label{118}
\sigma_k = -\; \frac{\Delta_k}{2\varepsilon_k} \; {\rm coth} \left (
\frac{\varepsilon_k}{2T} \right ) \; .
\end{equation}
Analyzing the behaviour of $n_k$ and $\sigma_k$ as functions of momentum
$k$, we find [15,32] that $|\sigma_k|\simeq n_k$ for $k\rightarrow 0$, while
$|\sigma_k|\gg n_k$ for large $k$. Therefore in no sense the anomalous
average $\sigma_k$ can be neglected, as compared to the normal average
$n_k$.

The derived equations simplify for the contact potential (83). Then
$\Phi_k=\Phi_0$, $\Delta_k=\Delta$,
$$
\Delta=(\rho_0 + \sigma_1) \Phi_0 \; ,
$$
and spectrum (106) takes the classical Bogolubov form
\begin{equation}
\label{119}
\varepsilon_k = \sqrt{(ck)^2 + \left ( \frac{k^2}{2m}\right )^2 } \; .
\end{equation}

The grand potential (36) in the HFB approximation is
\begin{equation}
\label{120}
\Omega = E_B + TV \int \ln \left ( 1  - e^{-\beta\varepsilon_k}
\right ) \; \frac{d{\bf k}}{(2\pi)^3 } \; .
\end{equation}
From this and other equations, obtained above, all thermodynamic
characteristics can be calculated.

\section{Condensate and Superfluid Fractions}

We shall concentrate our attention on the most interesting characteristics
of the Bose-condensed system, on the condensate and superfluid fractions.
The condensate fraction
\begin{equation}
\label{121}
n_0 =  1 \; - \; \frac{\rho_1}{\rho}
\end{equation}
can be found by calculating the density $\rho_1$ of uncondensed particles,
given by the integral
\begin{equation}
\label{122}
\rho_1 = \int n_k \; \frac{d{\bf k}}{(2\pi)^3} \; .
\end{equation}

To define the superfluid fraction, one considers the reaction of the
system to the boost with a velocity ${\bf v}$. In the laboratory frame,
the field operators of a moving system are represented by means of the
Galilean transformation
\begin{equation}
\label{123}
\hat\psi_v({\bf r},t) = \hat\psi({\bf r}-{\bf v} t,t) \exp\left \{ i \left (
m{\bf v} \cdot{\bf r} \; - \; \frac{mv^2}{2}\; t\right )\right \} \; .
\end{equation}
The Hamiltonian $H_v=H[\hat\psi_v]$, in terms of the new field operators
(123), is expressed through the Hamiltonian $H=H[\hat\psi]$, in terms of
the old operators $\hat\psi$, as
\begin{equation}
\label{124}
H_v = H + \int \hat\psi^\dagger({\bf r}) \left ( - i{\bf v}\cdot\nabla +
\frac{mv^2}{2} \right ) \hat\psi({\bf r})\; d{\bf r} \; .
\end{equation}
The total system momentum operator becomes
\begin{equation}
\label{125}
\hat{\bf P}_v = \hat{\bf P} + Nm{\bf v} \; ,
\end{equation}
where
\begin{equation}
\label{126}
\hat{\bf P}  = \int \hat\psi^\dagger({\bf r}) (-i\nabla)\hat\psi({\bf r})\; d{\bf r}
\end{equation}
is the old momentum operator. The average total momentum of the moving
system is
\begin{equation}
\label{127}
<\hat{\bf P}_v>_v \; \equiv \;
{\rm Tr}\hat\rho_v \hat{\bf P}_v \; ,
\end{equation}
with the statistical operator
\begin{equation}
\label{128}
\hat\rho_v \equiv \frac{e^{-\beta H_v}}{{\rm Tr}e^{-\beta H_v}} \; .
\end{equation}

The superfluid fraction is defined as
\begin{equation}
\label{129}
n_s \equiv \frac{1}{3mN}\; \lim_{v\rightarrow 0} \; \frac{\partial}{\partial{\bf v}}\;
\cdot <\hat{\bf P}_v>_v \; .
\end{equation}
It is possible to show (see, e.g., the detailed derivation in Refs.
[47,49]) that the superfluid fraction (129) for an arbitrary system can be
represented in the form
\begin{equation}
\label{130}
n_s =  1 \; - \; \frac{2Q}{3T} \; ,
\end{equation}
in which
\begin{equation}
\label{131}
Q \equiv \frac{\Delta^2(\hat{\bf P})}{2mN}
\end{equation}
is the {\it dispersed heat}, with the momentum dispersion
$$
\Delta^2(\hat{\bf P}) \equiv \; <\hat{\bf P}^2> - < \hat{\bf P}>^2 \; ,
$$
where the averages are with respect to the grand Hamiltonian $H[\hat\psi]$.

Formula (130) for the superfluid fraction is valid for any system,
equilibrium or nonequilibrium, uniform or not. For a system in equilibrium,
the total momentum is zero, $<\hat{\bf P}>=0$, so that the dispersed heat
(131) is
\begin{equation}
\label{132}
Q = \frac{<\hat{\bf P}^2>}{2mN} \; .
\end{equation}
For a uniform system,
\begin{equation}
\label{133}
<\hat{\bf P}^2> \; = \; \sum_{kp} \left ({\bf k} \cdot {\bf p} \right )
<\hat n_k \hat n_p> \; ,
\end{equation}
where $\hat n_k\equiv a_k^\dagger a_k$. In the HFB approximation,
\begin{equation}
\label{134}
<\hat n_k \hat n_p> \; = \; n_k n_p + \delta_{kp} n_k (1 + n_k) +
\delta_{-kp} \sigma_k^2 \; .
\end{equation}
Then the dispersed heat (132) becomes
\begin{equation}
\label{135}
Q = \frac{1}{\rho} \; \int \frac{k^2}{2m} \left (
n_k + n_k^2 - \sigma_k^2 \right ) \frac{d{\bf k}}{(2\pi)^3} \; .
\end{equation}
Substituting here expressions (117) and (118), and using the equality
${\rm cosh}x-1=2{\rm sinh}^2(x/2)$, so that
$$
n_k + n_k ^2 - \sigma_k^2 = \frac{1}{4{\rm sinh}^2(\beta\varepsilon_k/2)} \; ,
$$
we get from Eq. (135)
\begin{equation}
\label{136}
Q = \frac{1}{(4\pi)^2m\rho} \; \int_0^\infty \;
\frac{k^4\; dk}{{\sinh}^2(\beta\varepsilon_k/2)} \; .
\end{equation}
Note again the importance of the anomalous average $\sigma_k$. If one would
omit the latter in the dispersed heat (135), one would get a senseless
divergent quantity, while accurately taking account of $\sigma_k$ results in
the well-defined convergent integral (136).

To demonstrate explicitly the behaviour of the condensate and superfluid
fractions, let us resort to the contact potential (83). Then Eq. (99) yields
\begin{equation}
\label{137}
\Delta_k =\Delta = mc^2 \; ,
\end{equation}
and Eq. (116) gives
\begin{equation}
\label{138}
\omega_k = \frac{k^2}{2m} + mc^2 \; .
\end{equation}
Equation (115), defining the sound velocity, can be written as
\begin{equation}
\label{139}
mc^2 = (\rho_0 +\sigma_1) \Phi_0 \; ,
\end{equation}
where
\begin{equation}
\label{140}
\sigma_1 = \int \sigma_k \; \frac{d{\bf k}}{(2\pi)^3} \; .
\end{equation}

The density of uncondensed particles (122) can be represented in the form
\begin{equation}
\label{141}
\rho_1 = \frac{(mc)^3}{3\pi^2} \left \{ 1 + \frac{3}{2\sqrt{2}} \;
\int_0^\infty \; \left (\sqrt{1+x^2}-1\right )^{1/2} \left [
{\rm coth}\left ( \frac{mc^2}{2T}\; x\right ) -1 \right ]\; dx
\right \} \; ,
\end{equation}
which is a well-defined convergent integral.

The anomalous average (140) can be written as a sum of two parts
\begin{equation}
\label{142}
\sigma_1 = \sigma_0  - \int \; \frac{mc^2}{2\varepsilon_k} \left [
{\rm coth}\left ( \frac{\varepsilon_k}{2T} \right ) -1 \right ] \;
\frac{d{\bf k}}{(2\pi)^3} \; ,
\end{equation}
in which
\begin{equation}
\label{143}
\sigma_0 \equiv - \Delta \int \; \frac{1}{2\varepsilon_k} \;
\frac{d{\bf k}}{(2\pi)^3} \; .
\end{equation}
The second integral in Eq. (142) is convergent. But the integral in Eq. (143)
ultravioletly diverges. However, this divergence is known to be unphysical,
being simply caused by the usage of the contact interaction potential. This
and other similar divergences could be removed by employing more realistic
interaction potentials. For example, a common choice is a Gaussian type
potential [50]. Another known way of removing such divergences is by using
the analytic regularization in one of its variants, such as the subtraction
scheme, zeta regularization, or dimensional regularization (see details in
Ref. [31]). This kind of regularization is asymptotically exact for weak
interactions and is universal in the sense that it applies to any short-range
potential with a scattering length $a_s$ [31]. Using the dimensional
regularization in its region of validity, we obtain
$$
\int \frac{1}{2\varepsilon_k} \; \frac{d{\bf k}}{(2\pi)^3} =  - \;
\frac{m}{\pi^2}\; \sqrt{m\rho_0\Phi_0} \; .
$$
Then Eq. (143) results in
\begin{equation}
\label{144}
\sigma_0 = \frac{(mc)^2}{\pi^2} \; \sqrt{m\rho_0\Phi_0} \; .
\end{equation}
The anomalous average (142) can be represented as
\begin{equation}
\label{145}
\sigma_1 = \sigma_0 \; - \; \frac{(mc)^3}{2\sqrt{2}\;\pi^2} \;
\int_0^\infty \; \frac{(\sqrt{1+x^2}-1)^{1/2}}{\sqrt{1+x^2}} \left [
{\rm coth}\left (\frac{mc^2}{2T}\; x \right ) -1 \right ] \; dx \; .
\end{equation}
And for the dispersed heat (136), we obtain
\begin{equation}
\label{146}
Q = \frac{(mc)^5}{\sqrt{2}\; (2\pi)^2 m\rho} \; \int_0^\infty \;
\frac{(\sqrt{1+x^2}-1)^{3/2}\; xdx}
{\sqrt{1+x^2}\; {\rm sinh}^2(mc^2x/2T)} \; .
\end{equation}

At low temperatures, such that
$$
\frac{T}{mc^2} \ll 1 \; ,
$$
equations (141), (145), and (146) yield
$$
\rho_1 \simeq \frac{(mc)^3}{3\pi^2} + \frac{(mc)^3}{12} \; \left (
\frac{T}{mc^2}\right )^2 \; , \qquad
\sigma_1 \simeq \sigma_0 \; - \; \frac{(mc)^3}{12}\left (
\frac{T}{mc^2}\right )^2 \; ,
$$
$$
Q \simeq \frac{\pi^2(mc)^5}{15m\rho} \; \left (
\frac{T}{mc^2}\right )^5 \; .
$$
Therefore, the condensate fraction (121) behaves as
\begin{equation}
\label{147}
n_0 \simeq 1 \; - \; \frac{(mc)^3}{3\pi^2\rho} \; - \;
\frac{(mc)^3}{12\rho}\; \left ( \frac{T}{mc^2} \right )^2 \; .
\end{equation}
And the superfluid fraction (130) becomes
\begin{equation}
\label{148}
n_s \simeq 1 \; - \; \frac{2\pi^2(mc)^3}{45\rho} \;
\left ( \frac{T}{mc^2}\right )^4 \; .
\end{equation}
When the interaction is weak, or the condensate fraction $n_0$ tends to
zero at the critical temperature $T_c$, as a results of which, $c\rightarrow 0$,
so that
$$
\frac{mc^2}{T_c} \ll 1 \; ,
$$
then Eqs. (141), (145), and (146) lead to
$$
\rho_1 \simeq \rho \left ( \frac{T}{T_c} \right )^{3/2} +
\frac{(mc)^3}{3\pi^2} \; , \qquad \sigma_1 \simeq \sigma_0 \; - \;
\frac{m^2cT}{2\pi} \; ,
$$
$$
Q \simeq \frac{3}{2}\; T \; \left [ \left ( \frac{T}{T_c}
\right )^{3/2} - \; \frac{\zeta(1/2)}{\zeta(3/2)} \left (
\frac{T}{T_c} \right )^{1/2}
\frac{mc^2}{T_c} \right ] \; ,
$$
with the critical temperature
$$
T_c = \frac{2\pi}{m} \left [ \frac{\rho}{\zeta(3/2)}\right ]^{2/3}
$$
coinciding with that for the ideal Bose gas, as it should be in the case
of the mean-field approximation with the contact potential (83). Here
$\zeta(\cdot)$ is the Riemann zeta function. Then the condensate fraction
tends to zero at $T_c$ as
\begin{equation}
\label{149}
n_0 \simeq 1 \; - \; \left ( \frac{T}{T_c} \right )^{3/2} - \;
\frac{(mc)^3}{3\pi^2\rho} \; ,
\end{equation}
together with the superfluid fraction
\begin{equation}
\label{150}
n_s \simeq 1 \; - \; \left ( \frac{T}{T_c}\right )^{3/2} + \;
\frac{\zeta(1/2)}{\zeta(3/2)}\; \left ( \frac{T}{T_c}\right )^{1/2}
\frac{mc^2}{T_c} \; .
\end{equation}
Both, $n_0$ and $n_s$, tend to zero from the left, demonstrating the
second-order phase transition at $T_c$. The point why in the present case
$T_c$ has to be the same as the ideal Bose gas critical temperature is
explained in the Appendix C.

\section{Green Function Equations}

In order to conclude the proof of the complete self-consistency of the
developed approach, based on the introduced representative ensemble, we need
to show that for an arbitrary Bose system there exist the Green-function
equations of the usual form and that the Hugenholtz-Pines relation follows
for any uniform system.

The first-order Green function is defined in the standard way [3,4,7,9] as the
matrix $G(12)=[G_{\alpha\beta}(12)]$ with the elements
$$
G_{11}(12) = - i<\hat T\psi_1(1)\psi_1^\dagger(2)> \; , \qquad
G_{12}(12) = - i<\hat T\psi_1(1)\psi_1(2)> \; ,
$$
\begin{equation}
\label{151}
G_{21}(12) = - i<\hat T\psi_1^\dagger(1)\psi_1^\dagger(2)> \; , \qquad
G_{22}(12) = - i<\hat T\psi_1^\dagger(1)\psi_1(2)> \; ,
\end{equation}
in which the common abbreviation is employed, denoting the set $\{{\bf r}_j,t_j\}$
by a single number $j$, and $\hat T$ is the chronological operator.

For what follows, we shall need the notation for the triple field operator
$$
\Psi(123) \equiv \psi_1(1)\psi_1(2)\psi_1^\dagger(3) +
\eta(1)\psi_1(2)\psi_1^\dagger(3) + \psi_1(1)\eta(2)\psi_1^\dagger(3)
$$
\begin{equation}
\label{152}
+\; \psi_1(1)\psi_1(2)\eta^*(3) +
\eta(1)\eta(2)\psi_1^\dagger(3) + \eta(1)\psi_1(2)\eta^*(3) +
\psi_1(1)\eta(2)\eta^*(3) \; .
\end{equation}
One may notice that
$$
\Psi(123) = \hat\psi(1) \hat\psi(2)\hat\psi^\dagger(3)  -
\eta(1)\eta(2)\eta^*(3) \; ,
$$
where $\hat\psi=\eta+\psi_1$ is the shifted field operator. However,
for practical application, we need the explicit form (152).

The second-order Green function is a matrix
$B(1234)=[B_{\alpha\beta}(1234)]$ with the elements
$$
B_{11}(1234) = -<\hat T\Psi(123)\psi_1^\dagger(4)> \; , \qquad
B_{12}(1234) = -<\hat T\Psi(123)\psi_1(4)> \; ,
$$
\begin{equation}
\label{153}
B_{21}(1234) = -<\hat T\Psi^\dagger(123)\psi_1^\dagger(4)> \; , \qquad
B_{22}(1234) = -<\hat T\Psi^\dagger(123)\psi_1(4)> \; .
\end{equation}

For the general form of a retarded interaction potential
\begin{equation}
\label{154}
\Phi(12) \equiv \Phi({\bf r}_1 - {\bf r}_2) \delta(t_{12}+0) \; ,
\end{equation}
in which $t_{12}=t_1-t_2$, the self-energy is introduced through the equation
\begin{equation}
\label{155}
\int \Sigma(13) G(32)\; d(3) = i \int \Phi(13) B(1332)\; d(3) \; .
\end{equation}
The equation for $G(12)$, by defining the inverse propagator
\begin{equation}
\label{156}
G^{-1}(12) \equiv \left ( \hat\tau \; i \;
\frac{\partial}{\partial t_1} + \frac{\nabla_1^2}{2m} \; - \; U(1) +
\mu_1 \right ) \delta(12) \; \hat 1 - \Sigma(12) \; ,
\end{equation}
in which
\begin{eqnarray}
\nonumber
\hat\tau \equiv \left [ \begin{array}{rr}
1 & 0 \\
0 & -1 \end{array}  \right ]  \; , \qquad
\hat 1 \equiv \left [ \begin{array}{rr}
1 & 0 \\
0 & 1  \end{array} \right ] \; ,
\end{eqnarray}
can be represented in the form
\begin{equation}
\label{157}
\int G^{-1}(13) G(32) \; d(3) = \delta(12)\; \hat 1 \; .
\end{equation}

To solve Eq. (157), one may invoke perturbation theory, starting with
an available approximate Green function $G_{app}$, corresponding to an
approximate $\Sigma_{app}$, such that
\begin{equation}
\label{158}
\int G_{app}^{-1}(13) G_{app}(32)\; d(3) = \delta(12)\hat 1 \; .
\end{equation}
Then Eqs. (157) and (158) can be transformed to the Dyson representation
\begin{equation}
\label{159}
G(12) = G_{app}(12) + \int G_{app}(13) \left [
\Sigma(34) - \Sigma_{app}(34)\right ] G(42)\; d(34) \; ,
\end{equation}
which is convenient for using perturbation theory.

The above equations for the Green functions are valid for an arbitrary
nonequilibrium and nonuniform Bose system. If the considered system is
equilibrium and uniform, one passes to the Fourier transforms $G({\bf k},\omega)$
and $\Sigma({\bf k},\omega)$, employing the symmetry properties
$$
G_{11}({\bf k},-\omega) = G_{22}({\bf k},\omega) \; , \quad
G_{12}({\bf k},-\omega) = G_{21}({\bf k},\omega) = G_{12}({\bf k},\omega) \; ,
$$
\begin{equation}
\label{160}
G_{\alpha\beta}(-{\bf k},\omega) = G_{\alpha\beta}({\bf k},\omega)
\end{equation}
and, respectively,
$$
\Sigma_{11}({\bf k},-\omega) = \Sigma_{22}({\bf k},\omega) \; , \qquad
\Sigma_{12}({\bf k},-\omega) =\Sigma_{21}({\bf k},\omega) =\Sigma_{12}({\bf k},\omega) \; ,
$$
\begin{equation}
\label{161}
\Sigma_{\alpha\beta}(-{\bf k},\omega) =\Sigma_{\alpha\beta}({\bf k},\omega) \; .
\end{equation}
A detailed discussion of these relations can be found in Ref. [4].

Let us introduce the notation
\begin{equation}
\label{162}
D({\bf k},\omega) \equiv \Sigma_{12}^2 ({\bf k},\omega) -
G_{11}^{-1}({\bf k},\omega) G_{11}^{-1}({\bf k},-\omega) \; ,
\end{equation}
in which
$$
G_{11}^{-1}({\bf k},\omega) = \omega \; - \; \frac{k^2}{2m} \; - \;
\Sigma_{11}({\bf k},\omega) + \mu_1 \; .
$$
Then the solution to Eq. (157) can be written as
\begin{equation}
\label{163}
G_{11}({\bf k},\omega) =
\frac{\omega+k^2/2m+\Sigma_{11}({\bf k},\omega)-\mu_1}{D({\bf k},\omega)} \; , \qquad
G_{12}({\bf k},\omega) = -\; \frac{\Sigma_{12}({\bf k},\omega)}{D({\bf k},\omega)} \; ,
\end{equation}
with the denominator (162).

For the Green functions $G_{\alpha\beta}({\bf k},0)$ at zero energy there is the
Bogolubov theorem [4] rigorously proving the validity of the inequalities
\begin{equation}
\label{164}
| G_{11}({\bf k},0)| \geq \frac{mn_0}{2k^2} \; ,
\end{equation}
\begin{equation}
\label{165}
| G_{11}({\bf k},0) - G_{12}({\bf k},0) | \geq \frac{mn_0}{k^2} \; .
\end{equation}
From Eq. (163), we have
$$
G_{11}({\bf k},0) =
\frac{k^2/2m+\Sigma_{11}({\bf k},0)-\mu_1}{D({\bf k},0)} \; , \qquad
G_{12}({\bf k},0) = -\; \frac{\Sigma_{12}({\bf k},0)}{D({\bf k},0)} \; ,
$$
where
$$
D({\bf k},0) = \Sigma_{12}^2({\bf k},0) - \left [ \frac{k^2}{2m} +
\Sigma_{11}({\bf k},0) - \mu_1 \right ]^2 \; .
$$
Consequently,
$$
G_{11}({\bf k},0) - G_{12}({\bf k},0) = \left [ \mu_1 \; - \; \frac{k^2}{2m} +
\Sigma_{12}({\bf k},0) - \Sigma_{11}({\bf k},0) \right ]^{-1} \; .
$$
Using this in Eq. (165), we get
$$
\left | \mu_1 \; - \; \frac{k^2}{2m} + \Sigma_{12}({\bf k},0) -
\Sigma_{11}({\bf k},0) \right | \leq \frac{k^2}{mn_0} \; .
$$
Setting in the latter equation ${\bf k}\rightarrow 0$, we come to the Hugenholtz-Pines
relation
\begin{equation}
\label{166}
\mu_1 = \Sigma_{11}(0,0) - \Sigma_{12}(0,0) \; .
\end{equation}

Note that in the HFB approximation we have
$$
\Sigma_{11}(0,0) = (\rho+\rho_0)\Phi_0 + \frac{1}{V}\;
\sum_{p\neq 0} n_p \Phi_p \; , \qquad
\Sigma_{12}(0,0) = \rho_0 \Phi_0 + \frac{1}{V} \;
\sum_{p\neq 0} \sigma_p\Phi_p \; .
$$
Therefore Eq. (166) gives exactly the form of $\mu_1$ in Eq. (108).

Relation (166) guarantees that the system spectrum is gapless and of the
phonon character. The spectrum $\varepsilon_k$ is given by the zeros of the Green
functions (163), that is, by the equation
\begin{equation}
\label{167}
D({\bf k},\varepsilon_k) = 0 \; ,
\end{equation}
with Eq. (162). Equation (167) can be rewritten as
\begin{equation}
\label{168}
\varepsilon_k = \frac{1}{2} \; \left [ \Sigma_{11}({\bf k},\varepsilon_k) -
\Sigma_{22}({\bf k},\varepsilon_k) \right ] +
\sqrt{\omega_k^2- \Sigma_{12}^2({\bf k},\varepsilon_k)} \; ,
\end{equation}
where the notation
$$
\omega_k \equiv \frac{k^2}{2m} + \frac{1}{2} \; \left [
\Sigma_{11}({\bf k},\varepsilon_k) + \Sigma_{22}({\bf k},\varepsilon_k) \right ] - \mu_1
$$
is used. From Eq. (168), it follows that $\varepsilon_k\rightarrow 0$, as $k\rightarrow 0$. Moreover,
when the system is isotropic, its self-energy $\Sigma_{\alpha\beta}({\bf k},\varepsilon_k)$
depends only on the scalar $k^2$. This implies the asymptotic, as $k\rightarrow 0$,
expansion
$$
\Sigma_{\alpha\beta}({\bf k},\varepsilon_k) \simeq \Sigma_{\alpha\beta}(0,0) +
{\Sigma}'_{\alpha\beta}\;  k^2 \; ,
$$
in which
$$
{\Sigma'}_{\alpha\beta} \equiv \lim_{k\rightarrow 0} \; \frac{\partial}{\partial k^2} \;
\Sigma_{\alpha\beta}({\bf k},\varepsilon_k) \; .
$$
Using this expansion in spectrum (168) gives
\begin{equation}
\label{169}
\varepsilon_k \simeq ck \; , \qquad
c \equiv \sqrt{\frac{1}{m^*}\; \Sigma_{12}(0,0)} \; ,
\end{equation}
which is the phonon spectrum, with the sound velocity $c$, where the
effective mass is
\begin{equation}
\label{170}
m^* \equiv \frac{m}{1+m(\Sigma'_{11}+\Sigma'_{22}-2\Sigma'_{12})} \; .
\end{equation}

It is useful to compare the Lagrange multipliers $\mu_0$ and $\mu_1$
given by their general expressions (82) and (166), which are valid for
an arbitrary equilibrium uniform Bose system. These expressions are exact,
with no approximations being involved. As is seen, anomalous averages enter
$\mu_0$ with the sign plus, while the anomalous self-energy enters $\mu_1$
with the sign minus. This is why, in general, $\mu_0$ cannot coincide with
$\mu_1$. The HFB forms (108) and (109), or (110) and (111), are particular
illustrations. It is easy to check that $\mu_0$ coincides with $\mu_1$
solely in the Bogolubov limit of asymptotically weak interactions, when
$\Sigma_{11}(0,0)\rightarrow(\rho+\rho_0)\Phi_0$, $\Sigma_{12}(0,0)\rightarrow\rho_0\Phi_0$,
and $\mu_0\rightarrow\mu_1\rightarrow\mu=\rho\Phi_0$. But in any higher-order approximation,
$\mu_0\neq\mu_1$. The assumption that $\mu_0$ would be the same as $\mu_1$
would make the theory not self-consistent and would return us back to the
Hohenberg-Martin dilemma of conserving versus gapless theories.

\section{Discussion}

The main message of this paper is the necessity of employing representative
ensembles for correctly describing statistical systems. A representative
ensemble takes into account all imposed constraints and additional conditions
that allow for a unique description of the considered system. It is only using
a representative ensemble makes the theory self-consistent.

In the Bose system with broken global gauge symmetry, realized by the
Bogolubov shift, there are two particle components, corresponding to condensed
and uncondensed particles, with two related normalization conditions for $N_0$
and $N_1$. This requires to introduce two Lagrange multipliers, $\mu_0$ and
$\mu_1$, which makes the theory completely self-consistent in any approximation.

It is worth recalling that the introduction of several Lagrange multipliers
is rather common for spin systems. There, the order parameter is the average
spin, which, generally, is a three-component vector. The role of the effective
chemical potential for spin systems is played, as is well known, by an external
magnetic field, which is also a three-component vector. Hence, the number of
effective chemical potentials for spin systems is equal to the number of
components in the order parameter. Only then one is able to unambiguously
define the average spin.

The suggested approach, introducing two Lagrange multipliers, does not
contradict our physical understanding that the standard experiments fix, as
independent variables, temperature, volume, and the total number of particles.
To emphasize this once again, let us turn to the definition of the grand
thermodynamic potential (36), from which it follows that it is a function
$\Omega = \Omega(T,V,\mu_0,\mu_1)$, so that the free energy (42) is a function
$F=F(T,V,N_0,N_1)$. The Lagrange multiplier $\mu_0$ is defined from the
stability condition, yielding $\mu_0 = \mu_0(T,V,N_0,N_1)$, according to
Eqs. (81), (82), (84), (85), (109), and (111). The Lagrange multipler
$\mu_1$ satisfies the Hugenholtz-Pines theorem, which gives
$\mu_1=\mu_1(T,V,N_0,N_1)$, in agreement with Eqs. (108), (110), and (166).
The number of uncondensed particles is found from the direct calculation of
the average (21) expressed through Eqs. (90), (122), (141), and like that,
resulting in $N_1= N_1(T,V,N_0,N)$. From here, since $N_0=N-N_1$, it follows
that $N_0 = N_0(T,V,N)$. Substituting it back to $N_1$, one has
$N_1=N_1(T,V,N)$. Using these in the expressions for $\mu_0$ and $\mu_1$,
one gets $\mu_0=\mu_0(T,V,N)$ and $\mu_1=\mu_1(T,V,N)$. Then, from Eq. (44),
it is evident that $\mu=\mu(T,V,N)$, or, inverting the latter relation, one
has $N=N(T,V,\mu)$. Using this, the Lagrange multipliers $\mu_0$ and $\mu_1$
can be expressed as functions $\mu_0=\mu_0(T,V,\mu)$ and $\mu_1=\mu_1(T,V,\mu)$.
Substituting this into the grand potential, we have $\Omega = \Omega(T,V,\mu)$,
in line with Eq. (48). Respectively, the free energy becomes a function
$F=F(T,V,N)$, in accordance with Eq. (45).

Thus, at the end, we work with the standard variables $T,V$, and $N$,
which are usually fixed in experiments. All observable quantities are
also expressed through these variables. So that the suggested approach
is absolutely self-consistent, mathematically correct, and in agreement
with physics.

Several words are to be said with regard to the phase-transition order
of Bose-Einstein condensation. This transition is known to be of {\it
second-order}, which is firmly based on several facts. First, there exists
a general explanation, independent of the coupling strength, that this
transition is of second-order. This can be found in the book by Patashinsky
and Pokrovsky [51] (Chapter X, Section 2). As is also well known, the
Hamiltonian of Bose systems is mathematically equivalent to what in quantum
theory is termed the $\varphi^4$ model. The phase transition in this model has
been studied in numerous works using the renormalization group approach,
exhibiting the second-order transition [52]. The superfluid transition in
liquid $^4$He, which is believed to be accompanied by the Bose-Einstein
condensation, is also a second-order transition. A large body of experimental
data on measuring the continuous temperature dependence of the condensate
fraction in superfluid helium has been summarized by Wirth and Hallock [53].
There exists abundant literature, both theoretical and experimental, on
Bose-Einstein condensation in dilute trapped gases (see review works
[46--49]) and references therein) and there are several computer simulations
of this process [54--56]. Though Bose condensation in traps is smeared out
by finite-size effects, the subsequent increase of the number of particles
unambiguously demonstrates that the condensation approaches the standard
second-order phase transition. In addition, there have been many Monte Carlo
calculations for uniform Bose systems with various interaction potentials.
All these calculations, summarized in the review articles [57--59], clearly
prove that Bose-Einstein condensation is a second-order transition. So, the
second order of the Bose condensate transition has been established without
any doubt. This especially concerns the general theoretical investigations
[51,52] and rigorous Monte Carlo calculations [57,58].

The main idea of the present paper is the necessity of using representative
ensembles for describing Bose-condensed systems. This implies that proper
allowance must be made for all conditions which uniquely define the employed
field variables. Here the consideration has been based on the classical
Bogolubov approach [1--4] introducing two field variables, the condensate
function $\eta$ and the operator of uncondensed particles $\psi_1$. By their
definition, these variables are {\it independent} of each other. For a uniform
system, the condensate corresponds to the zero-momentum state, while this
state is excluded from the description of uncondensed particles. This is
evident from definition (86) of $\psi_1$, where ${\bf k}\neq 0$. The variables
$\eta$ and $\psi_1$ are also {\it orthogonal} to each other, in agreement with
definition (17), which becomes obvious from Eq. (86), since
\begin{equation}
\label{171}
\frac{1}{\sqrt{V}}\; \int \psi_1({\bf r})\; d{\bf r} = \sum_{k\neq 0}
a_k \delta_{k0} = 0 \; .
\end{equation}
For so introduced independent orthogonal variables, it is necessary to
define two normalization conditions and, respectively, two Lagrange
multipliers controlling these normalization conditions. Only then there can
be the assurance that the theory will be self-consistent in any calculations.

Of course, if the field variables are introduced in a different way, with
some other conditions, this would require to define another ensemble, with
the appropriate Lagrange multipliers, whose number could also be different.
For example, Hugenholtz and Pines [9], as well as later Gavoret and Nozieres
[60], when deriving the Hugenholtz-Pines relation on the basis of
thermodynamic properties, {\it did not use} the standard grand ensemble.
Any attentive reader can immediately infer from the original works [9,60]
that these authors have used a different ensemble. They treat a uniform
equilibrium system, defining the number of condensed particles $N_0$ by
minimizing the {\it internal energy} at zero temperature, which corresponds
to the minimization of {\it free energy} at finite temperatures. After
defining in that way $N_0=N_0(\rho,T)$, the latter is explicitly substituted
into the effective Hamiltonian $\hat H-\mu\hat N_1$, where $\hat N_1$ is the
operator of {\it uncondensed} particles, but not of the total number of
particles, as it would be in the standard grand Hamiltonian. Since
$N_0=N_0(\rho,T)$ has been explicitly substituted everywhere, one needs
the sole Lagrange multiplier for the normalization of {\it uncondensed}
particles. But mathematically this is absolutely equivalent to the
introduction of an additional Lagrange multiplier $\mu_0$, as is done
in the present paper, {\it before} substituting $N_0(\rho,T)$, which is
defined later from the corresponding normalization condition. These ways,
as is absolutely clear, are equivalent, but the latter method is more
convenient for more general cases of nonuniform or nonequilibrium systems.

Another example is given by the approach advanced by Faddeev and Popov [61]
and used later by Popov [12--14]. They define the field operator
\begin{equation}
\label{172}
\psi({\bf r}) =\sqrt{\rho_0} + \psi_{FP}({\bf r}) \; ,
\end{equation}
in which the second part in the right-hand side,
\begin{equation}
\label{173}
\psi_{FP}({\bf r}) = \frac{1}{\sqrt{V}}\; \sum_k a_k e^{i{\bf k}\cdot{\bf r}} \; ,
\end{equation}
includes the term $a_0\neq 0$ with ${\bf k}=0$. This is contrary to the Bogolubov
field operator of uncondensed particles (86), not containing this
zero-momentum term. The Fourier transform $c_k$ of the operator
$$
\psi({\bf r}) = \frac{1}{\sqrt{V}} \; \sum_k c_k e^{i{\bf k}\cdot{\bf r}}
$$
is connected with the Fourier transform of $\psi_{FP}$ by the equation
$$
c_k = \delta_{0k}\; \sqrt{N_0} + a_k \; .
$$
For the zero-momentum state, one has
$$
c_0 =\sqrt{N_0} + a_0 \; .
$$
Recall that in the Bogolubov case, one would have $c_0=\sqrt{N_0}$. Thus,
in the Faddeev-Popov approach, the condensate is not completely separated from
uncondensed particles, but $\psi_{FP}$ does contain the zero-momentum term
$a_0$. In other words, $\psi_{FP}$ is not independent from $\rho_0$. Moreover,
the Faddeev-Popov representation (172) for $\psi$ consists of two parts that
are not orthogonal to each other,
$$
\frac{1}{\sqrt{V}} \; \int \psi_{FP}({\bf r})\; d{\bf r} = a_0 \neq 0 \; ,
$$
which is contrary to the Bogolubov case (171). Hence, $\psi_{FP}$ is neither
independent of $\sqrt{\rho_0}$ no orthogonal to it, but both of them define
the sole variable (172). As far as, in the Faddeev-Popov approach, $a_0\neq
0$, the latter yields the interference terms in physical operators. For
instance, here the number-of-particle operator becomes
\begin{equation}
\label{174}
\hat N \equiv \int \psi^\dagger({\bf r})\psi({\bf r})\; d{\bf r} = N_0 +
\sum_k a_k^\dagger a_k +\sqrt{N_0}\left ( a_0^\dagger + a_0\right ) \; .
\end{equation}
This operator is normalized to the total number of particles, so that
\begin{equation}
\label{175}
N \;  =\; <\hat N> \; = \; N_0 + \sum_k <a_k^\dagger a_k> \; ,
\end{equation}
which requires $<a_0>=0$. Faddeev and Popov in their original paper [61]
emphasized that their representation (172) is principally different from the
Bogolubov shift (16) and discussed in detail the corresponding differences.
Because in the Faddeev-Popov approach there is a single independent variable
$\psi$, with the normalization condition (175), the appropriate representative
ensemble here is the standard grand ensemble, with the grand Hamiltonian $\hat
H-\mu\hat N$, though the number-of-particle operator (174) here is different
from the Bogolubov form (24). More complicated forms of physical operators and
the necessity to comply with relation (172) at each step of any calculational
procedure make the usage of the Faddeev-Popov approach more complicated and
quite inconvenient for mean-field type approximations. However, needless to say
that a theory with one independent variable and, respectively, with one
Lagrange multiplier, is mathematically equivalent to the theory with two
independent variables and two Lagrange multipliers. Which representation to
choose is rather a matter of convenience.

It is even admissible to introduce no Lagrange multipliers and to deal with
the canonical Gibbs ensemble. But the problem with the latter is that then one
has to keep the total number of particles fixed not merely on the average but
exactly at each step of any calculational procedure. This requires to work in
a restricted Fock space, where the number -of-particle operator degenerates to
the number $\hat N\equiv N$. To accomplish this, one may resort to the
Girardeau-Arnowitt approach [24,25] introducing the so-called
number-conserving field operators
$$
\alpha_k = \hat N_0^{-1/2} a_0^\dagger a_k\; , \qquad
\alpha_k^\dagger = a_k^\dagger a_0 \hat N_0^{-1/2} \; .
$$
This representation, however, is valid only when the number of condensed
particles is large, $N_0\gg 1$, so that it is not applicable in the vicinity
of the condensation point, when $N_0\rightarrow 0$. Also, such a canonical
representation in approximate calculations yields, as is known [24], a gap in
the spectrum. It was noticed by Girardeau [62] and showed by Takano [63] that
to get a self-consistent theory in the canonical ensemble requires to use the
whole Hamiltonian, without reducing it to approximate forms. But the theory
with a whole Hamiltonian for interacting particles has no exact solution. So,
the usage of the canonical ensemble is unpractical in analytic investigations,
though it may be employed for numerical computations [64,65].

Concluding, the choice of an ensemble is, generally speaking, a matter of
convenience. But in any case, the chosen ensemble must be representative,
which necessitates to accurately take into account all conditions uniquely
defining the considered system. The choice of an ensemble is intimately
connected with the choice of the appropriate field variables, which should
not be confused. Employing inappropriate variables, that is, using a
nonrepresentative ensemble may result in the appearance of inconsistences and
paradoxes. For example, resorting to the canonical ensemble, one should use
the Girardeau-Arnowitt representation with the number-conserving field
operators [24,25] or the Carusotto-Castin-Dalibard representation based on
stochastic fields [64,65]. If one prefers the standard grand ensemble, then
the Faddeev-Popov representation is appropriate [61]. A nonstandard variant of
the grand ensemble, due to Hugenholtz and Pines [9], can also be used. But
when the classical Bogolubov representation [1--4] is chosen, where the field
variables of condensed and uncondensed particles are separated from each
other, being independent and orthogonal, then the approach developed in the
present paper is the most convenient, being completely self-consistent.

\vskip 5mm

{\bf Acknowledgement}. I am grateful to the German Research Foundation
for financial support and to the Physics Departments of the Free University
of Berlin and University of Essen for hospitality. I appreciate useful
discussions with M. Girardeau, R. Graham, H. Kleinert, and E. Yukalova.

\newpage

{\Large{\bf Appendix A}}

\vskip 5mm

Here we prove relation (50). From definition (36) for the grand potential
with Hamiltonian (30), it immediately follows that
$$
N_0 = -\; \frac{\partial\Omega}{\partial\mu_0} \; , \qquad
N_1 = -\; \frac{\partial\Omega}{\partial\mu_1} \; .
$$
On the other hand,
$$
N= - \; \frac{\partial\Omega}{\partial\mu} = N_0 + N_1 \; .
$$
Comparing the above equations, we get
$$
\frac{\partial\Omega}{\partial\mu} = \frac{\partial\Omega}{\partial\mu_0} +
\frac{\partial\Omega}{\partial\mu_1} \; .
$$
Differentiating again the latter equation, we have
$$
\frac{\partial^2\Omega}{\partial\mu^2} =
\frac{\partial^2\Omega}{\partial\mu_0^2} + \frac{\partial^2\Omega}{\partial\mu_1^2} +
2\; \frac{\partial^2\Omega}{\partial\mu_0\partial\mu_1} \; .
$$
By direct calculations, we find that
$$
\frac{\partial^2\Omega}{\partial\mu_0^2} = -\beta\Delta^2(\hat N_0) \; , \qquad
\frac{\partial^2\Omega}{\partial\mu_1^2} = -\beta\Delta^2(\hat N_1) \; , \qquad
\frac{\partial^2\Omega}{\partial\mu_0\partial\mu_1} =
- \beta\; {\rm cov}(\hat N_0,\hat N_1) \; ,
$$
where the covariance of two operators, $\hat A$ and $\hat B$, is
$$
{\rm cov}(\hat A,\hat B) \equiv \frac{1}{2}<\hat A\hat B +
\hat B\hat A> - <\hat A><\hat B> \; .
$$
Summarizing these equations, and using the property of the dispersion
of a composite operator,
$$
\Delta^2(\hat A+\hat B) = \Delta^2(\hat A) + \Delta^2(\hat B) +
2{\rm cov}(\hat A,\hat B) \; ,
$$
we obtain
$$
\frac{\partial N}{\partial\mu} = -\; \frac{\partial^2\Omega}{\partial\mu^2} =
\beta\Delta^2(\hat N) \; ,
$$
which proves relation (50).

\newpage

{\Large{\bf Appendix B}}

\vskip 5mm

The fact that the terms linear in $\psi_1$ or $\psi_1^\dagger$ must be absent
in the Hamiltonian, in order to satisfy constraint (25), can be proved in
two ways.

We may, first, consider a quadratic Hamiltonian approximating the exact one.
Linear terms in such a Hamiltonian can also be present. Quadratic Hamiltonians
of this type can be diagonalized with the help of exact canonical
transformations [45]. Rigorous mathematical properties of the corresponding
transformations, called nonuniform, are expounded in the book by Berezin [42].
After the Hamiltonian is diagonalized, it is straightforward to calculate
explicitly all averages, which show that the linear terms in the Hamiltonian
induce nonzero $<\psi_1>$ and, vice versa, zero linear terms lead to
$<\psi_1>=0$. Then one should consider perturbation theory, starting with
the diagonalized quadratic from, and check that zero linear terms yield
$<\psi_1>=0$ in all orders of the theory. This way is rather cumbersome, and
below another method is presented.

Let us consider the general Hamiltonian (52). Its term, linear in $\psi_1$
and $\psi_1^\dagger$, is
$$
H^{(1)} = \int \psi_1^\dagger({\bf r}) \left ( - \; \frac{\nabla^2}{2m} + U
\right ) \eta({\bf r})\; d{\bf r} +
\int \psi_1({\bf r}) \left ( - \; \frac{\nabla^2}{2m} + U
\right ) \eta^*({\bf r})\; d{\bf r}
$$
$$
+\; \int \Phi({\bf r}-{\bf r}') \left [
\psi_1^\dagger({\bf r})|\eta({\bf r}')|^2 \eta({\bf r}) +
\eta^*({\bf r})|\eta({\bf r}')|^2 \psi_1({\bf r}) \right ]\; d{\bf r} d{\bf r}' \; .
$$
The equation of motion (35) can be represented as
$$
i\; \frac{\partial}{\partial t} \; \psi_1({\bf r},t) =
\frac{\delta[H - H^{(1)}]}{\delta\psi_1^\dagger({\bf r},t)} +
\frac{\delta H^{(1)}}{\delta\psi_1^\dagger({\bf r},t)} \; .
$$
The first term here gives the right-hand side of Eq. (60). The second term
results in the form
$$
\frac{\delta H^{(1)}}{\delta\psi_1^\dagger({\bf r},t)} =
C({\bf r},t)\; \hat 1_{\cal F} \; ,
$$
proportional to the unit operator $\hat 1_{\cal F}$ in the Fock space
${\cal F}(\psi_1)$ generated by the field operators $\psi_1^\dagger$, with the
nonoperator complex function
$$
C({\bf r},t) \equiv \left [ -\; \frac{\nabla^2}{2m} + U +
\int \Phi({\bf r}-{\bf r}') |\eta({\bf r}',t)|^2 d{\bf r}'\right ] \eta({\bf r},t) -
\lambda({\bf r},t) \; .
$$
The equation of motion, written above, is an operator equality defined on
${\cal F}(\psi_1)$. The operator equality assumes that it holds true at least in
the weak sense implying the equality of all matrix elements for the states
from the space the operators are defined on. In the present case, the latter
means the space ${\cal F}(\psi_1)$. The vacuum state of this space is defined in
the standard way as the vector $|0>$, for which $\psi_1({\bf r},t)|0>=0$.
Considering for the equation of motion, with respect to $\psi_1$, the matrix
element over the vacuum state, and, taking into account constraint (25), we
get
$$
C({\bf r},t)= 0 \; .
$$
The latter equation results in the Lagrange multiplier (54), which yields
$H^{(1)}=0$. Thus, in order to preserve constraint (25), the linear terms
in the Hamiltonian must be zero.

The physical meaning of the proved theorem, that $<\psi_1>=0$ necessarily
requires $H^{(1)}=0$, is quite evident. The terms in the Hamiltonian, linear
in $\psi_1$ and $\psi_1^\dagger$ describe the physical processes of annihilation
and creation of single particles. If such processes were permitted, then the
function $C({\bf r},t)$, introduced above, is not zero. Then the evolution
equation for the average $<\psi_1({\bf r},t)>$ contains the term
$$
\int_0^t  C({\bf r},t')\; dt'
$$
generating the nonzero value of $<\psi_1>$.

One should not confuse the considered situation with the often used method
of specially adding to the Hamiltonian the terms, linear in $\psi_1$ and
$\psi_1^\dagger$, in order to break the gauge symmetry. Then, of course, the
average $<\psi_1>$ is not zero. But at the end of calculations, one always
set the additional linear terms to become zero, hence restoring the property
$<\psi_1>=0$. This procedure is what is called the Bogolubov method of
infinitesimal sources [3,4]. These sources are usually lifted after the
thermodynamic limit, but can also be made infinitesimally small in the
process of taking the thermodynamic limit, provided this is done in the
appropriate way [5,66]. In the approach, followed in the present paper, we
do not need to introduce infinitesimal sources, since the gauge symmetry has
already been broken by the Bogolubov shift (16).

\newpage

{\Large{\bf Appendix C}}

\vskip 5mm

The critical temperature $T_c$ for the Bose system with contact interactions
has been obtained, in the HFB approximation, from expansions for the condensate
fraction (149) and superfluid fraction (150). This temperature was found to
coincide with the critical temperature of the ideal Bose gas. In order to
emphasize the correctness of this result, let us consider, first, the more
general case of an arbitrary interaction potential $\Phi({\bf r})$, provided it
possesses the standard property of symmetry, such that $\Phi(-{\bf r})=\Phi({\bf r})$,
and diminishes sufficiently fast with increasing $|{\bf r}|$.

At the critical temperature $T_c$, when $\rho_0=0$ and $\rho_1=\rho$, Eq.
(112) reduces to
$$
\omega_k = \frac{k^2}{2m} + \frac{1}{V}\; \sum_p n_p(\Phi_{k+p}-
\Phi_p ) \; .
$$
The Fourier transform
$$
\Phi_{k+p} = \int \Phi({\bf r}) e^{-i({\bf k}+{\bf p})\cdot{\bf r}} \; d{\bf r}
$$
can be simplfied remembering that, by assumption, the interaction potential
diminishes fast with increasing $r\equiv|{\bf r}|$. Then in the above integral,
one can expand $e^{-i{\bf k}\cdot{\bf r}}$ in powers of ${\bf k}\cdot{\bf r}$ up to the
second order. As a result, we get
$$
\Phi_{k+p} \cong \left ( 1 \; - \; \frac{1}{6}\; k^2 r_0^2 \right )
\Phi_p \; ,
$$
where the notation for the effective interaction radius $r_0$ is introduced,
defined by the equation
$$
r_0^2 \equiv \frac{\int r^2\Phi({\bf r})d{\bf r}}{\int\Phi({\bf r})d{\bf r}} \; .
$$
In this way, we find
$$
\omega_k = \frac{k^2}{2m} \; - \; \frac{1}{6V}\; \sum_p n_p
\Phi_p k^2 r_0^2 \; .
$$
The critical temperature is given by the equation
$$
\rho = \int n_k \; \frac{d{\bf k}}{(2\pi)^3} \; ,
$$
in which
$$
n_k =\left ( e^{\omega_k/T_c} - 1\right )^{-1} \; .
$$
From here we obtain
$$
T_c = \frac{2\pi}{m^*} \left [
\frac{\rho}{\zeta(3/2)} \right ]^{2/3} \; ,
$$
where the effective mass is
$$
m^* \equiv
\frac{m}{1-\frac{mr_0^2}{3}\int n_k\Phi_k\frac{d{\bf k}}{(2\pi)^3}} \; .
$$
Thus, for nonlocal interactions, with an interaction radius $r_0$, the
effective mass increases, so that the critical temperature diminishes, as
compared to the critical temperature of the ideal Bose gas.

However, for the contact interaction potential (83), we have $r_0=0$, hence
$m^*=m$, and $T_c$, in the frame of the HFB mean-field approximation,
coincides with the ideal gas condensation temperature.

This conclusion is in agreement with other studies of the critical
temperature for interacting Bose gas. There exists quite a number of such
investigations, as reviewed in Refs. [31,49]. The most accurate of these
calculations are those employing Monte Carlo simulations [67--71] and those
based on the optimized perturbation theory [72--74], as has been done in
Refs. [75--80]. These investigations show that the first correction to the
critical temperature comes from effects beyond the mean-field approximation.

\newpage

redpajama_set_name: RedPajamaArXiv