The time-dependent Hartree-Fock-Bogoliubov equations for Bosons

In this article, we use quasifree reduction to derive the time-dependent Hartree-Fock-Bogoliubov (HFB) equations describing the dynamics of quantum fluctuations around a Bose-Einstein condensate in $\mathbb R^d$. We prove global well-posedness for the HFB equations for sufficiently regular pair interaction potentials, and establish key conservation laws. Moreover, we show that the solutions to the HFB equations exhibit a symplectic structure, and have a form reminiscent of a Hamiltonian system. In particular, this is used to relate the HFB equations to the HFB eigenvalue equations encountered in the physics literature. Furthermore, we construct the Gibbs states at positive temperature associated with the HFB equations, and establish criteria for the emergence of Bose-Einstein condensation.


Introduction
In this paper, we study the time-dependent generalization of the Hartree-Fock-Bogoliubov (HFB) equations which describe the quantum fluctuations of the Bose field around a Bose-Einstein condensate.
We assume that (a) the external potential V is infinitesimally form bounded with respect to the Laplacian −∆ and (b) v 2 ≤ C(1 − ∆) in the sense of quadratic 1 forms, for some constant C < ∞. These conditions imply that H is selfadjoint on the domain of the operator H 0 := dx ψ * (x)(−∆)ψ(x) (see Appendix D). States of our system are positive linear ('expectation') functionals ω on the Weyl CCR algebra W over Schwartz space S(R d ). They can correspond either to nonzero densities and temperatures, or to a finite number of particles as in the case of BEC experiments in traps. In the latter case they are given by density operators on Fock space F , ω(A) = Tr(AD), for all observables A, where D is a positive, trace-class operator with unit trace on F .
It is convenient to define states ω on products ψ # (f 1 ) . . . ψ # (f n ) of creation and annihilation operators. This is done using derivatives ∂ s k of expectation values ω(W (s 1 f 1 ) . . . W (s n f n )) of the Weyl operators W (f ) := e iφ(f ) , with φ(f ) := ψ * (f ) + ψ(f ) (see [9], Section 5.2.3). We always assume that the states we consider are such that such derivatives (and therefore the corresponding correlation functions) exist, to arbitrary order. For n ≤ 4, which is our case, this is guaranteed by assuming that ω(N 2 ) < ∞, where N is the number operator N := dx ψ * (x)ψ(x). This implies, in particular, that ω is given by a density operator.
The evolution of states is given by the von Neumann-Landau equation ( [31,19], see also [33,5] and, for some history, [37]) for all observables A. We refer to [9,25] for a mathematical justification of this equation.
1.1. Quasifree States and Truncated expectations. As the evolution (4) is extremely complicated, one is interested in manageable approximations. The natural and most commonly used approximation is given in terms of quasifree states, the simplest class of states.
A state ω is called quasifree if the truncated expectations ω T (ψ 1 , . . . , ψ n ) vanish for all n > 2. We denote quasifree states by ω q and the set of quasifree states, by Q ⊆ S.
1.2. Quasifree reduction of the full dynamics. As was mentioned above, the detailed properties of the dynamics of the many-body system described by (4) are rather complicated, and an approximation is required to extract some key qualitative features. The main idea is to restrict the dynamics to quasifree states.
However, the property of being quasifree is not preserved by the dynamics given by (4) and the main question here is how to project the true quantum evolution onto the class of quasifree states.
The projection (or reduction) we propose is to map the solution ω t of (4), with a quasifree initial state ω 0 ∈ Q, to the family (ω q t ) t≥0 ∈ C 1 R + 0 ; Q of quasifree states satisfying for all observables A, which are at most quadratic in the creation and annihilation operators. (For the Hamiltonian H given by (1), the commutator [A, H] contains products of at most 4 creation and annihilation operators.) This is the quasifree reduction of equation (4).
As was mentioned above, a quasifree state ω q determines and is determined by the truncated expectations to the second order: Let γ and σ denote the operators with the integral kernels γ(x, y) and σ(x, y). This definition implies that γ = γ * ≥ 0 and σ * =σ, (9) whereσ = CσC with C being the complex conjugation. (More detailed characterizations of γ and σ are given in Proposition 3.1.) yields a system of coupled nonlinear PDE's for (φ t , γ t , σ t ), the Hartree-Fock-Bogoliubov (HFB) equations. Since quasifree states are characterized by their truncated expectations φ, γ, σ, this system of equations is equivalent to equation (7).
To give a flavor of the HFB equations at this point, we formally assume the pair interaction potential to be a delta distribution, v(x) = gδ(x) with coupling constant g ≥ 0. The HFB equations then have the form where [A, B] + = AB T + BA T , A T := CA * C and d(σ)(x) := σ(x, x) and Here and in what follows, we denote the multiplication operators and the functions by which they multiply by the same symbols. The difference is always clear from context.
The physical interpretation of the truncated expectations of ω q t is that φ t gives the quantum mechanical wave function for the Bose-Einstein condensate, while γ t and σ t describe the dynamics of sound waves in the quasifree approximation (in particular, γ t gives the density of the thermal cloud of atoms). (In the physics literature, n = d(γ) and m = d(σ) are called the non-condensate density and "anomalous" density, respectively.) The HFB equations provide a time-dependent extension of the standard stationary Hartree-Fock-Bogoliubov equations for a Bose gas appearing in the physics literature, see e.g. [14,15,32]. Related equations (with φ t = 0) appear in superconductivity (Bogoliubov-de Gennes equations), where they are equivalent to the BCS effective Hamiltonian description.
1.3. Main Results. The quasifree reduction of the dynamics (4) and evaluation of the resulting equations for the truncated expectations φ, γ, σ -the HFB equations -are among the main results of this paper (see Theorem 2.2).
We show that the HFB equations are equivalent to the self-consistent equation where H hfb (ω q ) is an explicitly constructed quadratic Hamiltonian (see (31)), for all observables A. See Theorem 2.4.
Equation (15) suggests to define the HFB stationary states as the quasifree states satisfying ω q ([A, H hfb (ω q )]) = 0, for all observables A. (If ω q is given by a density matrix, we can rewrite this equation as a fixed point problem, see (16) below.) The most interesting among such states are the ground or Gibbs states. These states are defined as where ω q L is the quasifree ground or Gibbs state of the Bose gas confined to a torus, Λ L = R d /2LZ d , i.e., the box [−L, L] d , with periodic boundary conditions. It satisfies the fixed point equation where β > 0 is the inverse temperature, µ is the chemical potential, and Ξ = Tr[exp(−β(H hf b (ω q L )) − µN)], and is a stationary solution to the equation (15) on Λ L . (To prove the second statement, , A]) = 0, for any observable A defined on some torus Λ L .) We also initiate a mathematical discussion of the HFB equations. In particular, for γ trace-class and σ Hilbert-Schmidt operators (for the precise restriction see Section 2), we show • Global well-posedness (Theorem 5.1) of the HFB equations; • Conservation of the total particle number where N is the number operator (see Corollary 2.7); • Existence and conservation (under certain conditions on v and ω q t ) of the energy: (see Corollary 2.7 and Theorem 2.8, or Prop 3.12).
More generally, any observable conserved by the von Neumann-Landau dynamics and which is at most quadratic in the creation and annihilation operators is also conserved by the quasifree dynamics. See Theorem 2.5. In the special case of the observable N, this gives the statement above.
Note that conservation of the total particle number is related to the U (1)-gauge invariance, i.e., invariance under the transformation ψ ♯ → (e iθ ψ) ♯ , of the Hamiltonian H.
In the thermodynamic limit, with V = 0, one should replace the energy and the number of particles by the energy and particle number densities.
For V = 0 and γ and σ translationally invariant, our full energy functional (see (35)) reduces to the energy density functional going back to the work [13] and considered in [28,29]. It is shown in the latter papers that this functional has minimizers under the constraint of constant entropy and particle densities. [28,29] show also the appearance of a Bogoliubov condensate for the corresponding minimizers. (One still has to show that states thus obtained are stationary solutions to the equation (15).) In this paper we do not consider the general problem of existence of static solutions. However, for V = 0, we present a result on the existence of the positive temperature, U (1)-gauge and translation invariant Gibbs HFB states and show that Bose-Einstein condensation (BEC) occurs above a critical density, see Theorem 6.3.
As was mentioned above, the U (1)-gauge invariant states have φ = 0 and σ = 0 and therefore are, in fact, stationary solutions of the bosonic Hartree-Fock equation. Moreover, as the results of [28,29] show, in the BEC regime, these states are not minimizers of the full HFB (mean) energy for the (mean) entropy and number of particles fixed. However, the existence of such states exhibiting Bose-Einstein condensation indicates that there are also symmetry breaking Gibbs HFB states, i.e. with φ = 0 and σ = 0.

