arXiv:0704.1089v1  [cond-mat.stat-mech]  9 Apr 2007
Representative Ensembles in Statistical Mechanics
V.I. Yukalov
Bogolubov Laboratory of Theoretical Physics,
Joint Institute for Nuclear Research, Dubna 141980, Russia
Abstract
The notion of representative statistical ensembles, correctly representing statistical
systems, is strictly formulated. This notion allows for a proper description of statis-
tical systems, avoiding inconsistencies in theory. As an illustration, a Bose-condensed
system is considered. It is shown that a self-consistent treatment of the latter, using a
representative ensemble, always yields a conserving and gapless theory.
Key words: Statistical systems; representative statistical ensembles; Bose-condensed
systems.
PACS numbers: 05.30.Ch, 05.30.Jp, 05.70.Ce, 64.10.+h, 67.40.Db
http://arxiv.org/abs/0704.1089v1
1 Introduction
Statistical systems are characterized by statistical ensembles. It is crucially important that
the given statistical system be correctly represented by a statistical ensemble. In other words,
the chosen statistical ensemble must be representative for the considered statistical system.
This necessitates a thorough definition of what, actually, a statistical system is, requiring an
accurate enumeration of all its basic features. The usage of a nonrepresentative ensemble,
incorrectly representing the considered statistical system, may lead, and often does lead, to
inconsistencies in the theoretical description of the system.
The necessity of defining a statistical ensemble that would correctly represent the given
statistical system was, first, emphasized already by Gibbs [1], who stressed that all additional
conditions and constraints, imposed on the system, must be taken into account. The prob-
lem of a proper representation of equilibrium statistical systems by equilibrium statistical
ensembles was discussed by ter Haar [2,3] and also analized in the review article [4].
The aim of the present paper is to formalize the notion of representative statistical en-
sembles by giving precise mathematical definitions and to generalize this notion for arbitrary
systems, whether equilibrium or nonequilibrium. The application of the notion is illustrated
by systems with Bose-Einstein condensate, when the global gauge symmetry is broken. It is
shown that employing a representative ensemble for a Bose-condensed system results in the
theory enjoying conservation laws and having no gap in the spectrum of collective excitations.
Throughout the paper, the system of units will be used with the Planck and Boltzmann
constants set to unity, h̄ = 1, kB = 1.
2 Representative Ensembles
Let us, first, recall several general preliminary definitions that are necessary for precisely
defining the basic notion of a representative statistical ensemble.
Physical system is a collection of objects characterized by their typical features distin-
guishing this collection from other systems.
For example, a collection of particles can be characterized by their Hamiltonian, that is,
by their energy operator.
Statistical system is a many-body physical system, whose typical features are compli-
mented by all additional constraints and conditions which are necessary for uniquely de-
scribing the statistical properties of the system.
Statistical systems are characterized by statistical ensembles.
Statistical ensemble is a pair {F , ρ̂(t)} composed by the space of microstates F and a
statistical operator ρ̂(t) on that space.
The space of microstates can be the Fock space or its appropriate subspace. The statis-
tical operator ρ̂(t), generally, is a function of time t. To give ρ̂(t) implies to define its form
ρ̂(0) at the initial time t = 0 and to specify the evolution operator Û(t) such that
ρ̂(t) = Û(t)ρ̂(0)Û+(t) . (1)
Therefore, a statistical ensemble can be defined as a triplet
{F , ρ̂(0), Û(t)} ←→ {F , ρ̂(t)} .
The knowledge of a statistical ensemble allows one to find statistical averages.
Statistical average for an operator Â(t) on F is
< Â(t) > = TrF ρ̂(t)Â(0) = TrF ρ̂(0)Â(t) . (2)
Here the Heisenberg representation of the operator Â(t) is assumed, for which
Â(t) = Û+(t)Â(0)Û(t) . (3)
Representative ensemble is a statistical ensemble equipped with all additional constraints
and conditions that are necessary for a unique representation of the given statistical system.
Additional constraints and conditions for statistical systems are usually formulated as
conditions on statistical averages for some specified condition operators Ĉi(t), where i =
1, 2, . . .. These operators do not need to be necessarily the integrals of motion, but they are
supposed to be Hermitian.
Statistical condition is a prescribed equality for the statistical average of a condition
operator,
Ci(t) = < Ĉi(t) > = Trρ̂(0)Ĉi(t) . (4)
Here and in what follows, the trace operation is assumed to be over the appropriate space
of microstates F .
Let us consider, first, an equilibrium statistical system, for which the statistical operator
does not depend on time,
ρ̂(t) = ρ̂(0) ≡ ρ̂ . (5)
The explicit form of the statistical operator follows from the principle of minimal information
[5]. The latter presumes the conditional maximization of the Gibbs entropy
S = −Trρ̂ ln ρ̂ (6)
under the statistical conditions (4), among which one usually distinguishes the definition of
the internal energy
E = Trρ̂Ĥ (7)
and the normalization condition
Trρ̂ = 1 . (8)
The information functional is
I[ρ̂] = −S + λ0 (Trρ̂− 1) + β
Trρ̂Ĥ −E
Trρ̂Ĉi − Ci
, (9)
where λ0 ≡ lnZ − 1 is the Lagrange multiplier preserving the normalization condition (8),
β is the inverse temperature, which is the Lagrange multiplier for condition (7), and βνi are
the Lagrange multipliers related to statistical conditions (4).
The minimization of the information functional (9) yields the statistical operator
e−βH , (10)
corresponding to the grand canonical ensemble with the grand Hamiltonian
H ≡ Ĥ +
νiĈi . (11)
The most customary expression for the grand Hamiltonian (11) is
H = Ĥ − µN̂ ,
where µ is the chemical potential and N̂ is the number-of-particle operator. However, the
general form of the grand Hamiltonian is given by Eq. (11), in which any condition operators
can be involved. Thus, an equilibrium representative ensemble is described by the statistical
operator (10) with the grand Hamiltonian (11). The evolution operator for an equilibrium
system is
Û(t) = e−iHt , (12)
which commutes with the statistical operator (10), because of which
ρ̂(t) = [H, ρ̂(t)] = 0 .
The general way of obtaining the evolution equations for arbitrary nonequilibrium sys-
tems is through the extremization of action functionals [6]. In our case, this extremization
has to be accomplished under the prescribed statistical conditions (4).
Let the system Hamiltonian be a functional of the field operators ψ(x, t) and ψ†(x, t),
that is, Ĥ = Ĥ [ψ], where ψ = ψ(x, t). The system Lagrangian is
L̂[ψ] ≡
ψ†(x, t) i
ψ(x, t) dx− Ĥ[ψ] . (13)
The action functional, or effective action, under the prescribed statistical conditions (4),
takes the form
A[ψ] ≡
L̂[ψ]−
νiĈi(t)
dt , (14)
where νi are the Lagrange multipliers guaranteeing the validity of the given statistical con-
ditions. The action functional is defined so that to be a self-adjoint operator,
A+[ψ] = A[ψ] . (15)
Similarly to Eq. (11), the grand Hamiltonian in the Heisenberg representation is
H [ψ] = Ĥ [ψ] +
νiĈi(t) . (16)
Then the effective action (14) can be rewritten as
A[ψ] =
ψ†(x, t) i
ψ(x, t) dx−H [ψ]
dt . (17)
The extremization of the action functional, requiring that
δA[ψ] = 0 ,
δA[ψ] =
δA[ψ]
δψ(x, t)
δψ(x, t) +
δA[ψ]
δψ†(x, t)
δψ†(x, t) ,
yields the evolution equations
δA[ψ]
δψ†(x, t)
= 0 ,
δA[ψ]
δψ(x, t)
= 0 . (18)
These equations are the Hermitian conjugated forms of each other.
From Eqs. (17) and (18), it is evident that the evolution equations for the field operators
can be represented as
ψ(x, t) =
δH [ψ]
δψ†(x, t)
and its Hermitian conjugated. This should be equivalent to the Heisenberg equation of
motion
ψ(x, t) = [ψ(x, t), H [ψ]] ,
that is, to the Heisenberg representation for the field operator
ψ(x, t) = Û+(t)ψ(x, 0)Û(t) .
Hence, the evolution operator satisfies the Schrödinger equation
Û(t) == H [ψ(x, 0)] Û(t) . (20)
In this way, a nonequilibrium representative ensemble is the set of the given space of
microstates F , initial statistical operator ρ̂(0), and of the evolution operator Û(t) defined
by Eq. (20). An equilibrium representative ensemble is, of course, just a particular case of
the general nonequilibrium ensemble.
3 Bose-Condensed Systems
To illustrate the explicit construction of a representative ensemble, let us consider a system
with Bose-Einstein condensate. Such systems possess a variety of interesting properties, as
can be inferred from review works [7–10]. Moreover, theoretical description of these systems
is known to confront the notorious difficulty of defining a self-consistent approach. The
theory of Bose-condensed systems is based on the Bogolubov idea [11-14] of breaking the
global gauge symmetry by means of the famous Bogolubov shift for field operators. The
condensate wave function, introduced in the course of this shift, has to satisfy the minimum
of the related thermodynamic potential, which is the stability condition necessary for mak-
ing the system stable and the theory conserving and self-consistent. At the same time, the
spectrum of collective excitations, according to the Hugenholtz-Pines theorem [15], has to
be gapless. The notorious problem is the appearance of the contradiction between the above
two requirements, when the theory is either nonconserving or gapful. This contradiction
does not arise only in the lowest orders with respect to particle interactions, when one uses
the Bogolubov approximation at low temperatures [11,12] or the quasiclassical approxima-
tion at high temperatures [16]. However, this contradiction immediately arises as soon as
the interaction strength is not asymptotically weak and one has to invoke a more elaborate
approximation. This problem of conserving versus gapless approximations was first empha-
sized by Hohenberg and Martin [17] and recently covered comprehensively by Andersen [9].
The problem is caused by the usage of nonrepresentative ensembles, which renders the sys-
tem unstable [18]. Here we show that employing a representative ensemble never yields the
above contradiction, always resulting in a self-consistent theory, being both conserving and
gapless.
We consider a system with the Hamiltonian
ψ†(r)
ψ(r) dr+
ψ†(r)ψ†(r′)Φ(r− r′)ψ(r′)ψ(r) drdr′ , (21)
in which the field operators ψ(r) = ψ(r, t) satisfy the Bose commutation relations, U =
U(r, t) is an external field, and Φ(r) = Φ(−r) is an interaction potential. For describing a
Bose-condensed system with broken global gauge symmetry, the Bogolubov shift [13,14] has
to be done through the replacement
ψ(r, t) −→ ψ̂(r, t) ≡ η(r, t) + ψ1(r, t) , (22)
where η(r, t) is the condensate wave function and ψ1(r, t) is the field operator of noncon-
densed particles. The latter field variables are assumed to be orthogonal to each other,
η∗(r, t)ψ1(r, t) dr = 0 . (23)
It is necessary to emphasize that the Bogolubov shift (22) realizes unitary nonequivalent
operator representations [18,19]. Accomplishing the Bogolubov shift (22) in Hamiltonian
(21), as well as in all operators of observables, we get the algebra of observables defined on
the Fock space F(ψ1) generated by the field operators ψ†1(r) (see details in Refs. [5,18,19]).
The condensate function is normalized to the number of condensed particles
|η(r, t)|2dr . (24)
The Bogolubov shift (22) is only rational when the number of condensed particles (24) is
macroscopic, which means that the limit
is not zero, where N is the total number of particles. The latter is given by the average
N = < N̂ > (25)
for the number-of-particle operator
ψ̂†(r)ψ̂(r) dr , (26)
in which the Bogolubov shift (22) is again assumed. The statistical averaging in Eq. (25)
and everywhere below is over the Fock space F(ψ1).
Substituting the Bogolubov shift (22) into Hamiltonian (21) gives in the latter the terms
linear in ψ1, because of which the average < ψ1 > can be nonzero. This, however, would
result in the nonconservation of quantum numbers, e.g., of spin or momentum. Therefore,
one has to impose the constraint for the conservation of quantum numbers,
< ψ1(r, t) > = 0 . (27)
Defining the self-adjoint condition operator
Λ̂(t) ≡
λ(r, t)ψ
1(r, t) + λ
∗(r, t)ψ1(r, t)
dr , (28)
in which λ(r, t) is a complex function, we may represent constraint (27) as the quantum
conservation condition
< Λ̂(t) > = 0 . (29)
In this way, there are three statistical conditions. The first condition is the normalization
(24) for the number of condensed particles. Condition (24) can be represented in the standard
form (4) by defining the operator
N̂0 ≡ 1̂
|η(r, t)|2dr , (30)
in which 1̂ is the unity operator in the Fock space F(ψ1). Then Eq. (24) reduces to the
statistical condition
N0 = < N̂0 > . (31)
The second condition is the normalization (25) for the total number of particles. Equivalently,
instead of normalization (25), we may consider the normalization condition
N1 = < N̂1 > , N̂1 ≡
1(r)ψ1(r) dr (32)
for the number of uncondensed particles N1 = N − N0. And the third condition is the
conservation condition (29). Respectively, the effective action, which is now a functional of
the two field variables, η(r, t) and ψ1(r, t), with taking account of the statistical conditions
(29), (31), and (32), becomes
A[η, ψ1] =
L̂+ µ0N̂0 + µ1N̂1 + Λ̂
dt . (33)
Here L̂ = L̂[ψ̂] is the Lagrangian (13) under the Bogolubov shift (22) and Λ̂ = Λ̂(t) from
Eq. (28). The quantities µ0, µ1, and λ(r, t) are the Lagrange multipliers guaranteeing the
validity of the corresponding statistical conditions. Introducing the grand Hamiltonian
H [η, ψ1] ≡ Ĥ − µ0N̂0 − µ1N̂1 − Λ̂ , (34)
in which Ĥ = Ĥ[ψ̂], with shift (22), and the effective Lagrangian
L[η, ψ1] ≡
η∗(r, t) i
η(r, t) + ψ
1(r, t) i
ψ1(r, t)
dr−H [η, ψ1] , (35)
for the action functional (33), we get
A[η, ψ1] =
L[η, ψ1] dt . (36)
The evolution equations follow from the extremization of the action functional (36), that
is, from the variations
δA[η, ψ1]
δη∗(r, t)
= 0 (37)
δA[η, ψ1]
1(r, t)
= 0 . (38)
These equations, as is clear from Eqs. (35) and (36), are equivalent to the equations of
motion
η(r, t) =
δH [η, ψ1]
δη∗(r, t)
ψ1(r, t) =
δH [η, ψ1]
1(r, t)
. (40)
One has to substitute here Hamiltonian (34) under the Bogolubov shift (22). Accomplishing
the variation in Eq. (39), we get
η(r, t) =
+ U − µ0
η(r)+
Φ(r− r′)
|η(r′)|2η(r) + X̂(r, r′)
dr′ , (41)
where the time dependence in the right-hand side, for short, is not explicitly shown, U =
U(r, t), and the notation for the correlation operator
X̂(r, r′) ≡ ψ†1(r′)ψ1(r′)η(r) + ψ
′)η(r′)ψ1(r)+
+η∗(r′)ψ1(r
′)ψ1(r) + ψ
′)ψ1(r
′)ψ1(r) (42)
is used. The variation in Eq. (40) gives
ψ1(r, t) =
+ U − µ1
ψ1(r)+
Φ(r− r′)
|η(r′)|2ψ1(r) + η∗(r′)η(r)ψ1(r′) + η(r′)η(r)ψ†1(r′) + X̂(r, r′)
dr′ . (43)
To get an equation for the condensate wave function, we have to take the statistical
average of Eq. (41). For this purpose, we introduce the normal density matrix
ρ1(r, r
′) ≡ < ψ†1(r′)ψ1(r) > , (44)
the anomalous density matrix
σ1(r, r
′) ≡ < ψ1(r′)ψ1(r) > , (45)
and their diagonal elements, giving the density of noncondensed particles
ρ1(r) ≡ ρ1(r, r) = < ψ†1(r)ψ1(r) > (46)
and the anomalous average
σ1(r) ≡ σ1(r, r) = < ψ1(r)ψ1(r) > . (47)
The quantity |σ1(r)| can be interpreted as the density of paired particles [19]. The total
density of particles
ρ(r) = ρ0(r) + ρ1(r) (48)
consists of the condensate density
ρ0(r) ≡ |η(r)|2 (49)
and the density of noncondensed particles (46). Averaging Eq. (41), we find the equation
for the condensate wave function
η(r, t) =
+ U − µ0
η(r)+
Φ(r− r′)
ρ(r′)η(r) + ρ1(r, r
′)η(r′) + σ1(r, r
′)η∗(r′)+ < ψ
′)ψ1(r
′)ψ1(r) >
dr′ . (50)
Equations (43) and (50) are the basic equations of motion for the field variables η(r, t) and
ψ1(r, t). These equations, according to Eqs. (39) and (40), are generated by the variation
of the grand Hamiltonian (34). The latter, in agreement with Eq. (20), also defines the
evolution operator Û(t), which satisfies the Schrödinger equation
Û(t) = H [η(r, 0), ψ1(r, 0)]Û(t) .
Thus, the representative ensemble for a Bose-condensed system is the triplet
{F(ψ1), ρ̂(0), Û(t)} .
It is important to stress that the so defined representative ensemble possesses a principal
feature making it different from the standardly used ensemble having the sole Lagrange
multiplier µ0 ≡ µ1. But then the normalization condition (24) cannot be guaranteed. Then
the evolution equation for the condensate wave function is not a result of a variational
procedure. For an equilibrium system, this means that the number of condensed particles
N0 does not provide the minimum of a thermodynamic potential, which implies the system
instability. All notorious inconsistencies in theory, manifesting themselves in the lack od
conservation laws or in the appearance of an unphysical gap in the spectrum, are caused by
the usage of nonrepresentative ensembles.
4 Green Functions
The equations of motion (43) and (50) allow us to derive the evolutional equations for the
Green functions. To this end, we shall use the compact notation denoting the set {rj, tj} by
the sole letter j, so that the dependence of functions on the spatial and temporal variables
looks like
f(12 . . . n) ≡ f(r1, t1, r2, t2, . . . , rn, tn) .
The product of the differentials drjdtj will be denoted as d(j), so that
d(12 . . . n) ≡
drj dtj .
We shall employ the Dirac delta function
δ(12) ≡ δ(r1 − r2) δ(t1 − t2) .
For the interaction potential, we shall use the retarded form
Φ(12) ≡ Φ(r1 − r2)δ(t1 − t2 + 0) . (51)
The matrix Green function G(12) = [Gαβ(12)] is a 2 × 2 matrix, with α, β = 1, 2, and
with the following elements:
G11(12) ≡ −i < T̂ψ1(1)ψ†1(2) > , G12(12) ≡ −i < T̂ψ1(1)ψ1(2) > ,
G21(12) ≡ −i < T̂ψ†1(1)ψ
1(2) > , G22(12) ≡ −i < T̂ψ
1(1)ψ1(2) > , (52)
where T̂ is the time-ordering operator.
Let us introduce the operator
K̂j ≡ −
+ U(j)− µ1 , (53)
the condensate effective potential
V (12) ≡ δ(12)
Φ(13)|η(3)|2d(3) + Φ(12)η(1)η∗(2) , (54)
and let us rewrite the correlation operator (42) in the form
X̂(12) = ψ
1(2)ψ1(2)η(1) + ψ
1(2)η(2)ψ1(1) + η
∗(2)ψ1(2)ψ1(1) + ψ
1(2)ψ1(2)ψ1(1) . (55)
We also define the matrix correlation function X(123) = [Xαβ(123)] with the elements:
X11(123) ≡ − < T̂ X̂(12)ψ†1(3) > , X12(123) ≡ − < T̂ X̂(12)ψ1(3) > ,
X21(123) ≡ − < T̂X̂+(12)ψ†1(3) > , X22(123) ≡ − < T̂X̂+(12)ψ1(3) > . (56)
From the equations of motion (43) and (50), we find the equations
− K̂1
G11(12)−
V (13)G11(32) d(3)−
Φ(13) [η(1)η(3)G21(32) + iX11(132)] d(3) = δ(12) ,
− K̂1
G12(12)−
V (13)G12(32) d(3)−
Φ(13) [η(1)η(3)G22(32) + iX12(132)] d(3) = 0 ,
− K̂1
G21(12)−
V ∗(13)G21(32) d(3)−
Φ(13) [η∗(1)η∗(3)G11(32) + iX21(132)] d(3) = 0 ,
− K̂1
G22(12)−
V ∗(13)G22(32) d(3)−
Φ(13) [η∗(1)η∗(3)G12(32) + iX22(132)] d(3) = δ(12) . (57)
The self-energy Σ(12) = [Σαβ(12)] is a matrix whose elements are defined by the relations
[Σ11(13)G11(32) + Σ12(13)G21(32)] d(3) =
V (13)G11(32) d(3) +
Φ(13) [η(1)η(3)G21(32) + iX11(132)] d(3) ,
[Σ11(13)G12(32) + Σ12(13)G22(32)] d(3) =
V (13)G12(32) d(3) +
Φ(13) [η(1)η(3)G22(32) + iX12(132)] d(3) ,
[Σ21(13)G11(32) + Σ22(13)G21(32)] d(3) =
V ∗(13)G21(32) d(3) +
Φ(13) [η∗(1)η∗(3)G11(32) + iX21(132)] d(3) ,
[Σ21(13)G12(32) + Σ22(13)G22(32)] d(3) =
V ∗(13)G22(32) d(3) +
Φ(13) [η∗(1)η∗(3)G12(32) + iX22(132)] d(3) . (58)
Let us introduce the matrix condensate propagator C(12) = [Cαβ(12)], with the elements
C11(12) ≡ −iη(1)η∗(2) , C12(12) ≡ −iη(1)η(2) ,
C21(12) ≡ −iη∗(1)η∗(2) , C22(12) ≡ −iη∗(1)η(2) , (59)
The latter have the properties
C11(21) = C22(12) , C11(11) = C22(11) , C12(21) = C12(12) ,
C21(21) = C21(12) , C
11(12) = −C22(12) , C∗12(12) = −C21(12) . (60)
The binary Green function is a matrix B(123) = [Bαβ(123)],
B(123) = C(12)G(23)− iρ0(2)G(13) +X(123) , (61)
whose elements are
B11(123) ≡ C11(12)G11(23) + C11(22)G11(13) + C12(12)G21(23) +X11(123) ,
B12(123) ≡ C11(12)G12(23) + C11(22)G12(13) + C12(12)G22(23) +X12(123) ,
B21(123) ≡ C22(12)G21(23) + C22(22)G21(13) + C21(12)G11(23) +X21(123) ,
B22(123) ≡ C22(12)G22(23) + C22(22)G22(13) + C21(12)G12(23) +X22(123) , (62)
and where ρ0(1) ≡ |η(1)|2.
With Eqs. (59) and (61), relations (58), defining the self-energy, can be rewritten in the
matrix form
Σ(13)G(32) d(3) = i
Φ(13)B(132) d(3) . (63)
Then the equations of motion (57) acquire the matrix representation
τ̂3 i
− K̂1
G(12)−
Σ(13)G(32) d(3) = δ(12) , (64)
in which the delta function δ(12) in the right-hand side is assumed to be factored with the
unity matrix 1̂ = [δαβ ] and
τ̂3 ≡
is a Pauli matrix.
Introducing the inverse propagator
G−1(12) ≡
τ̂3 i
− K̂1
δ(12)− Σ(12) (65)
allows us to transform Eq. (64) into
G−1(13)G(32) d(3) = δ(12) . (66)
An equivalent representation, following from Eq. (66), is
G(13)G−1(32) d(3) = δ(12) . (67)
For the self-energy, using Eq. (63), we have
Σ(12) = i
Φ(13)B(134)G−1(42) d(34) . (68)
The equations for the Green functions are to be complimented by the equation for the
condensate wave function (50), which, introducing one more anomalous average
ξ(12) ≡ < ψ†1(2)ψ1(2)ψ(1) > , (69)
can be represented as
η(1) =
+ U(1)− µ0
η(1)+
Φ(12) [ρ(2)η(1) + ρ1(12)η(2) + σ1(12)η
∗(2) + ξ(12)] d(2) . (70)
It is the equations for the Green functions and the equation for the condensate func-
tion, which become mutually incompatible in the standard approach, while employing the
representative ensemble renders the theory self-consistent in any approximation.
5 Theory Self-Consistency
One usually confronts inconsistencies in theory considering a uniform equilibrium Bose-
condensed system. Then, in any given approximation, one gets either a nonconserving
theory, that is, an unstable system, or one finds an unphysical gap in the spectrum, which,
actually, again corresponds to an unstable system [9,18]. To analyze this problem, we pass
now to the case of an equilibrium uniform system, when U = 0.
Then we use the Fourier transform for the Green function
G(12) =
G(k, ω)ei(k·r12−ωt12)
(2π)4
in which
r12 ≡ r1 − r2 , t12 ≡ t1 − t2 .
By their definition in Eq. (52), the Green function elements possess the properties
G11(21) = G22(12) , G12(21) = G12(12) , G21(21) = G21(12) . (71)
Therefore the corresponding Fourier transforms satisfy the relations
G11(−k,−ω) = G22(k, ω) , G12(−k,−ω) = G12(k, ω) ,
G21(−k,−ω) = G21(k, ω) . (72)
Assuming that the system is isotropic, one has
Gαβ(−k, ω) = Gαβ(k, ω) (73)
for all α, β. Combining Eqs. (72) and (73), we find
G11(k,−ω) = G22(k, ω) , G12(k,−ω) = G12(k, ω) ,
G21(k,−ω) = G21(k, ω) . (74)
Also, for a uniform equilibrium system, one has [14] the equality
G21(k, ω) = G12(k, ω) . (75)
Fourier-transforming the self-energy
Σ(12) =
Σ(k, ω)ei(k·r12−ωt12)
(2π)4
and, similarly, the inverse propagator (65), we have for the latter
G−1(k, ω) = τ̂3ω −
+ µ− Σ(k, ω) . (76)
Then Eq. (67) reduces to
G−1(k, ω)G(k, ω) = 1̂ , (77)
where 1̂ = [δαβ ].
From Eqs. (76) and (77), it follows that G−1(k, ω), hense, also Σ(k, ω), have the same
symmetry properties as G(k, ω). In particular,
Σαβ(−k, ω) = Σαβ(k, ω) , Σ12(k,−ω) = Σ12(k, ω) , Σ21(k,−ω) = Σ21(k, ω) ,
Σ21(k, ω) = Σ12(k, ω) , Σ11(k,−ω) = Σ22(k, ω) . (78)
The matrix equation (77), explicitly, is the system of equations
ω − k
+ µ1 − Σ11
G11 − Σ12G21 = 1 ,
+ µ1 − Σ11
G12 − Σ12G22 = 0 ,
−ω − k
+ µ1 − Σ22
G21 − Σ21G11 = 0 ,
−ω − k
+ µ1 − Σ22
G22 − Σ21G12 = 1 , (79)
where, for short, Gαβ = Gαβ(k, ω) and Σαβ = Σαβ(k, ω). The solutions to these equations
G11(k, ω) =
ω + k2/2m+ Σ11(k, ω)− µ1
D(k, ω)
G12(k, ω) = −
Σ12(k, ω)
D(k, ω)
, (80)
with the denominator
D(k, ω) ≡
ω − k
− Σ11(k, ω) + µ1
+ Σ22(k, ω)− µ1
+ Σ212(k, ω) . (81)
The solutions for Σ21(k, ω) and G22(k, ω) are defined by the symmetry properties (72) to
(75).
The excitation spectrum is given by the poles of the Green functions, that is, by the zero
of denominator (81),
D(k, εk) = 0 . (82)
Equation (82) can be represented as
[Σ11(k, εk)− Σ22(k, εk)]±
ω2k − Σ212(k, εk) , (83)
with the notation
[Σ11(k, εk) + Σ22(k, εk)]− µ1 . (84)
Denominator (81) enjoys the property
D(k,−ω) = D(k, ω) .
Consequently, if εk is a solution of Eq. (82), then −εk is also its solution, which is in
agreement with the form of Eq. (83).
For an equilibrium uniform system, the Bogolubov shift (22) is equivalent to the sepa-
ration of the zero-momentum term in the expansion of the field operator over plane waves.
The shift itself has meaning only under the normalization condition (24), in which N0 ∼ N ,
that is, the zero-momentum state is macroscopically occupied. The latter becomes possible
when the single particle spectrum touches zero. Therefore, the necessary condition for the
existence of Bose-Einstein condensate is
εk = 0 . (85)
This is to be complimented by the stability condition
Re εk ≥ 0 , Im εk ≤ 0 . (86)
This condition should be kept in mind when choosing the sign plus in front of the square
root in spectrum (83).
Taking limit (85) for spectrum (83), we notice that, according to properties (78),
Σ11(k, 0) = Σ22(k, 0) . (87)
By using perturbation theory for a stable system, one can show [15] that in all orders of the
theory
Σαβ(0, 0) ≥ 0 . (88)
Then the necessary condition (85) yields the expression for the chemical potential
µ1 = Σ11(0, 0)− Σ12(0, 0) , (89)
which is the Hugenholtz-Pines relation [15].
On the other hand, we have Eq. (70) for the condensate wave function. For an equilibrium
uniform system, with no external potential U , all densities do not depend on the spatial and
temporal variables,
ρ0(r) = ρ0 , ρ1(r) = ρ1 , σ1(r) = σ1 , ρ(r) = ρ . (90)
The condensate wave function reduces to the constant
η(r, t) = η =
ρ0 . (91)
Then we substitute into Eq. (70) the Fourier transforms for the interaction potential
Φ(r) =
ik·r dk
(2π)3
for the normal density matrix (44),
ρ1(r1, r2) =
ik·r12
(2π)3
and for the anomalous density matrix (45),
σ1(r1, r2) =
ik·r12
(2π)3
Similarly, the Fourier transform for the anomalous average (69) is
ξ1(r1, r2) =
ik·r12
(2π)3
As a result, Eq. (70) gives
µ0 = ρΦ0 +
nk + σk +
(2π)3
. (92)
Generally, expressions (92) and (89) do not coincide with each other, their difference
being
µ0 − µ1 = ρΦ0 +
nk + σk +
(2π)3
−Σ11(0, 0) + Σ12(0, 0) . (93)
This is the general expression for the difference between the Lagrange multipliers µ0 and µ1
for an arbitrary equilibrium uniform Bose-condensed system.
Usually, one does not distinguish between the Lagrange multipliers µ0 and µ1, which
implies setting µ0 − µ1 → 0. However, as is evident from Eq. (93), there is no any reason
for requiring that this quantity be zero. As an illustration, we may resort to the Hartree-
Fock-Bogolubov approximation, in which ξk = 0 and
Σ11(0, 0) = (ρ+ ρ0)Φ0 +
(2π)3
, Σ12(0, 0) = ρ0Φ0 +
(2π)3
Relation (89) then yields
µ1 = ρΦ0 +
(nk − σk)Φk
(2π)3
. (94)
The difference of the chemical potentials (93) becomes
µ0 − µ1 = 2
(2π)3
, (95)
which, certainly, is nonzero [20].
In this way, the introduction of the additional Lagrange multiplier makes the theory
completely self-consistent. All inconsistencies that often arise in other works, such as the
appearance of a gap in the spectrum, system instability or a distortion of the phase transition
order, are caused by neglecting the difference between the multiplier µ1 and the multiplier
µ0. It is worth emphasizing that the introduction of the Lagrange multiplier µ0 for preserving
the normalization condition (24), from the mathematical point of view is strictly necessary.
In other case, the employed ensemble would not be representative, hence, could not correctly
describe the Bose-condensed system with broken gauge symmetry.
Acknowledgement
I am grateful for many useful discussions to M. Girardeau, R. Graham, H. Kleinert, and
E.P. Yukalova.
References
[1] J.W. Gibbs, Collected Works (Longmans, New York, 1931), Vol. 2.
[2] D. ter Haar, Elements of Statistical Mechanics (Reinhart, New York, 1954).
[3] D. ter Haar, Rep. Prog. Phys. 24, 304 (1961).
[4] V.I. Yukalov, Phys. Rep. 208, 395 (1991).
[5] V.I. Yukalov, Statistical Green’s Functions (Queen’s University, Kingston, 1998).
[6] H. Kleinert, Path Integrals (World Scientific, Singapore, 2004).
[7] P.W. Courteille, V.S. Bagnato and V.I. Yukalov, Laser Phys. 11, 659 (2001).
[8] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation (Clarendon, Oxford, 2003).
[9] J.O. Andersen, Rev. Mod. Phys. 76, 599 (2004).
[10] K. Bongs and K. Sengstock, Rep. Prog. Phys. 67, 907 (2004).
[11] N.N. Bogolubov, J. Phys. (Moscow) 11, 23 (1947).
[12] N.N. Bogolubov, Moscow Univ. Phys. Bull. 7, 43 (1947).
[13] N.N. Bogolubov, Lectures on Quantum Statistics (Gordon and Breach, New York, 1967),
Vol. 1.
[14] N.N. Bogolubov, Lectures on Quantum Statistics (Gordon and Breach, New York, 1970),
Vol. 2.
[15] N.M. Hugenholtz and D. Pines, Phys. Rev. 116, 489 (1959).
[16] N. Prokofiev, O. Ruebenacker and B. Svistunov, Phys. Rev. A 69, 053625 (2004).
[17] P.C. Hohenberg and P.C. Martin, Ann. Phys. 34, 291 (1965).
[18] V.I. Yukalov, Phys. Rev. E 72, 066119 (2005).
[19] V.I. Yukalov, Laser Phys. 16, 511 (2006).
[20] V.I. Yukalov and E.P. Yukalova, Laser Phys. Lett. 2, 506 (2005).

  The notion of representative statistical ensembles, correctly representing