Fixed Point Equation.
Let U ω q (t, s) denote the unitary propagator on the Fock space solving Then, we can rewrite the equation (15) with an initial condition ω q 0 as an 'integral' equation where A is an arbitrary observable. This is the fixed point problem, 0)). Since U ω q (t, s) are generated by quadratic Hamiltonians, we have that ω q 0 (U ω q (t, 0) * A U ω q (t, 0)) is quasifree for any time t > 0. This formulation opens an opportunity to proving the existence of the quasifree dynamics directly, without going to the truncated expectations.
In this article, we do not show that the HFB equations provide an accurate approximation of the many-body dynamics (4) for finite times. There is a considerable literature on the derivation of the simpler Hartree and Hartree -Fock equations. Recently, the Hartree equation with linear fluctuations around the Hartree solutions (producing the linearized HFB equations with respect to γ and σ) were derived in [16,23,26,22,21], see [20] for a recent review.
1.5. Organization of the paper. In Section 2, we present the HFB equations, which we prove in Appendix B using the method of quasifree reduction, applied to the interacting Bose field in R d characterized by the Hamiltonian (1). We also show that conservation laws for the many-body problem imply conservation laws for the HFB equation. In Section 3, we show that the solutions to the HFB equations possess a symplectic structure, and that they have similarities with a Hamiltonian system. In Section 4 we show how the symplectic version of the HFB equations and the HFB eigenvalue equations found in the physics literature are related. In Section 5, we prove that the Cauchy problem for the HFB equations is globally well-posed in the energy space, provided that the pair interaction potential is sufficiently regular.
Our proof of global well-posedness is in part inspired by previous works on the Hartree-Fock equation [6,12,7,11,38]. In Section 6, we give a proof of the Bose-Einstein condensation in the stationary case. A brief account on quasifree states along with the proofs of various technical lemmata are collected in the Appendices.

The HFB equations and their basic properties
In this section, we formulate the HFB equations for a general pair potential v and prove the associated conservation laws. The derivation of the HFB equations is done in Appendix B by applying the quasifree reduction as in the introduction.
with ∆ x being the Laplacian in d dimensions. We denote by B the space of bounded operators on L 2 (R d ), with the operator norm denoted by · . For j ∈ N 0 we define the spaces with H j the Sobolev space H j (R d ) and B j = M −j BM −j . The norms on X j,∞ are given by the norms on the Banach spaces, H j × B j × B j , i.e., (The superindex ∞ indicates that we are dealing with bounded operators, as opposed to the trace-class and Hilbert-Schmidt ones appearing later.) Moreover, we let X ∞ T := C 0 ([0, T ); X j,∞ ) ∩ C 1 ([0, T ); X 0,∞ ), for a fixed j satisfying j > d/2 and j ≥ 2, and denote by X j,∞ qf and X ∞ T,qf spaces of quasifree states and families of quasifree states with the 1 st and 2 nd order truncated expectations from the spaces X j,∞ and X ∞ T , respectively.
The reason for introducing the spaces B j is the following elementary result.
Let v ∈ L 1 . Then the operator k : α → v ♯ α, defined through its integral kernel: Similarly, one shows the Hölder continuity.
For the second statement, the result above gives In what follows we use the same notation for functions and the operators of multiplication by these functions. Which one is meant in every instance is clear from the context.
with the Hamiltonian H defined in (1), if and only if the triple (φ t , γ t , σ t ) ∈ X ∞ T of the 1 st and 2 nd order truncated expectations of ω q t satisfies the time-dependent Hartree-Fock-Bogoliubov equations where If v = gδ, h(γ) agrees with h gδ (γ) in (13), and k(σ) agrees with the multiplication operator by g d(σ)(x) in (13).
Due to Lemma 2.1 and the fact that h(γ t ) is ∆−bounded, for each t > 0, the r.h.s. of (26) - (28) belongs to the space X 0,∞ . The proof of Theorem 2.2 is given in Appendix B.
Remark 2.3. If we base our spaces for γ and σ on the trace-class and Hilbert-Schmidt operators, instead of bounded ones, then we can relax the conditions on the potentials.
For a quasifree state ω q with 1 st and 2 nd order truncated expectations (φ, γ, σ) ∈ X 1 , we define the quadratic Hamiltonian parametrized by (φ, γ, σ) as (28) if and only if the corresponding quasifree state ω q t ∈ X ∞ T,qf satisfies the equation The proof of Theorem 2.4 is given in Appendix C.
We now prove the conservation laws for the number of particles (or more generally, for any observable commuting with the Hamiltonian H which is quadratic with respect to creation and annihilation operators), and for the energy.
Proof. This follows from (25) for A of order up to two, with [A, H] = 0.
To draw some consequences from this result we need to define additional spaces.
Remark 2.6. We identify Hilbert-Schmidt operators on L 2 (R d ), denoted by L 2 , with their kernels in L 2 (R 2d ), if no confusion may arise.
Notations. We denote by L 1 the space of trace-class operators on L 2 (R d ) endowed with the trace norm · L 1 . For j ∈ N 0 we define the spaces with H j being the Sobolev space . These are real Banach spaces, similar to the Sobolev spaces H j (R d ) in the scalar case, when endowed with the norms Furthermore, we let X T := C 0 ([0, T ); X 3 ) ∩ C 1 ([0, T ); X 1 ) and, as above, we denote by X j qf and X qf T the spaces of quasifree states and families of quasifree states with the 1 st and 2 nd order truncated expectations from the spaces X j and X T , respectively. (32)). Then the number of particles If, in addition, v is as in Theorem 2.4, then the energy ω q t (H) is conserved. Theorem 2.8. Let v be as in Theorem 2.4 and ω q ∈ X qf . Then the energy ω q (H) = E(φ, γ, σ) is given explicitly as Proof. We use where the Weyl operators are defined through W φ = exp ψ * (φ) − ψ(φ) and satisfy Note that the state ω q Ct is quasifree because ω q is quasifree. By construction ω q C (ψ(x)) = 0 and thus using (5) and the quasifreeness of ω q C one sees that ω q C vanishes on monomials of odd order in the creation and annihilation operators.
, hence using the vanishing on monomials of odd order in the creation and annihilation operators Then, using that ω q C is a quasifree state with expectations (0, γ, σ) yields which gives the expression of the energy in terms of φ, γ and σ.