statistical systems, is strictly formulated. This notion allows for a proper
description of statistical systems, avoiding inconsistencies in theory. As an
illustration, a Bose-condensed system is considered. It is shown that a
self-consistent treatment of the latter, using a representative ensemble,
always yields a conserving and gapless theory.

Introduction
Statistical systems are characterized by statistical ensembles. It is crucially important that
the given statistical system be correctly represented by a statistical ensemble. In other words,
the chosen statistical ensemble must be representative for the considered statistical system.
This necessitates a thorough definition of what, actually, a statistical system is, requiring an
accurate enumeration of all its basic features. The usage of a nonrepresentative ensemble,
incorrectly representing the considered statistical system, may lead, and often does lead, to
inconsistencies in the theoretical description of the system.
The necessity of defining a statistical ensemble that would correctly represent the given
statistical system was, first, emphasized already by Gibbs [1], who stressed that all additional
conditions and constraints, imposed on the system, must be taken into account. The prob-
lem of a proper representation of equilibrium statistical systems by equilibrium statistical
ensembles was discussed by ter Haar [2,3] and also analized in the review article [4].
The aim of the present paper is to formalize the notion of representative statistical en-
sembles by giving precise mathematical definitions and to generalize this notion for arbitrary
systems, whether equilibrium or nonequilibrium. The application of the notion is illustrated
by systems with Bose-Einstein condensate, when the global gauge symmetry is broken. It is
shown that employing a representative ensemble for a Bose-condensed system results in the
theory enjoying conservation laws and having no gap in the spectrum of collective excitations.
Throughout the paper, the system of units will be used with the Planck and Boltzmann
constants set to unity, h̄ = 1, kB = 1.
2 Representative Ensembles
Let us, first, recall several general preliminary definitions that are necessary for precisely
defining the basic notion of a representative statistical ensemble.
Physical system is a collection of objects characterized by their typical features distin-
guishing this collection from other systems.
For example, a collection of particles can be characterized by their Hamiltonian, that is,
by their energy operator.
Statistical system is a many-body physical system, whose typical features are compli-
mented by all additional constraints and conditions which are necessary for uniquely de-
scribing the statistical properties of the system.
Statistical systems are characterized by statistical ensembles.
Statistical ensemble is a pair {F , ρ̂(t)} composed by the space of microstates F and a
statistical operator ρ̂(t) on that space.
The space of microstates can be the Fock space or its appropriate subspace. The statis-
tical operator ρ̂(t), generally, is a function of time t. To give ρ̂(t) implies to define its form
ρ̂(0) at the initial time t = 0 and to specify the evolution operator Û(t) such that
ρ̂(t) = Û(t)ρ̂(0)Û+(t) . (1)
Therefore, a statistical ensemble can be defined as a triplet
{F , ρ̂(0), Û(t)} ←→ {F , ρ̂(t)} .
The knowledge of a statistical ensemble allows one to find statistical averages.
Statistical average for an operator Â(t) on F is
< Â(t) > = TrF ρ̂(t)Â(0) = TrF ρ̂(0)Â(t) . (2)
Here the Heisenberg representation of the operator Â(t) is assumed, for which
Â(t) = Û+(t)Â(0)Û(t) . (3)
Representative ensemble is a statistical ensemble equipped with all additional constraints
and conditions that are necessary for a unique representation of the given statistical system.
Additional constraints and conditions for statistical systems are usually formulated as
conditions on statistical averages for some specified condition operators Ĉi(t), where i =
1, 2, . . .. These operators do not need to be necessarily the integrals of motion, but they are
supposed to be Hermitian.
Statistical condition is a prescribed equality for the statistical average of a condition
operator,
Ci(t) = < Ĉi(t) > = Trρ̂(0)Ĉi(t) . (4)
Here and in what follows, the trace operation is assumed to be over the appropriate space
of microstates F .
Let us consider, first, an equilibrium statistical system, for which the statistical operator
does not depend on time,
ρ̂(t) = ρ̂(0) ≡ ρ̂ . (5)
The explicit form of the statistical operator follows from the principle of minimal information
[5]. The latter presumes the conditional maximization of the Gibbs entropy
S = −Trρ̂ ln ρ̂ (6)
under the statistical conditions (4), among which one usually distinguishes the definition of
the internal energy
E = Trρ̂Ĥ (7)
and the normalization condition
Trρ̂ = 1 . (8)
The information functional is
I[ρ̂] = −S + λ0 (Trρ̂− 1) + β
Trρ̂Ĥ −E
Trρ̂Ĉi − Ci
, (9)
where λ0 ≡ lnZ − 1 is the Lagrange multiplier preserving the normalization condition (8),
β is the inverse temperature, which is the Lagrange multiplier for condition (7), and βνi are
the Lagrange multipliers related to statistical conditions (4).
The minimization of the information functional (9) yields the statistical operator
e−βH , (10)
corresponding to the grand canonical ensemble with the grand Hamiltonian
H ≡ Ĥ +
νiĈi . (11)
The most customary expression for the grand Hamiltonian (11) is
H = Ĥ − µN̂ ,
where µ is the chemical potential and N̂ is the number-of-particle operator. However, the
general form of the grand Hamiltonian is given by Eq. (11), in which any condition operators
can be involved. Thus, an equilibrium representative ensemble is described by the statistical
operator (10) with the grand Hamiltonian (11). The evolution operator for an equilibrium
system is
Û(t) = e−iHt , (12)
which commutes with the statistical operator (10), because of which
ρ̂(t) = [H, ρ̂(t)] = 0 .
The general way of obtaining the evolution equations for arbitrary nonequilibrium sys-
tems is through the extremization of action functionals [6]. In our case, this extremization
has to be accomplished under the prescribed statistical conditions (4).
Let the system Hamiltonian be a functional of the field operators ψ(x, t) and ψ†(x, t),
that is, Ĥ = Ĥ [ψ], where ψ = ψ(x, t). The system Lagrangian is
L̂[ψ] ≡
ψ†(x, t) i
ψ(x, t) dx− Ĥ[ψ] . (13)
The action functional, or effective action, under the prescribed statistical conditions (4),
takes the form
A[ψ] ≡
L̂[ψ]−
νiĈi(t)
dt , (14)
where νi are the Lagrange multipliers guaranteeing the validity of the given statistical con-
ditions. The action functional is defined so that to be a self-adjoint operator,
A+[ψ] = A[ψ] . (15)
Similarly to Eq. (11), the grand Hamiltonian in the Heisenberg representation is
H [ψ] = Ĥ [ψ] +
νiĈi(t) . (16)
Then the effective action (14) can be rewritten as
A[ψ] =
ψ†(x, t) i
ψ(x, t) dx−H [ψ]
dt . (17)
The extremization of the action functional, requiring that
δA[ψ] = 0 ,
δA[ψ] =
δA[ψ]
δψ(x, t)
δψ(x, t) +
δA[ψ]
δψ†(x, t)
δψ†(x, t) ,
yields the evolution equations
δA[ψ]
δψ†(x, t)
= 0 ,
δA[ψ]
δψ(x, t)
= 0 . (18)
These equations are the Hermitian conjugated forms of each other.
From Eqs. (17) and (18), it is evident that the evolution equations for the field operators
can be represented as
ψ(x, t) =
δH [ψ]
δψ†(x, t)
and its Hermitian conjugated. This should be equivalent to the Heisenberg equation of
motion
ψ(x, t) = [ψ(x, t), H [ψ]] ,
that is, to the Heisenberg representation for the field operator
ψ(x, t) = Û+(t)ψ(x, 0)Û(t) .
Hence, the evolution operator satisfies the Schrödinger equation
Û(t) == H [ψ(x, 0)] Û(t) . (20)
In this way, a nonequilibrium representative ensemble is the set of the given space of
microstates F , initial statistical operator ρ̂(0), and of the evolution operator Û(t) defined
by Eq. (20). An equilibrium representative ensemble is, of course, just a particular case of
the general nonequilibrium ensemble.
3 Bose-Condensed Systems
To illustrate the explicit construction of a representative ensemble, let us consider a system
with Bose-Einstein condensate. Such systems possess a variety of interesting properties, as
can be inferred from review works [7–10]. Moreover, theoretical description of these systems
is known to confront the notorious difficulty of defining a self-consistent approach. The
theory of Bose-condensed systems is based on the Bogolubov idea [11-14] of breaking the
global gauge symmetry by means of the famous Bogolubov shift for field operators. The
condensate wave function, introduced in the course of this shift, has to satisfy the minimum
of the related thermodynamic potential, which is the stability condition necessary for mak-
ing the system stable and the theory conserving and self-consistent. At the same time, the
spectrum of collective excitations, according to the Hugenholtz-Pines theorem [15], has to
be gapless. The notorious problem is the appearance of the contradiction between the above
two requirements, when the theory is either nonconserving or gapful. This contradiction
does not arise only in the lowest orders with respect to particle interactions, when one uses
the Bogolubov approximation at low temperatures [11,12] or the quasiclassical approxima-
tion at high temperatures [16]. However, this contradiction immediately arises as soon as
the interaction strength is not asymptotically weak and one has to invoke a more elaborate
approximation. This problem of conserving versus gapless approximations was first empha-
sized by Hohenberg and Martin [17] and recently covered comprehensively by Andersen [9].
The problem is caused by the usage of nonrepresentative ensembles, which renders the sys-
tem unstable [18]. Here we show that employing a representative ensemble never yields the
above contradiction, always resulting in a self-consistent theory, being both conserving and
gapless.
We consider a system with the Hamiltonian
ψ†(r)
ψ(r) dr+
ψ†(r)ψ†(r′)Φ(r− r′)ψ(r′)ψ(r) drdr′ , (21)
in which the field operators ψ(r) = ψ(r, t) satisfy the Bose commutation relations, U =
U(r, t) is an external field, and Φ(r) = Φ(−r) is an interaction potential. For describing a
Bose-condensed system with broken global gauge symmetry, the Bogolubov shift [13,14] has
to be done through the replacement
ψ(r, t) −→ ψ̂(r, t) ≡ η(r, t) + ψ1(r, t) , (22)
where η(r, t) is the condensate wave function and ψ1(r, t) is the field operator of noncon-
densed particles. The latter field variables are assumed to be orthogonal to each other,
η∗(r, t)ψ1(r, t) dr = 0 . (23)
It is necessary to emphasize that the Bogolubov shift (22) realizes unitary nonequivalent
operator representations [18,19]. Accomplishing the Bogolubov shift (22) in Hamiltonian
(21), as well as in all operators of observables, we get the algebra of observables defined on
the Fock space F(ψ1) generated by the field operators ψ†1(r) (see details in Refs. [5,18,19]).
The condensate function is normalized to the number of condensed particles
|η(r, t)|2dr . (24)
The Bogolubov shift (22) is only rational when the number of condensed particles (24) is
macroscopic, which means that the limit
is not zero, where N is the total number of particles. The latter is given by the average
N = < N̂ > (25)
for the number-of-particle operator
ψ̂†(r)ψ̂(r) dr , (26)
in which the Bogolubov shift (22) is again assumed. The statistical averaging in Eq. (25)
and everywhere below is over the Fock space F(ψ1).
Substituting the Bogolubov shift (22) into Hamiltonian (21) gives in the latter the terms
linear in ψ1, because of which the average < ψ1 > can be nonzero. This, however, would
result in the nonconservation of quantum numbers, e.g., of spin or momentum. Therefore,
one has to impose the constraint for the conservation of quantum numbers,
< ψ1(r, t) > = 0 . (27)
Defining the self-adjoint condition operator
Λ̂(t) ≡
λ(r, t)ψ
1(r, t) + λ
∗(r, t)ψ1(r, t)
dr , (28)
in which λ(r, t) is a complex function, we may represent constraint (27) as the quantum
conservation condition
< Λ̂(t) > = 0 . (29)
In this way, there are three statistical conditions. The first condition is the normalization
(24) for the number of condensed particles. Condition (24) can be represented in the standard
form (4) by defining the operator
N̂0 ≡ 1̂
|η(r, t)|2dr , (30)
in which 1̂ is the unity operator in the Fock space F(ψ1). Then Eq. (24) reduces to the
statistical condition
N0 = < N̂0 > . (31)
The second condition is the normalization (25) for the total number of particles. Equivalently,
instead of normalization (25), we may consider the normalization condition
N1 = < N̂1 > , N̂1 ≡
1(r)ψ1(r) dr (32)
for the number of uncondensed particles N1 = N − N0. And the third condition is the
conservation condition (29). Respectively, the effective action, which is now a functional of
the two field variables, η(r, t) and ψ1(r, t), with taking account of the statistical conditions
(29), (31), and (32), becomes
A[η, ψ1] =
L̂+ µ0N̂0 + µ1N̂1 + Λ̂
dt . (33)
Here L̂ = L̂[ψ̂] is the Lagrangian (13) under the Bogolubov shift (22) and Λ̂ = Λ̂(t) from
Eq. (28). The quantities µ0, µ1, and λ(r, t) are the Lagrange multipliers guaranteeing the
validity of the corresponding statistical conditions. Introducing the grand Hamiltonian
H [η, ψ1] ≡ Ĥ − µ0N̂0 − µ1N̂1 − Λ̂ , (34)
in which Ĥ = Ĥ[ψ̂], with shift (22), and the effective Lagrangian
L[η, ψ1] ≡
η∗(r, t) i
η(r, t) + ψ
1(r, t) i
ψ1(r, t)
dr−H [η, ψ1] , (35)
for the action functional (33), we get
A[η, ψ1] =
L[η, ψ1] dt . (36)
The evolution equations follow from the extremization of the action functional (36), that
is, from the variations
δA[η, ψ1]
δη∗(r, t)
= 0 (37)
δA[η, ψ1]
1(r, t)
= 0 . (38)
These equations, as is clear from Eqs. (35) and (36), are equivalent to the equations of
motion
η(r, t) =
δH [η, ψ1]
δη∗(r, t)
ψ1(r, t) =
δH [η, ψ1]
1(r, t)
. (40)
One has to substitute here Hamiltonian (34) under the Bogolubov shift (22). Accomplishing
the variation in Eq. (39), we get
η(r, t) =
+ U − µ0
η(r)+
Φ(r− r′)
|η(r′)|2η(r) + X̂(r, r′)
dr′ , (41)
where the time dependence in the right-hand side, for short, is not explicitly shown, U =
U(r, t), and the notation for the correlation operator
X̂(r, r′) ≡ ψ†1(r′)ψ1(r′)η(r) + ψ
′)η(r′)ψ1(r)+
+η∗(r′)ψ1(r
′)ψ1(r) + ψ
′)ψ1(r
′)ψ1(r) (42)
is used. The variation in Eq. (40) gives
ψ1(r, t) =
+ U − µ1
ψ1(r)+
Φ(r− r′)
|η(r′)|2ψ1(r) + η∗(r′)η(r)ψ1(r′) + η(r′)η(r)ψ†1(r′) + X̂(r, r′)
dr′ . (43)
To get an equation for the condensate wave function, we have to take the statistical
average of Eq. (41). For this purpose, we introduce the normal density matrix
ρ1(r, r
′) ≡ < ψ†1(r′)ψ1(r) > , (44)
the anomalous density matrix
σ1(r, r
′) ≡ < ψ1(r′)ψ1(r) > , (45)
and their diagonal elements, giving the density of noncondensed particles
ρ1(r) ≡ ρ1(r, r) = < ψ†1(r)ψ1(r) > (46)
and the anomalous average
σ1(r) ≡ σ1(r, r) = < ψ1(r)ψ1(r) > . (47)
The quantity |σ1(r)| can be interpreted as the density of paired particles [19]. The total
density of particles
ρ(r) = ρ0(r) + ρ1(r) (48)
consists of the condensate density
ρ0(r) ≡ |η(r)|2 (49)
and the density of noncondensed particles (46). Averaging Eq. (41), we find the equation
for the condensate wave function
η(r, t) =
+ U − µ0
η(r)+
Φ(r− r′)
ρ(r′)η(r) + ρ1(r, r
′)η(r′) + σ1(r, r
′)η∗(r′)+ < ψ
′)ψ1(r
′)ψ1(r) >
dr′ . (50)
Equations (43) and (50) are the basic equations of motion for the field variables η(r, t) and
ψ1(r, t). These equations, according to Eqs. (39) and (40), are generated by the variation
of the grand Hamiltonian (34). The latter, in agreement with Eq. (20), also defines the
evolution operator Û(t), which satisfies the Schrödinger equation
Û(t) = H [η(r, 0), ψ1(r, 0)]Û(t) .
Thus, the representative ensemble for a Bose-condensed system is the triplet
{F(ψ1), ρ̂(0), Û(t)} .
It is important to stress that the so defined representative ensemble possesses a principal
feature making it different from the standardly used ensemble having the sole Lagrange
multiplier µ0 ≡ µ1. But then the normalization condition (24) cannot be guaranteed. Then
the evolution equation for the condensate wave function is not a result of a variational
procedure. For an equilibrium system, this means that the number of condensed particles
N0 does not provide the minimum of a thermodynamic potential, which implies the system
instability. All notorious inconsistencies in theory, manifesting themselves in the lack od
conservation laws or in the appearance of an unphysical gap in the spectrum, are caused by
the usage of nonrepresentative ensembles.
4 Green Functions
The equations of motion (43) and (50) allow us to derive the evolutional equations for the
Green functions. To this end, we shall use the compact notation denoting the set {rj, tj} by
the sole letter j, so that the dependence of functions on the spatial and temporal variables
looks like
f(12 . . . n) ≡ f(r1, t1, r2, t2, . . . , rn, tn) .
The product of the differentials drjdtj will be denoted as d(j), so that
d(12 . . . n) ≡
drj dtj .
We shall employ the Dirac delta function
δ(12) ≡ δ(r1 − r2) δ(t1 − t2) .
For the interaction potential, we shall use the retarded form
Φ(12) ≡ Φ(r1 − r2)δ(t1 − t2 + 0) . (51)
The matrix Green function G(12) = [Gαβ(12)] is a 2 × 2 matrix, with α, β = 1, 2, and
with the following elements:
G11(12) ≡ −i < T̂ψ1(1)ψ†1(2) > , G12(12) ≡ −i < T̂ψ1(1)ψ1(2) > ,
G21(12) ≡ −i < T̂ψ†1(1)ψ
1(2) > , G22(12) ≡ −i < T̂ψ
1(1)ψ1(2) > , (52)
where T̂ is the time-ordering operator.
Let us introduce the operator
K̂j ≡ −
+ U(j)− µ1 , (53)
the condensate effective potential
V (12) ≡ δ(12)
Φ(13)|η(3)|2d(3) + Φ(12)η(1)η∗(2) , (54)
and let us rewrite the correlation operator (42) in the form
X̂(12) = ψ
1(2)ψ1(2)η(1) + ψ
1(2)η(2)ψ1(1) + η
∗(2)ψ1(2)ψ1(1) + ψ
1(2)ψ1(2)ψ1(1) . (55)
We also define the matrix correlation function X(123) = [Xαβ(123)] with the elements:
X11(123) ≡ − < T̂ X̂(12)ψ†1(3) > , X12(123) ≡ − < T̂ X̂(12)ψ1(3) > ,
X21(123) ≡ − < T̂X̂+(12)ψ†1(3) > , X22(123) ≡ − < T̂X̂+(12)ψ1(3) > . (56)
From the equations of motion (43) and (50), we find the equations
− K̂1
G11(12)−
V (13)G11(32) d(3)−
Φ(13) [η(1)η(3)G21(32) + iX11(132)] d(3) = δ(12) ,
− K̂1
G12(12)−
V (13)G12(32) d(3)−
Φ(13) [η(1)η(3)G22(32) + iX12(132)] d(3) = 0 ,
− K̂1
G21(12)−
V ∗(13)G21(32) d(3)−
Φ(13) [η∗(1)η∗(3)G11(32) + iX21(132)] d(3) = 0 ,
− K̂1
G22(12)−
V ∗(13)G22(32) d(3)−
Φ(13) [η∗(1)η∗(3)G12(32) + iX22(132)] d(3) = δ(12) . (57)
The self-energy Σ(12) = [Σαβ(12)] is a matrix whose elements are defined by the relations
[Σ11(13)G11(32) + Σ12(13)G21(32)] d(3) =
V (13)G11(32) d(3) +
Φ(13) [η(1)η(3)G21(32) + iX11(132)] d(3) ,
[Σ11(13)G12(32) + Σ12(13)G22(32)] d(3) =
V (13)G12(32) d(3) +
Φ(13) [η(1)η(3)G22(32) + iX12(132)] d(3) ,
[Σ21(13)G11(32) + Σ22(13)G21(32)] d(3) =
V ∗(13)G21(32) d(3) +
Φ(13) [η∗(1)η∗(3)G11(32) + iX21(132)] d(3) ,
[Σ21(13)G12(32) + Σ22(13)G22(32)] d(3) =
V ∗(13)G22(32) d(3) +
Φ(13) [η∗(1)η∗(3)G12(32) + iX22(132)] d(3) . (58)
Let us introduce the matrix condensate propagator C(12) = [Cαβ(12)], with the elements
C11(12) ≡ −iη(1)η∗(2) , C12(12) ≡ −iη(1)η(2) ,
C21(12) ≡ −iη∗(1)η∗(2) , C22(12) ≡ −iη∗(1)η(2) , (59)
The latter have the properties
C11(21) = C22(12) , C11(11) = C22(11) , C12(21) = C12(12) ,
C21(21) = C21(12) , C
11(12) = −C22(12) , C∗12(12) = −C21(12) . (60)
The binary Green function is a matrix B(123) = [Bαβ(123)],
B(123) = C(12)G(23)− iρ0(2)G(13) +X(123) , (61)
whose elements are
B11(123) ≡ C11(12)G11(23) + C11(22)G11(13) + C12(12)G21(23) +X11(123) ,
B12(123) ≡ C11(12)G12(23) + C11(22)G12(13) + C12(12)G22(23) +X12(123) ,
B21(123) ≡ C22(12)G21(23) + C22(22)G21(13) + C21(12)G11(23) +X21(123) ,
B22(123) ≡ C22(12)G22(23) + C22(22)G22(13) + C21(12)G12(23) +X22(123) , (62)
and where ρ0(1) ≡ |η(1)|2.
With Eqs. (59) and (61), relations (58), defining the self-energy, can be rewritten in the
matrix form
Σ(13)G(32) d(3) = i
Φ(13)B(132) d(3) . (63)
Then the equations of motion (57) acquire the matrix representation
τ̂3 i
− K̂1
G(12)−
Σ(13)G(32) d(3) = δ(12) , (64)
in which the delta function δ(12) in the right-hand side is assumed to be factored with the
unity matrix 1̂ = [δαβ ] and
τ̂3 ≡
is a Pauli matrix.
Introducing the inverse propagator
G−1(12) ≡
τ̂3 i
− K̂1
δ(12)− Σ(12) (65)
allows us to transform Eq. (64) into
G−1(13)G(32) d(3) = δ(12) . (66)
An equivalent representation, following from Eq. (66), is
G(13)G−1(32) d(3) = δ(12) . (67)
For the self-energy, using Eq. (63), we have
Σ(12) = i
Φ(13)B(134)G−1(42) d(34) . (68)
The equations for the Green functions are to be complimented by the equation for the
condensate wave function (50), which, introducing one more anomalous average
ξ(12) ≡ < ψ†1(2)ψ1(2)ψ(1) > , (69)
can be represented as
η(1) =
+ U(1)− µ0
η(1)+
Φ(12) [ρ(2)η(1) + ρ1(12)η(2) + σ1(12)η
∗(2) + ξ(12)] d(2) . (70)
It is the equations for the Green functions and the equation for the condensate func-
tion, which become mutually incompatible in the standard approach, while employing the
representative ensemble renders the theory self-consistent in any approximation.
5 Theory Self-Consistency
One usually confronts inconsistencies in theory considering a uniform equilibrium Bose-
condensed system. Then, in any given approximation, one gets either a nonconserving
theory, that is, an unstable system, or one finds an unphysical gap in the spectrum, which,
actually, again corresponds to an unstable system [9,18]. To analyze this problem, we pass
now to the case of an equilibrium uniform system, when U = 0.
Then we use the Fourier transform for the Green function
G(12) =
G(k, ω)ei(k·r12−ωt12)
(2π)4
in which
r12 ≡ r1 − r2 , t12 ≡ t1 − t2 .
By their definition in Eq. (52), the Green function elements possess the properties
G11(21) = G22(12) , G12(21) = G12(12) , G21(21) = G21(12) . (71)
Therefore the corresponding Fourier transforms satisfy the relations
G11(−k,−ω) = G22(k, ω) , G12(−k,−ω) = G12(k, ω) ,
G21(−k,−ω) = G21(k, ω) . (72)
Assuming that the system is isotropic, one has
Gαβ(−k, ω) = Gαβ(k, ω) (73)
for all α, β. Combining Eqs. (72) and (73), we find
G11(k,−ω) = G22(k, ω) , G12(k,−ω) = G12(k, ω) ,
G21(k,−ω) = G21(k, ω) . (74)
Also, for a uniform equilibrium system, one has [14] the equality
G21(k, ω) = G12(k, ω) . (75)
Fourier-transforming the self-energy
Σ(12) =
Σ(k, ω)ei(k·r12−ωt12)
(2π)4
and, similarly, the inverse propagator (65), we have for the latter
G−1(k, ω) = τ̂3ω −
+ µ− Σ(k, ω) . (76)
Then Eq. (67) reduces to
G−1(k, ω)G(k, ω) = 1̂ , (77)
where 1̂ = [δαβ ].
From Eqs. (76) and (77), it follows that G−1(k, ω), hense, also Σ(k, ω), have the same
symmetry properties as G(k, ω). In particular,
Σαβ(−k, ω) = Σαβ(k, ω) , Σ12(k,−ω) = Σ12(k, ω) , Σ21(k,−ω) = Σ21(k, ω) ,
Σ21(k, ω) = Σ12(k, ω) , Σ11(k,−ω) = Σ22(k, ω) . (78)
The matrix equation (77), explicitly, is the system of equations
ω − k
+ µ1 − Σ11
G11 − Σ12G21 = 1 ,
+ µ1 − Σ11
G12 − Σ12G22 = 0 ,
−ω − k
+ µ1 − Σ22
G21 − Σ21G11 = 0 ,
−ω − k
+ µ1 − Σ22
G22 − Σ21G12 = 1 , (79)
where, for short, Gαβ = Gαβ(k, ω) and Σαβ = Σαβ(k, ω). The solutions to these equations
G11(k, ω) =
ω + k2/2m+ Σ11(k, ω)− µ1
D(k, ω)
G12(k, ω) = −
Σ12(k, ω)
D(k, ω)
, (80)
with the denominator
D(k, ω) ≡
ω − k
− Σ11(k, ω) + µ1
+ Σ22(k, ω)− µ1
+ Σ212(k, ω) . (81)
The solutions for Σ21(k, ω) and G22(k, ω) are defined by the symmetry properties (72) to
(75).
The excitation spectrum is given by the poles of the Green functions, that is, by the zero
of denominator (81),
D(k, εk) = 0 . (82)
Equation (82) can be represented as
[Σ11(k, εk)− Σ22(k, εk)]±
ω2k − Σ212(k, εk) , (83)
with the notation
[Σ11(k, εk) + Σ22(k, εk)]− µ1 . (84)
Denominator (81) enjoys the property
D(k,−ω) = D(k, ω) .
Consequently, if εk is a solution of Eq. (82), then −εk is also its solution, which is in
agreement with the form of Eq. (83).
For an equilibrium uniform system, the Bogolubov shift (22) is equivalent to the sepa-
ration of the zero-momentum term in the expansion of the field operator over plane waves.
The shift itself has meaning only under the normalization condition (24), in which N0 ∼ N ,
that is, the zero-momentum state is macroscopically occupied. The latter becomes possible
when the single particle spectrum touches zero. Therefore, the necessary condition for the
existence of Bose-Einstein condensate is
εk = 0 . (85)
This is to be complimented by the stability condition
Re εk ≥ 0 , Im εk ≤ 0 . (86)
This condition should be kept in mind when choosing the sign plus in front of the square
root in spectrum (83).
Taking limit (85) for spectrum (83), we notice that, according to properties (78),
Σ11(k, 0) = Σ22(k, 0) . (87)
By using perturbation theory for a stable system, one can show [15] that in all orders of the
theory
Σαβ(0, 0) ≥ 0 . (88)
Then the necessary condition (85) yields the expression for the chemical potential
µ1 = Σ11(0, 0)− Σ12(0, 0) , (89)
which is the Hugenholtz-Pines relation [15].
On the other hand, we have Eq. (70) for the condensate wave function. For an equilibrium
uniform system, with no external potential U , all densities do not depend on the spatial and
temporal variables,
ρ0(r) = ρ0 , ρ1(r) = ρ1 , σ1(r) = σ1 , ρ(r) = ρ . (90)
The condensate wave function reduces to the constant
η(r, t) = η =
ρ0 . (91)
Then we substitute into Eq. (70) the Fourier transforms for the interaction potential
Φ(r) =
ik·r dk
(2π)3
for the normal density matrix (44),
ρ1(r1, r2) =
ik·r12
(2π)3
and for the anomalous density matrix (45),
σ1(r1, r2) =
ik·r12
(2π)3
Similarly, the Fourier transform for the anomalous average (69) is
ξ1(r1, r2) =
ik·r12
(2π)3
As a result, Eq. (70) gives
µ0 = ρΦ0 +
nk + σk +
(2π)3
. (92)
Generally, expressions (92) and (89) do not coincide with each other, their difference
being
µ0 − µ1 = ρΦ0 +
nk + σk +
(2π)3
−Σ11(0, 0) + Σ12(0, 0) . (93)
This is the general expression for the difference between the Lagrange multipliers µ0 and µ1
for an arbitrary equilibrium uniform Bose-condensed system.
Usually, one does not distinguish between the Lagrange multipliers µ0 and µ1, which
implies setting µ0 − µ1 → 0. However, as is evident from Eq. (93), there is no any reason
for requiring that this quantity be zero. As an illustration, we may resort to the Hartree-
Fock-Bogolubov approximation, in which ξk = 0 and
Σ11(0, 0) = (ρ+ ρ0)Φ0 +
(2π)3
, Σ12(0, 0) = ρ0Φ0 +
(2π)3
Relation (89) then yields
µ1 = ρΦ0 +
(nk − σk)Φk
(2π)3
. (94)
The difference of the chemical potentials (93) becomes
µ0 − µ1 = 2
(2π)3
, (95)
which, certainly, is nonzero [20].
In this way, the introduction of the additional Lagrange multiplier makes the theory
completely self-consistent. All inconsistencies that often arise in other works, such as the
appearance of a gap in the spectrum, system instability or a distortion of the phase transition
order, are caused by neglecting the difference between the multiplier µ1 and the multiplier
µ0. It is worth emphasizing that the introduction of the Lagrange multiplier µ0 for preserving
the normalization condition (24), from the mathematical point of view is strictly necessary.
In other case, the employed ensemble would not be representative, hence, could not correctly
describe the Bose-condensed system with broken gauge symmetry.
Acknowledgement
I am grateful for many useful discussions to M. Girardeau, R. Graham, H. Kleinert, and
E.P. Yukalova.
References
[1] J.W. Gibbs, Collected Works (Longmans, New York, 1931), Vol. 2.
[2] D. ter Haar, Elements of Statistical Mechanics (Reinhart, New York, 1954).
[3] D. ter Haar, Rep. Prog. Phys. 24, 304 (1961).
[4] V.I. Yukalov, Phys. Rep. 208, 395 (1991).
[5] V.I. Yukalov, Statistical Green’s Functions (Queen’s University, Kingston, 1998).
[6] H. Kleinert, Path Integrals (World Scientific, Singapore, 2004).
[7] P.W. Courteille, V.S. Bagnato and V.I. Yukalov, Laser Phys. 11, 659 (2001).
[8] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation (Clarendon, Oxford, 2003).
[9] J.O. Andersen, Rev. Mod. Phys. 76, 599 (2004).
[10] K. Bongs and K. Sengstock, Rep. Prog. Phys. 67, 907 (2004).
[11] N.N. Bogolubov, J. Phys. (Moscow) 11, 23 (1947).
[12] N.N. Bogolubov, Moscow Univ. Phys. Bull. 7, 43 (1947).
[13] N.N. Bogolubov, Lectures on Quantum Statistics (Gordon and Breach, New York, 1967),
Vol. 1.
[14] N.N. Bogolubov, Lectures on Quantum Statistics (Gordon and Breach, New York, 1970),
Vol. 2.
[15] N.M. Hugenholtz and D. Pines, Phys. Rev. 116, 489 (1959).
[16] N. Prokofiev, O. Ruebenacker and B. Svistunov, Phys. Rev. A 69, 053625 (2004).
[17] P.C. Hohenberg and P.C. Martin, Ann. Phys. 34, 291 (1965).
[18] V.I. Yukalov, Phys. Rev. E 72, 066119 (2005).
[19] V.I. Yukalov, Laser Phys. 16, 511 (2006).
[20] V.I. Yukalov and E.P. Yukalova, Laser Phys. Lett. 2, 506 (2005).