Generalized One-particle Density Matrix and Bogolubov Transforms
In this section, we consider the HFB equations (27) -(28) for γ t and σ t and reformulate them in terms the generalized one-particle density matrix Γ t = γt σt σt 1+γt .
We show that the diagonalizing maps for Γ t are symplectomorphisms (see below for the definition) and that the resulting equation for Γ t is equivalent to the evolution equation for these symplectomorphisms. The latter will allow us to (a) give another proof of the conservation of energy without using to the second quantization framework and (b) connect the time-dependent HFB equations (27) -(28) to the time-independent HFB equations used in the physics literature. See Section 4.
Proposition 3.1. The generalized one-particle density matrix, Γ, satisfies: This property is equivalent to the following statements: (Statement (4) follows from (1) and (3) and is given here for later convenience of references.) Proof. We remark that the truncated expectations γ and σ are the expectations of the state Statements (1) and (2) are obvious. The inequality in Point (3) follows from the Schur complement argument: Finally, we observe that (1) and (3) and the inequality γ ≤ Tr[γ]1 imply the following bound on σσ * , Inserting M = √ 1 − ∆ x on both sides and taking the trace yields (4).
Notations. For j ∈ N 0 we define the spaces where the spaces H j and H j s are defined after (34). The norms on Y j are given by In what follows we fix a number T > 0 and a family φ (26)) and do not display it in our notation. A simple computation yields the first result of this section: with (29) and (30), To formulate the next result we introduce some definitions.
If, moreover, there exists a unitary transformation U on the Fock space, sometimes called implementation of U, such that then the symplectomorphism U is said to be implementable.
Remark 3.4. The operator U is a symplectomorphism in the sense that it preserves the symplectic form ℑ · , S · on h × h (i.e. is a canonical map). (In fact, U preserves · , S · .) Remark 3.5. The operator U is a symplectomorphism if and only if the operator f → uf + vf is a symplectomorphism on (h, ℑ ·, · ) in the usual sense (i.e., it preserves the symplectic form ℑ ·, · ) Remark 3.6. The conditions in (42) are equivalent to satisfying the four equations Remark 3.7. The transformation is called the Bogoliubov transformation. It is easy to check that it preserves the CCR iff the operator U = u v vū satisfies (42).
If v is Hilbert-Schmidt, then the Bogoliubov transformation (44) is implementable. This condition is referred to as the Shale condition; see [36].
For later use, we introduce the Banach space using the same identification between operators and kernels as before.
We begin with an auxiliary result: where 0 ≤ γ ′ ≤ γ. The operator γ ′ is unique up to conjugation by a unitary operator.
This result is related to Theorem 1 of [27], which is stronger. See also [2,3]. As the relation between the two results is not obvious, we give a direct proof of Proposition 3.8 after the proof of Proposition 3.9.
The next result relates the evolution of Γ t to the evolution of implementable symplectomorphisms U t ∈ H ∞,2 (h × h), diagonalizing Γ t . Proposition 3.9. (i) For any Γ t ∈Ỹ T and any implementable symplectomorphism U 0 ∈ H ∞,2 , the initial value problem has a unique solution in H ∞,2 , which is a symplectomorphism for every t.
(ii) Let Γ t ∈Ỹ T solve the equation (41), with an initial condition Γ 0 ∈Ỹ 3 , s.t. Γ 0 ≥ 0. Let U 0 be an implementable symplectomorphism diagonalizing Γ 0 : Then the continuous family of implementable symplectomorphisms U t in H ∞,2 (h×h) satisfying (45), with the above U 0 , diagonalizes Γ t : Proof of Prop. 3.9. The operator Λ t can be decomposed as Λ t = Λ 1 + Λ 2,t with The first operator, Λ 1 , is the generator of a continuous one-parameter group in H ∞,2 . As for the second one, using the continuity of t → ρ t ∈ X 1 , and Prop. E.1, we get the continuity of t → Λ 2,t ∈ H ∞,2 . We can thus use classical results of functional analysis (see, e.g., [18]) to obtain the existence and uniqueness of U t and its regularity.
The same arguments as in the next lemma prove that U t is a symplectomorphism.
Finally, Γ t and U * t Γ 0 U t satisfy the same differential equation, and the uniqueness of a solution to (45) proves the last equality.
Proof of existence in Prop. 3.8. We split the proof into two lemmas, Lemmas 3.10 and 3.11 below. The strategy is to construct Γ t and symplectomorphisms U t such that U t Γ t U * t = Γ 0 , for all t, and in the limit t → ∞, Γ ∞ has the desired form. The key step will be to use a differential equation for Γ t implying σ t H 1 s ց 0. Lemma 3.10. Let T > 0 and Λ t = at bt btāt ∈ C([0, T ); H ∞,2 ). Then, the ordinary differential equation with initial data U * 0 = 1 0 0 1 , has a unique global solution U t ∈ C 1 ([0, T ); H ∞,2 ), and U t is a symplectomorphism for all time.
with initial data σ 0 = σ, γ 0 = γ given in Prop 3.9(i), then, for all time t, Proof. The existence and uniqueness of U * t follows from the theory of time-dependent linear ordinary differential equations once one observes that H 1 and H 1 s are continuously embedded in B(H 1 ) and We choose a t and b t in (48) and (49) such that σ t vanishes in the limit t → ∞. Let L 1 (h) and L 2 (h) denote the spaces of trace-class and Hilbert -Schmidt operators on the space h.
Lemma 3.11. The ordinary differential equation with initial data σ 0 = σ, γ 0 = γ given in Prop. 3.9(i), has a unique global solution −iσt 0 , and U t = ut vt vtūt and Γ t = γt σt σt 1+γt as in Lemma 3.10: Proof. The existence of maximal solutions to (51) -(52) follows from the Picard-Lindelöf theorem. Now using the U t constructed in Lemma (3.10), one gets that which implies that Γ t ≥ 0 and thus γ t ≥ 0. It then follows from (51) that γ t is decreasing in the sense of quadratic forms and γ t H 1 ≤ γ 0 H 1 . Using (52) and γ t ≥ 0 yields s and the maximal time of the solution is T = ∞. We also get that γ t → γ ∞ in H 1 as t → ∞ as γ t is decreasing and bounded by below, and σ t → 0.
Integrating the derivative of U * t and taking the norm of both sides yields The Grönwall lemma, combined with U * 0 H ∞,2 = 1 and the estimate on σ t H 1 Thus, the integral ∞ 0 SΛ s U * s ds is absolutely convergent and in H ∞,2 , and the limit U * ∞ is still an implementable symplectomorphism. Hence, and we want to prove that γ ′ and γ ′′ are unitarily equivalent in L 2 . The offdiagonal entries in (54) yield u * (γ ′ + 1/2)v + v T (γ ′ + 1/2)ū = 0 and as U is a symplectomorphism, we get from (43) that u is invertible and vū −1 = u * −1 v T . Thus, We can now use a known method to solve the Lyapunov (or Sylvester) equations: where we used that γ + 1/2 ≥ 1/2, so that there is no problem in handling the integrals. Hence v = 0, and, using (43) again, u is a unitary operator. And thus γ ′′ = u * γ ′ u which proves the result.
We now write the HFB equations in a form that is reminiscent of a Hamiltonian structure, and use it to give a direct proof of the conservation of the energy.

Proof. Eq. (46) is equivalent to
Hence, we can rewrite the expression of the energy in terms of φ t , u t , and v t as We then compute the derivatives of H γ ′ 0 : which are in fact (55), (56), (57) using the HFB equations. Hence, using first the chain rule, then (55), (56), and (57), We can now use that the evolution equation (45) on U t is equivalent to  Table 1. Correspondence between the notations of this article and some notations common in the physics literature [15].
which then vanishes since v T tū t = u * t v t for a symplectomorphism (see (43)), and the terms give real traces.

Relation with the HFB eigenvalue equations
In this section, we link our work with the HFB eigenvalue equations often encountered in the physics literature [15,14,32].
To be explicit, we give, in Table 1, the correspondence between the notations of this article and those of an article of Griffin [15]. We note that the setting in [15] is not exactly the same as ours, since the class of external potentials V that we consider excludes trapping potentials, and the solutions Φ(r) considered in [15] are time-independent. Moreover, we note that in this paper, we give rigorous proofs in the case of a two-body interaction potential v such that v 2 is relatively formbounded with respect to the Laplacian, which excludes potentials as singular as gδ; hence, the correspondence we establish in this section is only formal. Nevertheless, we believe that pointing out this relationship is useful.
Moreover, we note that in the physics literature (see e.g., [15, (23)]), the HFB eigenvalue equations are often investigated using a generalized eigenbasis decomposition (using vectors often denoted by u j , v j which play the same role as below), which we can relate to our approach in the following manner, based on our discussion from Section 3.
Let U t = ut vt utvt , and let γ ′ 0 ≥ 0 be a trace class operator as in Prop. 3.9, with the orthonormal decomposition γ ′ 0 = j≥0 N j |ζ j ζ j |. Let u j,t := u * t ζ j and v j,t := −v * t ζ j . Then (46) yields which yield [15, (25)] by evaluation on the diagonal: We now consider a pair interaction potential v = gδ. We assume that φ is independent of time and u j,t , v j,t have the simple form We also distinguish the quantities corresponding to v = gδ by the index gδ. Then (45) formally yields the HFB eigenvalue equations as presented in the work of Griffin [15, (23)]. Note that (60), (61), and (62) imply that γ t (x; x) and σ t (x; x) are time independent, since the phases simplify.
We conclude that the HFB eigenvalue equations are the stationary version of our equation (45). It amounts to finding eigenvalues and eigenvectors for the matrix ΛS in (45), which is a nonlinear problem since Λ depends on γ and σ (that is, on u, v and γ ′ 0 ). Furthermore, the decomposition in functions u j and v j corresponds to a "diagonalization" of the generalized one-particle density matrix Γ in the sense of Proposition 3.8.

Existence and Uniqueness of Solutions to the HFB Equations
We prove the global in time existence and uniqueness of mild solutions to the time-dependent Hartree-Fock-Bogoliubov equations in the H 1 -setting, and for a class of interaction potentials including the Coulomb potential. The proof is based on a standard fixed point argument (through an application of the Cauchy-Lipschitz and Picard-Lindelöf theorem in the proof of Theorem 5.1).
Let X be a Banach space, f ∈ C(X) a continuous function on X, and −iA the infinitesimal generator of a strongly continuous semigroup G(t) such that G(t) ≤ exp(κt), for all t ≥ 0 with a fixed κ ∈ R. We say that a continuous function ρ : [0, T ) → X is a mild solution of the problem if ρ t solves the fixed point equation in integral form (with the integral in Bochner's sense).
In what follows we will use the space X = X 1 . Moreover, A determines the linear part in the HFB equations and f the nonlinear part. (The explicit form of A and f is given in Eq. (66) and Eq. (67).) Theorem 5.1. Let d ≤ 3 and ρ 0 = (φ 0 , γ 0 , σ 0 ) ∈ X 1 . Assume that the potentials V and v are such that • V is infinitesimally form-bounded with respect to the Laplacian, • v is symmetric, v(x) = v(−x), and such that v 2 ≤ CM 2 for some C > 0.
Then the following hold: (1) Existence and uniqueness of a local mild solution: There exists a unique maximal solution to the HBF equations (26) to (28) in the mild sense, for some 0 < T ≤ ∞. (2) Existence and uniqueness of a local classical solution: If ρ 0 ∈ X 3 , then and ρ t satisfies the HBF equations (26) to (28) in the classical sense.
where ρ := (φ, γ, σ) ∈ X 2 . Then the linear part in the HFB equations is given by and the nonlinear part f := (f 1 , f 2 , f 3 ) by From Lemma 5.4, below, we obtain that f is continuously Fréchet differentiable in X 1 and therefore is locally Lipschitz, and from Lemma 5.3, we obtain that G(t) = exp(itA) defines a strongly continuous uniformly bounded semigroup on X 1 .

Consequently, we can rewrite the HFB equations (26) -(28) as a fixed point problem
and use the Banach contraction theorem to show that (26) -(28) have the unique local mild solution to in X 1 for the given initial data. (For the details for this standard argument, see [24], Section 9.2e, Theorem 3.) We will now prove our main Lemmata on G(t) = exp(itA) and f . First, we introduce norms used below: If V is infinitesimally form-bounded with respect to the Laplacian, and v is symmetric, v(x) = v(−x), and such that v 2 ≤ CM 2 for some C > 0, then G(t) = exp(itA) defines a strongly continuous, uniformly bounded semigroup on X 1 , i.e., the map t → G(t) B(X 1 ) is bounded.
Note that the proof Lemma 5.3 uses that −∆ is h-bounded, and h is −∆bounded. The latter condition does not hold for confining potentials. The proof also uses the translation invariance of M .
Proof. The letter C denotes a constant, which changes along the computations below. For (φ, γ, σ) ∈ X 1 , and thus exp(−ith)σ H 1 s ≤ C σ H 1 s which completes the proof.
The estimates on the operators b and k of Lemma E.1 allow us to control the nonlinear term f in the HFB equations.
Lemma 5.4. Assume that the pair interaction potential v is symmetric, v(x) = v(−x), and that v 2 ≤ CM 2 for some C > 0.
Proof of Lemma 5.4. Each f j is a linear combination of multi-linear maps. It is thus enough to prove continuity estimates for each of those multi-linear terms to prove that f is both well defined and Fréchet differentiable. It is sufficient to prove that, for the quadratic part , γ] + X 1 ≤ C ρ 2 X 1 , and, for the cubic part All the cubic estimates can be deduced from their quadratic counterparts using We thus only consider the quadratic terms. We estimate [b[γ], γ] using Lemma. E.1.(1) , γ] is controlled with the same method. Let us give the estimates for the first term b[γ]φ in full detail. Using Lemma E.1.(1) For the other terms we only give the main steps, using Lemma. E.1 each time. For k[σ]φ, using Lemma E.1. (2), In other words, ρ t , if it exists at all, then satisfies the differential equation (63) in the obvious sense.

Proof of Theorem 5.1.(3) [Conservation Laws].
For classical solutions, the conservation of the number particle and of the energy were proven as consequences of the same conservation laws for the many body system in Theorem 2.8 and 2.5. Another proof of the conservation law for the energy using only the HFB equations (independently from the many body problem) was given in Prop. 3.12, and the conservation of the particle number could also be proven directly from (27). We can now use those results since we proved the local existence of a classical solution. The conservation laws then extend to mild solutions by approximation.
Proof of Theorem 5.1.(4) [Global Solution]. We recall that for a maximal solution ρ t of the mild problem (64) defined on an interval [0, T ), we have that either T = ∞ or sup t∈[0,T ) ρ t X 1 = ∞ (see, e.g., Theorem 4.3.4 in [10]). It is thus enough to prove that sup to show that the solutions are global. Let Because V is infinitesimally form bounded with respect to the Laplacian, holds. And, because v 2 ≤ CM 2 , it follows that, for any ε ∈ (0, 1], (We write C for constants which depend on v, d and change along the estimates.) Applying this with ε = 1/(3(n − 1)), Then, summing the n(n − 1)/2 terms of this form on each n-particles subspace of the Fock space, we obtain that Hence, from the definition of H, (71) and (74) we get We now take the expectation value of ω q t . Using the conservation of the particle number and of the energy, Combined with the conservation of the particle number, this estimate provides bounds on γ t H 1 and φ t H 1 that are uniform in t. Moreover, uniform bounds on σ t H 1 s are then obtained from Proposition 3.1. It thus follows that the solution is global, as claimed.

Gibbs states and Bose-Einstein condensation
In this section, we determine translation-and U (1) gauge-invariant Gibbs states for the HFB equations without an external potential, and with an interaction potential gδ, and discuss the emergence of a Bose-Einstein condensate at positive temperature. (Recall from the introduction that U (1) gauge-invariant Gibbs states for the HFB equations are, in fact, Gibbs states for the Hartree-Fock equations.) We consider the system on a torus, Λ L = R d /2LZ d , i.e., [−L, L] d with periodic boundary conditions. Accordingly, we denote Λ * L := π L Z d the lattice reciprocal to 2LZ d . We will eventually take the thermodynamic limit, L → ∞, and discuss the emergence of a Bose-Einstein condensate.
The Hamiltonian H of the Bose gas is U (1) gauge-invariant (that is, invariant under the transformation ψ ♯ → (e iθ ψ) ♯ ), and, as we consider the case with no external potential, translation invariant. On a compact torus, where the volume is finite, these symmetries are also present in the Gibbs states of system (the notion of translation invariance should be, of course, appropriately modified). We are interested in quasifree states ω q L which on the one hand satisfy both the U (1) gauge invariance and the translation invariance, and, on the other hand satisfy a fixed point equation corresponding to the consistency condition (32) in the dynamical case: where β > 0 is the inverse temperature, µ is the chemical potential, and Ξ = Tr[exp(−β(H HF B (ω q L )) − µN)]. The U (1) gauge-invariance of ω q L then implies that the truncated expectations φ ω q L and σ ω q L vanish. Indeed, if one of them was nonzero, then the HFB Hamiltonian H HF B would include terms which would break U (1) gauge invariance, such as dx m(x) ψ * (x)ψ * (x) + h.c. . The quasifree states we consider are thus characterized by their truncated expectation γ L , and we will replace the variable ω q L by γ L in the sequel of this section. We use the expression of the HFB Hamiltonian (31) with v = gδ (and φ = 0, σ = 0), although this expression was derived for more regular interaction potentials v's: with n = n(x) = γ L (x; x). The translation invariance implies that the kernel γ L (x; y) is a function of x − y, that we still denote by γ L , and therefore n = n(x) = γ L (x; x) is independent of x.
Applying the fixed point equation (78) with A = ψ * (y)ψ(x) one can express it equivalently in the variable γ L : for n ∈ [0, ∞). The operator γ L is a pseudodifferential operator with symbol of γ L . Thus As the Fourier coefficients of γ L depend only of the number n, we obtain from (80), (81) and (82) a nonlinear fixed point equation for n: Note that the knowledge of n satisfying (83), or of γ L satisfying (80) or of ω q L satisfying (78) are equivalent.
From a physical point of view, it is natural to fix the density n, which can be tuned in an experiment and to compute µ. So n will be a parameter and we will solve (83) with the unknown µ. Lemma 6.1. Let g, β, n > 0, and, for d ≥ 3. Let n c be the critical density where ζ(x) = n≥1 n −x and Γ(x) = ∞ 0 t x−1 e −t dt. We define S L : (−∞, gn) → R and S ∞ : (−∞, gn] → R through Then: • There exists a unique µ L (n) < gn such that (83) holds, i.e., • If n < n c , there exists a unique µ ∞ (n) < gn such that We extend the function µ ∞ to (0, ∞) by setting µ ∞ (n) = gn for n ≥ n c .
Remark 6.2. The critical density n c can be explicitly computed.
Proof. In the discrete case, the existence follows from the intermediate value theorem because the map S L is continuous with limits 0 at −∞ and ∞ at gn. The map S L is strictly increasing and thus there exists a unique µ L (n) such that n = S L (µ L (n)).
In Theorem 6.3, we prove that the thermodynamic limit γ ∞ of the self-consistent equation (80) for γ L is well defined and exhibits the so called Bose-Einstein condensation. Theorem 6.3. Let g, β, n > 0 and d ≥ 3. Let γ L , n c , µ L and µ ∞ as defined in (80) and Lemmata 6.1. Then whereγ Remark 6.4. The presence of the δ(k) term is interpreted as the existence of Bose-Einstein condensation, because there is an accumulation of particles in the zero mode. It occurs when β d/2 n ≥ C d with C d a constant depending only on the dimension.
Proof of Theorem 6.3. First we prove the convergence of µ L (n) towards a µ ∞ (n). We first remark that µ L (n) ≥ −C for some constant C > 0 independent of L.
(Otherwise one could extract a subsequence such that n = S Lj (µ Lj (n)) → 0 < n.) Thus the accumulation points of µ L (n) are contained in [−C, gn]. Let µ Lj (n) denote an extracted sequence converging to an accumulation point µ ′ .
In the case n < n c : If µ ′ = gn then for j large enough µ Lj (n) ≥ (µ ∞ (n) + gn)/2, thus and which would lead to a contradiction. Note that it is crucial that µ ∞ (n) < gn for n < n c to get the convergence to the integral S ∞ gn+µ∞(n) 2 . It thus follows that µ ′ < gn. Then S Lj (µ Lj (n)) converges to n, because by definition of µ L (n) this sum is equal to n, and also to S ∞ (µ ′ ). (One has to control the dependency in µ Lj (n) in the Riemann sums.) Hence µ ′ = µ ∞ (n) and the unique accumulation point is µ ∞ (n). We thus proved the convergence of µ L (n) to µ ∞ (n).
In the case n ≥ n c , we sketch an argument similar to the one above. If an accumulation point µ ′ was such that µ ′ < gn, then the sums S Lj (µ Lj (n)) would converge to integrals with a value strictly smaller than n c and thus strictly smaller than n. This would lead to a contradiction. Thus the only possible accumulation point is gn and µ L (n) → gn = µ ∞ (n).
We now prove the convergence of γ L towards γ ∞ . Let ϕ ∈ C ∞ 0 (R d ). For L large enough the support of ϕ is included in Λ L , and On the other hand γ ∞ , ϕ D ′ = γ ∞ , ϕ S ′ = γ ∞ ,φ S ′ (Note that in the normalization we choose, the Fourier coefficients of ϕ on Λ L and the Fourier transform coincide, there is thus no need to specify the hat notation.) The convergence of γ L to γ ∞ is thus equivalent to for all ϕ.
In the case n < n c the convergence is thus just a convergence of Riemann sums of the integral (with the small additional difficulty that µ L (n) depends on L in the sum) because there is no singularity in the function k → (exp(β(k 2 +gn−µ ∞ (n)))− 1) −1 .
In the case n ≥ n c : Let ε > 0. First note that, for any fixed η > 0 as L → ∞. We choose η > 0 small enough so that The first condition on η yields then, the second condition on η implies lim sup Hence lim sup and as this holds for any ε > 0, we get the result.
Appendix A. Definition of quasifree states For brevity, we write ψ ♯ j := ψ ♯ (x j ). We recall that the truncated expectations are defined via where P n are partitions of the ordered set {1, ..., n} into ordered subsets.
In particular the equations of the dynamics for γ t and σ t can be rewritten We compute W * φt HW φt modulo terms of odd degree and of degree 0 in the creation and annihilation operators: Because ω q C,t vanishes on monomials of odd order in the fields and using the commutator, the knowledge of W * φt HW φt modulo terms of odd degree and of degree 0 in the creation and annihilation operators is sufficient to compute the time derivative (109) of γ t . Thus using the CCR we get i∂ t γ t (x; y) = ω q + v(z − x)φ t (z)φ t (x)ψ * (y)ψ * (z) − v(z − y)φ t (z)φ t (y)ψ(z)ψ(x) + v(z − x)ψ * (y)ψ * (z)ψ(x)ψ(z) − v(z − y)ψ * (z)ψ * (y)ψ(z)ψ(x) dz .
Appendix D. Self-adjointness of the Hamiltonian H The assumption v 2 ≤ C(1 − ∆) ensures that the original Hamiltonian H in (1) is self-adjoint, although not necessarily semi-bounded. In fact the same arguments as those used to prove (74) allow to deduce with V defined in (74), T defined in (70), from for some C j > 0. One can then use the KLMN theorem and the Nelson theorem (see [34,30]) to prove the self-adjointness of H. (Details can be adapted from, e.g., [1,Section 3].) Appendix E. Operators b and k Lemma E.1. Assume that the pair interaction potential v is symmetric, v(x) = v(−x), and such that v 2 ≤ CM 2 for some C > 0.
Then, the operators b and k defined in (29) and (30) possess the following properties: (1) b is linear and continuous from H 1 to B(H 1 ) ≃ M B(L 2 )M −1 .
(2) k is linear and continuous from H 1 s to M L 2 (L 2 )M −1 .
Proof. The linearity is clear. For the detailed proof of statement (1), we refer to [7]. For the reader's convenience, we recall here the main arguments.
For statement (1), we first consider the direct term, i.e., the first term in the definition of b. It is sufficient to prove that v * n (with n(x) = γ(x; x)) and ∇v * n are functions uniformly bounded by γ H 1 . As those two bounds are very similar, we focus on the more difficult one, ∇v * n. Note that γ can be decomposed as γ = ∞ j=1 λ j |ϕ j ϕ j | with λ j ≥ 0. Using this decomposition, combined with the Cauchy-Schwarz inequality, we find The estimates for the exchange term (the second term in the definition of B) are similar.
This yields the result as the last line is controlled by 4C σ 2