The Dirac-Frenkel Principle for Reduced Density Matrices, and the Bogoliubov-de-Gennes Equations

The derivation of effective evolution equations is central to the study of non-stationary quantum many-body sytems, and widely used in contexts such as superconductivity, nuclear physics, Bose-Einstein condensation and quantum chemistry. We reformulate the Dirac-Frenkel approximation principle in terms of reduced density matrices, and apply it to fermionic and bosonic many-body systems. We obtain the Bogoliubov-de-Gennes and Hartree-Fock-Bogoliubov equations, respectively. While we do not prove quantitative error estimates, our formulation does show that the approximation is optimal within the class of quasifree states. Furthermore, we prove well-posedness of the Bogoliubov-de-Gennes equations in energy space and discuss conserved quantities.

The time evolution of the state ψ t of a system of N quantum particles is described by the time-dependent Schrödinger equation where ψ t ∈ L 2 a (R 3N ) for spinless fermions (where the wave function is totally antisymmetric under any exchange of particles) and ψ t ∈ L 2 s (R 3N ) for bosons (where the wave function is totally symmetric under any exchange of particles). In this generality, the Schrödinger equation models a vast range of physical systems, starting from nucleons in the atomic nucleus over electrons in semiconductors to stars for the fermionic theory, or Bose-Einstein condensates for the bosonic theory, depending on the choice of the one-particle Hamiltonian h (we think of h = −∆ + V ext (x) with some external potential V ext : R 3 → R; h i denotes this operator as acting in the variable x i ∈ R 3 ) and the pair interaction V : R 3 → R. Unfortunately, these systems also have an enormous number of degrees of freedom, making analytical and numerical solutions generally impossible. For this reason, there is a lot of interest in approximate theories (also called effective evolution equations), which contain fewer degrees of freedom and make analytical and numerical treatments possible. Of course, such theories do not achieve the broad validity of the Schrödinger equation and provide a good approximation only in specific physical regimes, which are mathematically modeled as scaling limits. In this paper, we discuss a geometric method for the derivation of effective evolution equations. This method, even though not the most convenient for proving quantitative error estimates for the obtained approximation, directly shows that the obtained equations are optimal as far as the available degrees of freedom permit. The method is also independent of any choice of scaling limit.
For this introduction we focus on (for simplicity of notation spinless) fermionic systems. The corresponding bosonic notions will be introduced in Sect. 4. Here we deal only with pure quasifree states. Only in Sect. 5 for the topic of well-posedness we also consider mixed states; we will always highlight explicitly when we talk about mixed states.
The most basic approximate theory of fermionic systems is obtained by restricting the Schrödinger equation to wave functions that are Slater determinants, also denoted as the antisymmetrized tensor product ψ = (N !) −1/2 f 1 ∧· · ·∧f N , where the one-particle wave functions f j constitute an orthonormal set in L 2 (R 3 ). The corresponding effective evolution equation for the Slater determinant is given by the Hartree-Fock system of N nonlinear coupled PDEs for the one-particle wave functions: (1.2) More conveniently, the Hartree-Fock equations can be formulated in terms of the oneparticle reduced density matrix. The one-particle reduced density matrix of a state ψ ∈ L 2 (R 3N ) is defined as the non-negative trace-class operator γ on L 2 (R 3 ) obtained by taking the partial trace over N − 1 particles of the many-body density matrix |ψ ψ|, γ = N tr 2,...N |ψ ψ|, where we have chosen to normalize the one-particle reduced density matrix such that tr γ = N . If ψ = (N !) −1/2 f 1 ∧ . . . ∧ f N , we find the rank-N projection γ = N j=1 |f j f j |. Conversely, every rank-N orthogonal projection specifies (uniquely up to a phase) a Slater determinant (just take the spectral decomposition of the projection to find the one-particle wave functions f j ).
If the one-particle wave functions have time dependence given by the Hartree-Fock equations (1.2), then γ t = N j=1 |f j,t f j,t | satisfies the equivalent equation where V * ρ γt (x) = dy V (x − y)γ t (y, y) is a multiplication operator called the direct term, and X V (γ t )(x, y) = V (x − y)γ t (x, y) is the integral kernel of an operator called the exchange term.
They satisfy the canonical anti-commutation relations (CAR), i. e. The vector Ω = (1, 0, 0, . . .) ∈ F a is called the vacuum state and is in the kernel of all annihilation operators, a(f )Ω = 0 for all f ∈ L 2 (R 3 ); it describes a system not containing any particles. It is convenient to introduce the operator-valued distributions a * x and a x with the defining property that (in a weak sense, within expectation values) They satisfy the formal canonical anti-commutation relations {a x , a y } = 0, {a * x , a * y } = 0, and {a x , a * y } = δ(x − y). The Hamiltonian is generalized to Fock space as , (1.4) where H n denotes the first quantized Hamiltonian as given in (1.1). The Hamiltonian H can also be represented in terms of creation and annihilation operators by 1 H = dΓ(h) + 1 2 dxdy V (x − y)a * x a * y a y a x . (1.5) Restricted to N -particle states, this Hamiltonian agrees with H N . For the geometric considerations in this paper, the explicit form of H does not need to be specified, as long as it is a self-adjoint operator and conserves the number of particles. (A Hamiltonian conserves the number of particles if it commutes with the particle number operator N defined by N ψ = nψ (n) ∞ n=0 . In particular all operators of the form (1.4) conserve the particle number.) Only Sect. 5 refers to the particular Hamiltonian (1.5).
Using the creation and annihilation operators, the definition of the one-particle reduced density matrix γ is extended to states ψ ∈ F a by defining it to have the integral kernel γ(x, y) := ψ, a * y a x ψ . (1.6) We can now give a simple proof that γ ≤ 1: For all f ∈ L 2 (R 3 ), using the CAR, Slater determinants are a special case of a class of more general states in Fock space called quasifree states (see e. g. [57] for a very readable introduction). The defining property of quasifree states is that they are exactly those states ψ ∈ F a for which the Wick theorem holds, i. e. expectation values of creation and annihilation operators can be reduced to the sum of the expectation values of all possible pairings of just two operators (with the sign being the sign of the corresponding pairing); for example (Here we used the notation a ♮ j to denote an operator without specifying whether it is a creation or annihilation operator.) Notice that the Wick theorem allows us to express any expectation value of creation and annihilation operators in a quasifree state purely in terms of the one-particle reduced density matrix γ and the pairing density α(x, y) := ψ, a y a x ψ . (1.7) Slater determinants are exactly those quasifree states for which the pairing density identically vanishes, α = 0. Notice that for a general quasifree state γ is not a rank-N projection; instead a quasifree state always satisfies γ 2 − γ = αα and αγ = γα. (1.8) Given γ and α satisfying (1.8), there is a (up to a phase unique) quasifree state ψ ∈ F a such that (1.6) and (1.7) hold. There is a more compact way of writing the equations (1.8) by introducing the generalized one-particle reduced density matrix Γ. The generalized one-particle reduced density matrix Γ is an operator on L 2 (R 3 ) ⊕ L 2 (R 3 ) given by The characterization (1.8) of quasifree states is equivalent to Γ being an orthogonal projection, Γ 2 = Γ = Γ * . Any (not necessarily quasifree) generalized one-particle reduced density matrix has the property 0 ≤ Γ ≤ 1, and thus Γ 2 ≤ Γ. (1.10) In the theory of superconductivity, the pairing density is interpreted as the wave function of electrons that have formed Cooper pairs, which in many ways behave like bosons. These Cooper pairs are seen as the carriers of the superconducting current that has attracted so much attention for its technological applicability in the dissipationless transport of electricity.
In this paper, our focus lies on the effective evolution equation obtained by restriction of the many-body evolution to quasifree states with pairing. This system of effective evolution equations is known as the Bogoliubov-de Gennes equations ). More compactly, γ t and α t satisfy (1.11) if and only if the generalized one-particle density matrix (1.12) as in [37,36] we use the generalized Hartree-Fock operator In the present paper, our goal is to formulate a systematic approximation principle by which we can obtain the Bogoliubov-de Gennes equations from many-body quantum theory. The approximation principle we establish is a reformulation of the Dirac-Frenkel principle in the space of reduced density matrices, and it yields the equations sometimes called the quasifree reduction principle. Applying the quasifree reduction principle to the Hamiltonian (1.5), one obtains the Bogoliubov-de Gennes equations. Afterward we study well-posedness and conserved quantities for the Bogoliubov-de Gennes equations.
While it is general knowledge that the quasifree reduction principle should be a consequence of the Dirac-Frenkel principle, we are not aware of a direct proof having appeared before; in particular the formulation of the Dirac-Frenkel principle in terms of reduced density matrices has not been given before. Among the advantages of our approach is that it shows that the obtained effective equations describe the optimal evolution possible within the approximation manifold.
Earlier results. The derivation of effective evolution equations for many-body systems has attracted a lot of attention in the community of mathematical physics and can be seen as a cornerstone of non-equilibrium statistical mechanics. Consequently, the literature is vast and we cannot claim to provide a complete overview. Let us say so much, that the geometric approximation principle on which we build in this paper goes back to the founding fathers of quantum mechanics [18,27]. A rigorous mathematical discussion and highly valuable presentation has been given in [42].
The next step after the geometric derivation of the correct effective equation lies in the proof of convergence toward the effective equation and then the derivation of quantitative error bounds in given physical regimes modeled as scaling limits. This topic has attracted a lot of attention in recent years. For bosonic systems, many such results on the approximation of reduced density matrices have been proven for the meanfield model [6,26,53,16,48,40,56,1,39,29,2] and the Gross-Pitaevskii model [7,21,23,22,24,25,49,47]. For fermionic systems, most results have only appeared in the last few years. The main regimes treated here are the mean-field regime on short time scales [5,28], the mean-field regime with slow variation of the effective interaction potential [4,46,45], and the combined mean-field/semi-classical limit for high-density systems [44,58,20,10,9,8,50]. The papers cited in this paragraph do not use the Dirac-Frenkel principle but other methods that have been specifically developed for many-body systems, like the BBGKY hierarchy, coherent states, a Schwinger-Dyson expansion or counting the number of particles well-described by the effective evolution equation. Some applications of the Dirac-Frenkel principle with explicit error estimates can be found in [42]; another example is [30].
In the context of proving convergence toward effective equations for bosonic systems in mean-field and Gross-Pitaevskii scaling limits, equations of the Hartree-Fock-Bogoliubov-type appear when considering second-order corrections beyond the above results on approximation of reduced density matrices, that is approximations in the norm of the many-body Hilbert space [32,31,34,33,43,11]. The difference to our discussion here is that we are just interested in the approximation of reduced density matrices, not in the norm of the many-body Hilbert space, but instead we focus on ensuring optimality of the derived effective equations. Moreover, while the scaling regimes are crucial for obtaining convergence and quantitative error estimates, in the qualitative geometric approach that we use, it is not necessary to specify a particular scaling limit.
For fermionic systems, the derivation of quantitative error bounds in appropriate scaling limits for the Bogoliubov-de Gennes equations or even the BCS equations remains an open problem.
Organization of the paper. In Sect. 2.1 we shall recall the variational principle of Dirac and Frenkel, which is at the base of our paper. In Sect. 2.2 we recall the principle of quasifree reduction, which is less fundamental and less general than the Dirac-Frenkel principle, but more convenient for explicit calculations. In Sect. 3.1 we give a re-derivation of the Hartree-Fock equation as the simplest example to introduce our formulation of the Dirac-Frenkel principle in the space of one-particle reduced density matrices. In Sect. 3.2 we use our formulation of the Dirac-Frenkel principle for systems with pairing, yielding the time-dependent Bogoliubov-de Gennes equations. In Sect. 4, we present the analogous formulation of the Dirac-Frenkel principle for bosonic systems. Finally, in Sect. 5 we discuss well-posedness of the time-dependent fermionic Bogoliubov-de Gennes equations.

Approximation Principles
In this section we first introduce the Dirac-Frenkel principle, which can be seen as the fundamental principle for deriving optimal effective evolution equations. Afterwards we introduce the principle of quasifree reduction which does not exhibit the optimality but leads to the same results as the Dirac-Frenkel principle and is calculationally simpler.

The Dirac-Frenkel Variational Principle
In this section we discuss the abstract formulation of the Dirac-Frenkel variational principle for the approximation of the time-dependent Schrödinger equation by projection onto a submanifold. The Dirac-Frenkel variational principle is particularly interesting for its clear geometrical content, which shows that the obtained equation on the submanifold is the optimal choice. We follow [42,Chapter II].
We consider the Schrödinger equation as an evolution equation in a complex Hilbert space H. We use the convention that the scalar product is anti-linear in the first and linear in the second argument. The Hamiltonian H is a self-adjoint operator on H. The Schrödinger equation reads (2.14) Now consider a smooth (typically infinite dimensional) submanifold M of H. The tangent space of M in the point u ∈ M is denoted by T u M; it consists of the derivatives u ′ 0 of all differentiable paths t → u t passing through u 0 = u. We are interested in approximating the solution ψ t of the Schrödinger equation with a path u t on the manifold M, assuming that initially ψ 0 = u 0 ∈ M. As pointed out already by Dirac, the path t → u t is to be chosen such that at every time t the derivative u ′ t ∈ T ut M is as close as possible to 1 i Hu t ; in other words, the path is determined by choosing its derivative as the orthogonal projection of 1 i Hu t onto the tangent space:

15)
P (u t ) being the orthogonal projection from H onto T ut M. From this formulation it is clear that the Dirac-Frenkel principle yields the effective evolution equation which at every infinitesimal time step is optimal.
In the case of fermionic many-body systems, one typically chooses H = L 2 a (R 3N ) (for the case of no pairing) and M as the set of N -particle Slater determinants. While this approach does yield the time-dependent Hartree-Fock equations (1.2), it is not expected to ever do so with controllable errors: it is a general fact that in many-body systems, the norm of many-body wave functions as a measure of distance has unfortunate behavior as the number of particles grows. In fact, the quantitative derivation of effective equations for many-body systems in appropriate scaling limits is typically proven in terms of the trace norm or Hilbert-Schmidt norm of reduced density matrices (see the overview of results in Sect. 1). For this reason, and to make the connection to the principle of quasifree reduction, in this paper we formulate the Dirac-Frenkel principle in the space of one-particle density matrices. This formulation will be calculationally very convenient when we consider quasifree states with pairing.

The Principle of Quasifree Reduction
The principle of quasifree reduction appears to be the computationally most accessible principle for deriving effective equations for many-body quantum systems. It applies to the particular case were the approximation manifold is given by a class of quasifree states. Typically it is formulated directly in the language of reduced density matrices.
The principle of quasifree reduction asserts that for fermionic systems the effective evolution equations are where ψ qf t ∈ F a is the quasifree state uniquely (up to a phase) assigned to γ t and α t . For bosonic systems the equations proposed as the quasifree reduction principle, including a condensate ϕ t ∈ L 2 (R 3 ), are the following [3] where ψ bog t in the bosonic Fock space is the Bogoliubov state associated with (ϕ t , γ t , α t ) (see (4.38) for the definition).
While this principle is easy to formulate, calculationally accessible and has been frequently used in many contexts (e. g. very recently to derive the Hartree-Fock-Bogoliubov equations for bosonic systems [3] and the Bogoliubov-de Gennes equations for fermionic systems [17]), it is not completely obvious that it is a consequence of the more fundamental Dirac-Frenkel principle. Maybe more severely, it is not at all clear whether the quasifree reduction principle yields the optimal approximation possible within the manifold of quasifree states.
In the present paper, we prove that the principle of quasifree reduction does follow directly from the Dirac-Frenkel principle, in particular showing that it yields the optimal effective evolution equations.

Derivation of the Quasifree Reduction Principle for Fermions
In this section we derive the principle of quasifree reduction from a re-formulation of the Dirac-Frenkel principle. We first sketch the instructive case of systems without pairing before generalizing to fermionic systems that also exhibit pairing.

Fermionic Systems without Pairing
Here we shall warm up with the case of no pairing (α = 0), i. e. giving a derivation of the standard Hartree-Fock equation. To this end we shall formulate the Dirac-Frenkel principle in terms of the one-particle reduced density matrix and from there, derive the principle of quasifree reduction. To our knowledge, this is the first formulation of the Dirac-Frenkel principle in terms of reduced density matrices. Notice that the derivation of the quasifree reduction principle does not make use of any particular form of the Hamiltonian; we assume only that it commutes with the particle number operator (i. e. the number of particles is conserved along the many-body evolution).
We consider the Schrödinger equation in the space L 2 a (R 3N ), i. e. describing a fermionic system of N -particles. The many-body evolution t → ψ t induces an evolution of the associated one-particle reduced density matrix γ ψt , which satisfies (3.18) The one-particle reduced density matrix is a non-negative operator on L 2 (R 3 ). Since our system has a finite number of particles, tr γ ψt = N , the one-particle reduced density matrix is a Hilbert-Schmidt operator. We thus choose the ambient Hilbert space in which γ ψt lives to be Due to the condition of self-adjointness this is only a real-linear (instead of complexlinear) space; in the following all spaces are real-linear only. Corresponding to Slater determinants in the wave function picture, we choose our approximation manifold to be given by orthogonal projections, This is an infinite dimensional Hilbert submanifold of the self-adjoint Hilbert-Schmidt operators. The effective evolution equation is to be found in M. This is achieved in the optimal way by applying the Dirac-Frenkel principle reformulated in terms of the one-particle reduced density matrix.
Dirac-Frenkel Principle for Reduced Density Matrices. The effective evolution equation within the submanifold M is given by where ψ qf t is the state uniquely (up to the phase) associated with γ t , and proj(γ t ) : H → T γt M is the projection onto the tangent space in the point γ t .
We start by determining the tangent space and the projection onto the tangent space.
The orthogonal projection from H onto T γ M is given by Proof. Let A ∈ T γ M. By definition there exists a differentiable curve t → γ t in M such that γ 0 = γ and γ ′ 0 = A. (By definition of differentiability in the norm of the ambient Hilbert space H, A is a Hilbert-Schmidt operator.) Taking the derivative of the projection condition γ 2 t = γ t , we find Aγ + γA = A. Multiplying from the left and right by γ, we get 2γAγ = γAγ, so γAγ = 0. Furthermore, by multiplying it from the left and the right by ( is anti-self-adjoint and Hilbert-Schmidt, so e tB is a unitary. Now let γ t := e tB γe −tB . This is a curve of orthogonal projections with γ 0 = γ. Its derivative at zero is It is easy to see that proj(γ) is a projection and has the claimed image.
The manifold of orthogonal projections M has several connected components, corresponding to the value tr γ ∈ N. Clearly, any differentiable curve always stays within the same connected component, so we do not have to worry about this. , with H and M as chosen above, is equivalent to the quasifree reduction principle without pairing (α = 0): where ψ qf t is the Slater determinant uniquely (up to the phase) associated with γ t . (Equation (3.20) is the same as equation (2.16) but written without the use of operator-valued distributions.) We refer the reader to the more general proof of Theorem 3.4.
One may convince oneself that (3.20) indeed yields the Hartree-Fock equation (1.3) when evaluating the expectation value on the r. h. s. with the many-body Hamiltonian from (1.5), using the canonical anti-commutation relations and the Wick theorem.
Having (1.3) at hand, we can also obtain the quasifree reduction principle from the Dirac-Frenkel principle for reduced densities as follows. Clearly (1.3) implies that γ 2 t also satisfies the Hartree-Fock equation, So if we have a projection as initial data, γ 2 0 = γ 0 , assuming uniqueness, we conclude that γ 2 t = γ t for all times t. Alternatively, we could argue that the Hartree-Fock equation preserves the spectrum of γ t , which also implies γ 2 t = γ t for all times. Either way, we conclude that the derivative is in the tangent space of M, which makes the projection in the Dirac-Frenkel principle trivial and yields the quasifree reduction principle. However: • This argument uses (1.3) which is obtained by explicitly evaluating the quasifree reduction principle. Using only the equations of the quasifree reduction principle (2.16), there is no easy way to formulate (3.21); in fact, a direct verification that (2.16) stays within M seems complicated to us.
• The argument depends on the choice of the one-particle Hamiltonian h and regularity and decay of the interaction potential V and the initial data. For the initial value problem with pairing (α 0 = 0), uniqueness or conservation of the spectrum are by themselves non-trivial problems, see Sect. 5.
Our derivation does not require any specification of the Hamiltonian beyond its existence as a self-adjoint, particle number conserving, operator. Furthermore, our geometric approach makes it clear that the quasifree reduction principle is the optimal approximation within the set of quasifree states.
or ∂ t γ t = 0 would be an evolution in the manifold of quasifree states-but far from being the optimal approximation to the many-body problem.)

Fermionic Systems with Pairing
We now extend our formulation of the Dirac-Frenkel principle to derive the approximation of fermionic many-body systems by quasifree states with pairing, namely the Bogoliubov-de Gennes equations. As before, the derivation of the quasifree reduction principle does not require any particular form of the Hamiltonian; our only assumption is that it conserves the number of particles.
The geometry becomes very similar to the case of no pairing by using the generalized one-particle reduced density matrix Γ. Thus the manifold of quasifree states with pairing can be described by Γ 2 = Γ and the block structure (1.9), which however comes 'for free' since it is present in all generalized one-particle reduced density matrices (in particular also in the one derived from the many-body Schrödinger equation).
Let us be a bit more precise and define the involved spaces. First of all notice that, due to the form (1.9), we have the generalized one-particle density matrix is not a Hilbert-Schmidt operator. We remedy this problem by considering the generalized one-particle reduced density matrix as a point in an affine space, and the approximation manifold as a submanifold of this affine space. Let us denote by Γ vac := 0 0 0 1 the generalized one-particle reduced density matrix of the vacuum Ω ∈ F. Then any generalized one-particle reduced density matrix can be written as Every generalized one-particle reduced density matrix satisfies Γ 2 ≤ Γ, which implies γ 2 − αα ≤ γ (only for quasifree states we had equality here); thus which is twice the expected number of particles and as such assumed to be finite. The expected number of particles is trivially conserved along the many-body evolution since we assume the Hamiltonian to commute with the particle number operator; it is typically also conserved along the effective evolution, c. f. Lemma 5.4, so it is justified to take Γ as a Hilbert-Schmidt operator. Let us therefore introduce the affine space Similar to the no-pairing case, A is a real-linear space.
Now notice that the requirement of having the block structure of Γ in terms of γ and α as in (1.9) can be rewritten 2 as the condition being the anti-linear operator of complex conjugation. So we can think of the evolution of the many-body generalized one-particle reduced density matrix as living in the affine subspace of A given by (But not every Γ ∈ A − is the generalized one-particle reduced density matrix of a Fock space vector.) The approximation manifold is again given by the generalized one-particle density matrices corresponding to quasifree states: So compared to the no-pairing case not much has changed-the only additional complication is that we have to impose the block structure of Γ in terms of γ and α. Luckily, this block structure is present in any generalized one-particle reduced density matrix including the one of the many-body evolution. So the many-body evolution describes a curve in the affine subspace A − , of which M is a submanifold.
To provide a characterization of the tangent space, we also introduce as an auxiliary space the manifold of projections which do not necessarily have the block structure Notice that, since A is an affine space, T Γ A = A for any Γ ∈ A. (i) onto the tangent space of projection operators (3.25) (ii) onto the tangent space of the affine subspace with the block structure There is a subtlety here: The latter one however is not a Hilbert-Schmidt perturbation of Γvac and thus not a solution within A; in fact it corresponds to a state formally obtained from Γ by a particle-hole transformation replacing the vacuum by an infinite number of fermions filling up all the Hilbert space.
(iii) and onto the tangent space of quasifree states proj(Γ) : Then, for Γ ∈ M, we have Proof. The projection onto the tangent space of projection operators (3.25) is known from Lemma 3.1.
Since A − is an affine subspace, we can simply take the derivative of the defining equation to find It is easy to check that the formula (3.26) defines an orthogonal projection, maps into . Thus, using Lemma 3.1, Obviously P is an orthogonal projection. It is easy to verify that it maps into T Γ M and is surjective onto We can now derive the quasifree reduction principle from the Dirac-Frenkel principle.
Theorem 3.4 (Derivation of the Quasifree Reduction Principle, With Pairing). The effective equation for the generalized one-particle reduced density matrix Γ t obtained by applying the Dirac-Frenkel principle to the many-body evolution with M and A as chosen above yields the principle of quasifree reduction where ψ qf t is the quasifree state uniquely (up to its phase) assigned to Γ t . (Equation (3.28) is a compact way of writing (2.16), avoiding the use of operator-valued distributions by testing against F 1 and F 2 .) Proof. The proof uses some theory of Bogoliubov transformations, for which we recommend [57, Chapters 9 and 10] as a reference. (For the no-pairing case, the Bogoliubov transformation is a simple particle-hole transformations, see, e. g., [10,9].) Recall that the generalized one-particle reduced density matrix Γ ψ of a Fock space vector ψ is, avoiding the use of operator-valued distributions by testing against where the generalized creation and annihilation operators are and So for ψ t being the solution of the many-body Schrödinger equation, the associated generalized one-particle reduced density matrix satisfies Notice that the r. h. s., like the derivative of any differentiable curve of generalized oneparticle reduced density matrices, lies in T Γ A − . According to the Dirac-Frenkel principle, we have to project it onto the tangent space of quasifree states. We apply the projection as given by (3.27) and (3.25) to get where ψ qf t is the quasifree state uniquely assigned to Γ t . Comparing to the quasifree reduction principle (3.28), we see that we simply have to show that Since ψ qf t is a quasifree state, it can be written in terms of an implementable Bogoliubov Using this last identity we calculate that (Here we made use of the fact that g 2 ,g 1 as a complex number commutes with everything, and of the fact that any annihilation operator applied to the vacuum gives zero.) Similarly, we find for the other diagonal block as well that it vanishes, Using the Wick theorem and the CAR, it is a simple calculation that the quasifree reduction principle (3.28), applied to the Hamiltonian (1.5), yields the time-dependent Bogoliubov-de Gennes equations (1.11).
Remark. The reader may wonder how it is possible that the many-body evolution gives rise to the equation and the effective evolution solves the seemingly identical equation yet the two evolutions in general differ even if they both start from quasifree initial data. The answer is that (3.30) is not a well-posed initial value problem, simply because a general Fock space state has many more degrees of freedom than just the generalized oneparticle reduced density matrix; the r. h. s. is not a function of only Γ t . The equation (3.30) only makes sense if the r. h. s. is already prescribed by the Schrödinger equation On the other hand, (3.31) is a well-defined initial value problem because quasifree states in Fock space are (up to a phase) one-to-one with their generalized one-particle reduced density matrix. So the r. h. s. is a function only of Γ t here (alternatively think of the Wick rule which also shows that the r. h. s. can be expressed in terms of only Γ t ).
We provide the rigorous proof of well-posedness for a main class of physically relevant Hamiltonians and initial data in Sect. 5.

Derivation of the Quasifree Reduction Principle for Bosons
In this section we present the formulation of the Dirac-Frenkel principle for one-particle reduced density matrices of bosonic systems. This is slightly more complicated than for fermionic systems because the simple projection condition has to be replaced, and because we include a condensate, but can be treated by modifications of the previously developed geometric notions.
We start by reviewing some definitions for bosonic systems where they differ from the corresponding fermionic formulas. For a comprehensive introduction we refer to [57]. Bosonic Fock space is defined in the same way as for fermionic systems, simply replacing antisymmetric by symmetric wave functions: Creation and annihilation operators a * (f ) and a(f ) (where f ∈ L 2 (R 3 ), a one-particle wave function) are defined as The bosonic creation and annihilation operators satisfy the canonical commutation relations (CCR), i. e.
[ . Quasifree states are defined as those states for which the Wick theorem holds, which only differs from the fermionic case by having all positive signs, e. g., The r. h. s. can be expressed in terms of γ(x, y) = ψ, a * y a x ψ and α(x, y) = ψ, a y a x ψ .
For any bosonic quasifree state, γ and α are related by conversely all γ and α satisfying these two equations define a (up to a phase) unique quasifree state in bosonic Fock space. The generalized one-particle reduced density matrix is defined as The relations (4.32) characterizing it as belonging to a quasifree state can be rewritten The generalized one-particle reduced density matrix Γ is a non-negative operator on As for fermionic systems, also bosonic quasifree pure states can be written in terms of a Bogoliubov transformation [57]: If ψ ∈ F s is quasifree, then there exists an implementable Bogoliubov map V : where the generalized creation/annihilation operators A(F ), A * (F ), F ∈ L 2 (R 3 ) ⊕ L 2 (R 3 ) are defined exactly the same way as for fermions.
With Γ vac = 0 0 0 1 the generalized one-particle reduced density matrix of the vacuum (identical to the fermionic case), we define the spaces (4.37) These take the role of: A the ambient affine space defining the scalar product, A + the affine subspace in which the many-body evolution can be found, and M the approximation manifold of generalized one-particle reduced density matrices of quasifree states. To see that M is indeed a submanifold of A + , notice that by (4.32) (or by (2Γ + S)S(2Γ + S) = S, which is equivalent to (4.34)), we can write every Γ ∈ M as and thus M as a graph. Alas! The computation of the tangent spaces of M from its graph representation involves the derivative of the operator square root around 1, which leads to Lyapunov equations of type {X, A} = B with A = √ 1 + αα and B given. There is no simple closed formula for the solution to this equation in operator form (one can only express X as a function of the eigenvalues and the eigenfunctions of A). We overcome this problem by noticing that it is sufficient to have a parametrization of the orthogonal complement of the tangent space.
Proof. Let us define the map ψ : Γ → P = −ΓS. We can explicitly write down its inverse: Γ = −P S since S 2 = 1. Let us specify domains and codomains. In parallel to the spaces A, A + and M we introduce (notice that P vac := ψ(Γ vac ) = Γ vac ) It is easy to check that ψ is an isometric isomorphism A → G (both sides with the Hilbert-Schmidt scalar product) and is also an isometric isomorphism A + → G + . Furthermore, it is a diffeomorphism M → M G .
Following the strategy of Lemma 3.1 and Lemma 3.3, with the condition SP * S = P taking the place of self-adjointness everywhere, we obtain T P M G = B ∈ G : P BP = 0 = (1 − P )B(1 − P ) and B + J BJ = 0 .
The differential of ψ is given by D Γ ψB = −BS, which is also an isometric isomorphism T Γ A + → T P G + . In particular it conserves orthogonality, so it is also an isomorphism where the orthogonal complement is defined by the decomposition T P G + = T P M G ⊕ (T P M G ) ⊥ . Rewriting Noticing that T P G + = {B ∈ G : B + J BJ = 0}, the proof is complete.
Unlike fermionic states, bosonic states can exhibit condensation, so that for ψ ∈ F s it is possible that for some f ∈ L 2 (R 3 ) we have the additional degree of freedom ψ, a(f )ψ = 0.
(For any quasifree state this is vanishing.) Let us define the Weyl operator The Weyl operator is unitary and W(ϕ) * = W(−ϕ); furthermore they also satisfy W(ϕ 1 )W(ϕ 2 ) = W(ϕ 1 + ϕ 2 )e −i Im ϕ 1 ,ϕ 2 for all ϕ 1 , ϕ 2 ∈ L 2 (R 3 ). An ideal condensate is described by a coherent state Ψ = W(ϕ)Ω; we have and consequently Ψ, a(f )Ψ = f, ϕ . The expected number of particles in the coherent state is Ψ, N Ψ = ϕ 2 L 2 . Using the BCH formula 3 together with the CCR we find From the last formula we see that a coherent state is a linear combination of different particle numbers, where the probability to measure n particles is given by a Poisson distribution peaked at the value ϕ 2 L 2 . We now enlarge the class of quasifree states to the class of Bogoliubov states by including a condensate; more precisely, a Bogoliubov state 4 is any state of the form where ϕ ∈ L 2 (R 3 ) (typically not normalized) and V is any implementable Bogoliubov map. Using the fact that any expectation value of an odd number of creation and annihilation operators in a quasifree state vanishes, Furthermore we find that the one-particle reduced density matrix is given by (4.40) Similarly, we find the pairing density to be In other words, γ =γ + |ϕ ϕ| and α =α + ϕ ⊗ ϕ. Theγ andα so introduced are called the truncated expectations. They clearly satisfy the quasifree-propertỹ So by first obtaining ϕ through (4.39) and then solving (4.40) and (4.41) forα,γ, we have a natural way of assigning a unique (ϕ,γ,α) to every quasifree state; conversely every triple (ϕ,γ,α) satisfying (4.42) defines a (up to a phase) unique Bogoliubov state in Fock space through (4.38). So as we just argued, Bogoliubov states are characterized by independently the condensate wave function ϕ ∈ L 2 (R 3 ) and the truncated expectations, i. e.Γ. We therefore introduce the manifold where M is the manifold of quasifree generalized one-particle reduced density matrices as determined before. Of course, the tangent space is given by (4.43) So we can now formulate the Dirac-Frenkel principle for the condensate wave function and the generalized reduced density matrix of bosonic Bogoliubov states: Calculate the derivative of the condensate wave function evolving by the many-body Hamiltonian H in the Bogoliubov state associated with ϕ t andΓ qf t , then apply the projection onto the tangent space to ∂ t ϕ t . Calculate the derivative of the generalized one-particle reduced density matrix evolving by the many-body Hamiltonian H in the quasifree state associated withΓ qf t , then apply the projection onto the tangent space to ∂ tΓ qf t . The projected derivatives describe the effective evolution.
Theorem 4.2 (The Quasifree Reduction Principle for Bogoliubov States). The Dirac-Frenkel principle, applied by projecting the curve of the many-body evolution from the space L 2 (R 3 )×A to the approximation manifold of Bogoliubov states M bog = L 2 (R 3 )×M yields the equations of the quasifree reduction principle (2.17).
As it was already the case for fermionic systems, the only difference between the Dirac-Frenkel principle and the principle of quasifree reduction is the projection onto the tangent space. So instead of really doing the projection onto the tangent space, we simply check that the right-hand sides of (2.17) are orthogonal to the orthogonal complement of the tangent space.
Proof. Recall (4.43): as far as ϕ is concerned, the tangent space is given by all of L 2 (R 3 ); i. e. the projection onto the tangent space is just the identity. Therefore we only have to take care of projecting the evolution of γ and α; more precisely we will check that the derivative ofΓ qf t already lives in the tangent space. For this, it is sufficient to show that A, ∂ tΓ qf t S 2 = 0 for all operators A ∈ (T Γ M) ⊥ . Notice that a priori ∂ tΓ qf t ∈ T Γ A + , so it is sufficient to consider the orthogonal complement as a subspace of T Γ A + instead of all of A. So by Lemma 4.1, we can write A = − (P * BP * + (1 − P * )B(1 − P * )) S for some operator B.
Since B is Hilbert-Schmidt, it has a singular value decomposition B = j λ j |ξ j ϕ j |, So it suffices that every such expectation value vanishes individually. Recall that we have . Recall also that we can write ψ qf t = U V Ω, where the Bogoliubov map V satisfies Using these facts we find We now check that the second kind of expectation value vanishes. For all A similar calculation shows that also F 1 , P (∂ tΓ qf t )SP F 2 = 0. So we have shown that ∂ tΓ qf t ∈ TΓqf t M, thus the projection on the tangent space is trivial, and the Dirac-Frenkel principle becomes the quasifree reduction principle.
The calculation to obtain the explicit evolution equations from the bosonic principle of quasifree reduction (2.17) was sketched in [3], based on the canonical commutation relations and the Wick theorem. The explicit equations are known as the Hartree-Fock-Bogoliubov equations, here written in terms of the truncated expectationsα t andγ t and the condensate wave function ϕ t , where h HFB differs only by the sign of the exchange term from the fermionic h HF : More compactly, (4.45) can be written symplectically [3, Eq. (41)] whereΓ t is the truncated generalized one-particle density matrix (4.42), and Γ t the (non-truncated) generalized one-particle density matrix. The generalized Hartree-Fock-Bogoliubov operator is (4.46)

Well-Posedness of the Fermionic Bogoliubov-de Gennes Equations
The well-posedness of the effective equation obtained from the Dirac-Frenkel principle is not automatic; however under reasonable assumptions on the interaction potential and the initial data it can be established by standard methods. Since to our knowledge there is no proof completely spelled out in the literature, in this section we give a detailed proof that the time-dependent fermionic Bogoliubov-de Gennes equations are well-posed. We consider only the case h = −∆ (in particular no external potential V ext is included), and we are interested in interaction potentials V including the Coulomb potential.
In this section we also consider mixed states as initial data for the Bogoliubov-de Gennes equation, i. e. generalized one-particle density matrices satisfying 0 ≤ Γ 0 ≤ 1 (whereas before in the derivation we considered only pure states, Γ 2 0 = Γ 0 ). Well-posedness for similar equations has been discussed before, e. g. in [36,41] for a relativistic system (which generally exhibits finite-time blow-up, and global well-posedness only for small initial data) and in [3] for the bosonic Hartree-Fock-Bogoliubov equations. They are generalizations of earlier work on the Hartree-Fock equations without pairing (α = 0) [14,15,12,38]; see also [2]. All these works are applications of the abstract formalism developed by Segal [55].

Duhamel Formula and Integral Form
The standard approach to show local well-posedness is through an application of the Banach fixed-point theorem to the integral equation obtained from the Duhamel formula (in the spirit of the Picard-Lindelöff Theorem). This is the strategy we also follow here.
There is a small complication as it does not seem possible 5 to apply the Banach fixedpoint theorem in energy space (denoted Y below) when V has a Coulomb singularity. Instead we introduce an additional Banach space Z tailored to the Banach fixed-point theorem. We find local solutions in Z and we show afterward that these give rise to global energy space solutions. 5 Let V (x) be the Coulomb potential 1 |x| and M the Fourier multiplier √ 1 − ∆. In energy space (5.50) the norm is given by γ Y 1 = M γM S 1 and α Y 2 = α(·, ·) H 1 (R 3 ×R 3 ) . Consider the mild equation for γt: while we can control the term with [V(γs), γs] with sup 0≤s≤t γs 2 Y 1 , a homogeneity argument shows that we cannot control those with −ΠV (αs)αs and αsΠV (αs) with sup 0≤s≤t αs 2 Y 2 . Indeed we have tr ∂j This expression has contributions scaling like 1/[Length] 3 (the Coulomb potential, and one derivative contributing from each ∂j ) and is of degree 2 in αs. However, a bound in terms of αs 2 Y 2 can accommodate at most two derivatives, i. e. has contributions scaling at most like 1/[Length] 2 . The same holds for the term corresponding to αsΠV (αs). This provides us with reasonable doubts about the validity of the fixed-point scheme in energy space. To be controllable the two terms −ΠV (αs)αs and αsΠV (αs) cannot be separated, and the smoothing effect from the conjugation by e i∆(t−s) does not seem to help to estimate trace-class norms.
Recall the Bogoliubov-de Gennes equations from (1.11) (with h = −∆): Using Duhamel's formula, we obtain the integral form 6 of the Bogoliubov-de Gennes equations: where in the equation for α t , we are using the identification of the Hilbert-Schmidt operator α t with a two-particle wave function in L 2 (R 3 × R 3 ), denoting action of an operator on the i-th variable (i ∈ {1, 2}) by (·) i , and reading h := −∆ 1 −∆ 2 +V (x 1 −x 2 ) as a two-body Schrödinger operator.

Choice of Banach Spaces
For short we write S 1 = S 1 (L 2 (R 3 )) for the space of trace-class operators on L 2 (R 3 ), S 2 = S 2 (L 2 (R 3 )) for that of Hilbert-Schmidt operators, and B = B(L 2 (R 3 )) for that of linear bounded operators (equipped with the operator norm · B ). For another Banach space X , the space of linear bounded operators on X will be written B(X ) with norm · B(X ) . Let us introduce the Fourier multiplier M := (1 − ∆) 1/2 . We define the Banach space , usually called the energy space, by Here α T denotes the operator with integral kernel α T (x, y) = α(y, x). Here and below α(·, ·) H 1 refers to the norm of α in H 1 (R 3 × R 3 ). We define the Banach space Z := Z 1 × Z 2 for the purpose of constructing local solutions by the Banach fixed-point theorem by For any operator A we have By (5.51), since the Hilbert-Schmidt norm is smaller than the trace norm, we also have γ(·, ·) H 1 ≤ γ Z 1 . To shorten notations, we sometimes identify an operator with its integral kernel and write α H 1 , hα and e ith α.
Then Y and Z are invariant under the group action where e ith acts on the integral kernel α(·, ·) ∈ L 2 (R 3 × R 3 ).
Proof. The action on γ is given by conjugation with the unitary Fourier multipliers e is∆ , conserving the self-adjointness and the S 1 -,

Results on Well-Posedness
We now study the existence of solutions to the time-dependent Bogoliubov-de Gennes equations. As usual we expect them to conserve the number of particles tr(γ) and the energy of the system We need the following assumptions on the potential V : Observe that in fact the condition where C 3 > 0 only depends on the dimension. This notation of C V will be used throughout this paper.

Lemma 5.2 (Local Well-Posedness in Z)
. Assume that V satisfies (5.53). Consider a pair of initial data (γ 0 , α 0 ) ∈ Z. Then there exists a unique Z-continuous solution to the mild equations (5.48).
If [0, T ) is the maximal interval of existence of the solution, we have the usual blow-up alternative: either T = +∞ or lim t→T − ω t Z = ∞. Lemma 5.3 (Regularity of the Solution). If the pair of initial data (γ 0 , α 0 ) ∈ Z satisfies [−∆, γ 0 ] ∈ S 1 and α 0 (·, ·) ∈ H 2 , then the mild solution of the previous lemma is a strong The reader might wonder why we do not simply refer to [55,Lemma 3.1]. Transposed here, it states that if [−∆, γ 0 ] ∈ Z 1 and hα 0 ∈ Z 2 , that is α 0 ∈ dom S 2 (|h| 3/2 ), then the mild solution is a strong solution with time derivative in Z. When V is the Coulomb potential, even if a kernel α(·, ·) is of Schwartz class, hα is generally not in H 1 (R 3 × R 3 ), unless the diagonal α(x, x) identically vanishes. Here, thanks to the transpose symmetry α T 0 = −α 0 , the result still has an important value, but Lemma 5.3 is in some sense optimal: it states regularity when the minimal requirements are satisfied.
Remark. In the previous Lemma we have seen that the expected number of particles tr γ t is always a conserved quantity. Now let us introduce the quantity It coincides with 2 tr γ t if and only if Γ t describes a pure quasifree state, and it is also a conserved quantity because Γ t is unitarily equivalent to Γ 0 and because there holds We can interpret (5.54) as a measure of the deviation from Γ t being pure quasifree.
In the bosonic case, the same role is played by where P t = −Γ t S, P vac = −Γ vac S,γ t is the truncated one-particle density matrix, and Γ t is the truncated generalized one-particle density matrix. To see that this is conserved, simply notice that tr P 2 t − P t = tr Γ t SΓ t +Γ t S, and by [3, Lemma 3.10]Γ 0 = U tΓt U * t for some symplectomorphism U t (meaning that U * t SU t = S = U t SU * t ). As a word of caution: for bosonic systems, trγ t is not conserved by itself, only the particle number tr γ t = trγ t + ϕ t 2 is conserved.
The following lemma is used, together with conservation of the energy and the particle number, to globally control the Y-norm of a solution by the Y-norm of the initial data and thus to ensure that it does not blow-up.
Consider ω ∈ Y satisfying 0 ≤ γ 2 − αα ≤ γ (or equivalently 0 ≤ Γ ≤ 1 for the associated Γ). Then the following (crude) estimates hold: • for any δ > 0 there exists a C δ > 0, depending only on C V and δ, such that The proofs of these lemmas are postponed to the following subsections. We now state the main result of this section.

Estimates on the Non-Linearities
We state here results needed to control the operators V(γ) = V * ρ γ − X V (γ) and Π V (α) in the nonlinearities.
. For γ ∈ S 1 and α ∈ S 2 we have and for the multiplication operator V * ρ γ we have Proof of Lemma 5.7. 1. The estimate on Π V (α) S 2 is a simple application of the operator inequality we use the following trick explained in [13,Section 6]. We decompose γ = a + − a − , where a ± := M −1/2 (M 1/2 γM 1/2 ) ± M −1/2 , and (M 1/2 γM 1/2 ) ± ≥ 0 are, respectively, the positive and negative part of M 1/2 γM 1/2 in its spectral decomposition. By monotonicity of the square root, and writing down the spectral decomposition of a ± , we obtain For γ ∈ Z 1 self-adjoint, by splitting γ = γ + − γ − (positive and negative part), we find Now with (λ i ) i denoting the eigenvalues of M 1/2 γM 1/2 and (ϕ i ) i a corresponding orthonormal basis of eigenvectors, we thus get Using cyclicity of the trace in the last estimate, and afterward is an application of the Cauchy-Schwarz inequality and Fubini-Tonelli theorem: for ψ ∈ L 2 (R 3 ), let φ := M −1 ψ. Then 4. Let us now prove the estimate on X V (γ)M −1 S 2 ( Π V (α)M −1 S 2 is estimated in the same way). We use the idea of [3, Lemma E.1]. By cyclicity of the trace, for A ∈ S 2 , we have A S 2 = A * S 2 , and it suffices to show the boundedness of As V (x) = V (−x), the trace is equal to where M −2 (y − z) denotes the Yukawa potential at point y − z. For almost all x, we consider the L 2 -function g x (y) := γ(x, y)V (x − y). By operator monotonicity of the inverse, we have M −2 We obtain the upper bound: We introduce the "polarization" of K 1 and K 2 as a bilinear form (it is not necessary to give it a symmetrized form): Lemma 5.8. Let V satisfy (5.53). Then for the nonlinearities K 1 and K 2 seen as bilinear maps (as defined below in (5.55) and (5.56)) we have the estimates where the constant C depends only on C V .
Proof. The estimates follow directly from Lemma 5.7.

Local Well-Posedness in Z (Proof of Lemma 5.2)
The norm on C(I, Z) is given by (γ, α) := sup t∈I (γ t , α t ) Z . Denoting the initial values by γ 0 = γ and α 0 = α, we define the Picard operator ♣ : C(I, Z) → C(I, Z) on an interval I containing t = 0 by setting (with K 1 and K 2 as before) To establish local existence we show that the nonlinearities (when simply taking the norm inside the integral and neglecting the unitaries-for the second equation this is possible because h is infinitesimally operator bounded w. r. t. the Laplacian) are locally Lipschitz; then ♣ on a sufficiently short time interval I is a contraction, and a solution is found by the Banach fixed-point theorem. We refer to [55, Theorem 1] for details. To prove local Lipschitz continuity of the quadratic terms K 1 and K 2 , it suffices to show continuity of the corresponding bilinear maps K 1 (·, ·) and K 2 (·, ·) introduced in (5.55)-(5.56). Indeed, we can write down the difference of the quadratic terms in terms of the polarization according to the following formula: Lemma 5.9 (Continuity of the Polarized Non-Linearities). The bilinear forms K 1 (·, ·) and K 2 (·, ·) (5.55)-(5.56) are continuous from Z 2 to Z 1 and from Z 2 to Z 2 respectively. Their norms are bounded by a constant depending only on C V . Proof. We have to estimate K i (ω 1 , ω 2 ) Z i (i = 1, 2) in terms of ω 1 Z and ω 2 Z .
First, let us establish a formula allowing us to pull derivatives through the interaction operators. The crucial fact is that the multiplication operator V (x − y) commutes with [∇, ·] = ∇ x − ∇ y on S 2 , that is: [∇, X V (γ)] = X V ([∇, γ]) and the analogous for Π V (α). Similarly we have ∇(V * ρ γ ) = V * ρ [∇,γ] (with V * ρ γ read as a function; to prove this identity we use the spectral decomposition of γ) and [∇, V * ρ γ ] = V * ρ [∇,γ] (with V * ρ γ read as a multiplication operator on L 2 (R 3 )). We obtain the pull-through identity which enables us to transfer an M from the left side of V(γ) to its right side. The same calculation holds for Π V (α).
We can now prove continuity of the nonlinearity K 1 . Recall that by definition From (5.55) we get Now we employ the pull-through formula (5.58) and afterward Lemma 5.7, as well as the fact that M −1 and ∂ j M −1 are bounded operators: Similarly, using the pull-through formula for Π V (α), we obtain Combining everything we get M K 1 (ω 1 , ω 2 ) S 1 ≤ C ω 1 Z ω 2 Z . We can estimate the term K 1 (ω 1 , ω 2 )M S 1 in the same way. We conclude that To prove continuity of the nonlinearity K 2 , we use (5.51) as an upper bound for K 2 2 H 1 , and then proceed by the same method as for K 1 . We obtain

Regularity of the Solution (Proof of Lemma 5.3)
We now assume that the initial data satisfy also the additional regularity conditions [γ 0 , −∆] ∈ S 1 and α 0 (·, ·) ∈ H 2 (or equivalently α 0 (·, ·) ∈ dom S 2 (h)). It suffices to adapt Segal's result [55,Lemma 3.1] to ensure that the solution is a strong solution: instead of the Z-norm, we apply the same argument to the S 1 × S 2 -norm.
Identifying γ t and α t with their integral kernel, the derivativesγ t ,α t are well-defined as bounded linear operators from H 2 (R 3 ) to its dual H −2 (R 3 ), and they are equal to −i times the r. h. s. of (5.47). This establishes the equations, and it just remains to prove that ω t is Fréchet-differentiable in S 1 × S 2 .
We construct the putative derivatives by the Banach fixed-point theorem. We start by formally differentiating the Bogoliubov-de Gennes equations to see that the derivatives should satisfy the equations where ∂ ∂ω K 1 (ω t ) and ∂ ∂ω K 2 (ω t ) denotes the Fréchet derivative of the nonlinearities. As initial data for the fixed-point problem of the derivatives we have (as given by the Bogoliubov-de Gennes equations) We now write the equations forγ t andα t in mild form: (5.59) As K 1 and K 2 are quadratic functions of ω, we have for all δω ∈ Z ∂ ∂ω K j (ω t )δω = K j (ω t , δω) + K j (δω, ω t ), (j = 1, 2) where the K j (·, ·)'s are the (un-symmetrized) polarizations as defined in (5.55) and (5.56). By Lemma 5.8, their extensions ∂ ∂ω K j (ω t ) : S 1 × S 2 → S j to S 1 × S 2 ⊃ Z are continuous with norm controlled by C V ω t Z . So we can apply the Banach fixed-point theorem to (5.59) in the Banach space S 1 × S 2 . We obtain a unique local solution v t = (g t , a t ) t∈[0,T ) .
It is easy to see through a simple Grönwall argument that v t S 1 ×S 2 growths at most like exp Ct sup s∈[0,t] ω s Z , hence the maximal interval of existence of v t is the same as that of ω t .

Existence of Unitary Propagator and Conservation of Spectrum
In this section we prove that the solution Γ t at time t of the Bogoliubov-de Gennes equation is related to Γ 0 by conjugation with a unitary propagator. This implies that the spectrum of Γ t is time-independent.
For regular initial data We start with the case where [−∆, γ t ] ∈ S 1 and α(·, ·) ∈ H 2 (R 6 ). We split F Γt into unbounded time-independent and bounded time-dependent part as Let us now verify the C 1 -condition. Since we have assumed regular initial data, the mild solution Γ t is continuously differentiable. Since ϕ ∈ H 2 (R 3 ), we can insert 1 = M −1 M and obtain d dt Using the estimates of Lemma 5.7, we find that t → B(∂ t Γ t )M −1 is continuous (w. r. t. the operator norm).
Consider the evolution S t := U (t, 0)Γ 0 U (0, t). If we prove that it satisfies the integral form of the equation , then it follows by Grönwall's uniqueness argument applied in the space To verify that S t satisfies (5.61), it suffices (by density) to show that the expectation value φ, S t ψ tested with functions ψ, φ ∈ H 2 (R 3 ) 2 satisfies the integral equation (i. e. a weak formulation of the integral equation). Since H 2 (R 3 ) 2 = D(A+B t ), the expectation value φ, S t ψ is differentiable with derivative −i φ, [F Γt , S t ]ψ , and by the standard Duhamel trick we find that it satisfies the integral version of (5.61), Extension to non-regular initial data It remains to extend to the case when the initial data have less regularity. Given arbitrary (γ, α) ∈ Z, we regularize them by setting γ n := P −∆≤n γP −∆≤n and α n := P −∆≤n αP −∆≤n .
Consider the Z-solutions Γ (n) t resp. Γ t of the Bogoliubov-de Gennes equation with initial data (γ n , α n ) resp. (γ, α). As (γ n , α n ) Z ≤ (γ, α) Z , they are all defined at least on a common interval [0, T 1 ] (the interval used for the Banach fixed-point scheme depends only on the norm of the initial data and on C V ). By a simple Grönwall argument they converge to Γ t in C([0, T 1 ], Z).
By the argument we gave above, for the solution Γ Both equations are solvable on [0, T 1 ] by applying the Banach fixed-point theorem (hence the subscript "BFP") in the Banach space of bounded operators, and the obtained solutions are as usual unique. However, a priori we do not know that these solutions are unitaries. But the solution U for every fixed t, in operator norm; this will imply the unitarity of U BFP (t).
The convergence U Taking the operator norm and using Lemma 5.8, we obtain The first summand converges to zero as n → ∞. By Grönwall's method, we now have U (n) (t, 0) = U Finally, this implies that, as n → ∞, U (n) (t, s) = U

Conservation of the Particle Number tr(γ)
The conservation is easy to establish for strong solutions by differentiating the particle number tr γ t . In fact, consider regular initial data (γ 0 , α 0 ), i. e. [−∆, γ 0 ] ∈ S 1 and α 0 (·, ·) ∈ H 2 . Then by Lemma 5.3, we can freely differentiate, and find The first trace vanishes since it can be written as the derivative of a function which is constant due to cyclicity of the trace, i. e. .
The second trace vanishes by cyclicity (note that V(γ t ) is bounded). The third trace vanishes since we can write it out as an integral and use V (x) = V (−x). We now turn to arbitrary initial data in Z. Since we have existence of solutions due to a Banach fixed-point argument in Z, the solutions are continuous in Z-norm, w. r. t. initial data in Z. The number of particles tr γ t is obviously Z-continuous, and so by approximating Z-initial data by regular initial data, tr γ t is constant again.
Regularization Since the kinetic part of the energy functional is not Z-continuous, we cannot use the same strategy as for particle number conservation. Instead, we introduce a regularization for which the conservation of energy holds. As before P Λ denotes 1 (−∆<Λ) , and to shorten notations, we also denote P Λ ⊗ 1 C 2 (acting on L 2 (R 3 ) 2 ) by P Λ . We regularize both the equation and the functions: for any Λ > 0, we consider the solution (Γ (Λ) Above, P Λ F St P Λ denotes the bounded operator: The are obtained from the original ones K 1 , K 2 by replacing V(γ) and Π V (α) by P Λ V(γ)P Λ and P Λ Π V (α)P Λ respectively.
We can apply the Banach fixed-point theorem to the regularized equations with the · Z -norm: the estimates are the same, and the interval of existence [0, T ) only depends on the initial data Γ 0 , it does not depend on the cutoff Λ > 0. For any Λ > 0, we thus obtain a solution (Γ for the corresponding objects of the regularized solution.
As the operator (5.64) Since the regularized equation has the same structure as the original one, conservation of the spectrum still holds (by the argument we gave before). In particular we have: 0 ≤ Γ (Λ) t ≤ 1. We show consecutively the following four points.
Let us show how we can then establish the conservation of the energy. The conservation of the spectrum gives 0 ≤ g t ≤ 1 for g t = γ (Λ) t and g t = γ t . Together with point (iv) it ensures that (ω t ) 0≤t<T is Y-valued. Indeed for any Λ 0 > 0, we have by Fatou's lemma: Taking the limit Λ 0 → ∞ yields tr(−∆γ t ) < ∞ by monotone convergence, and we obtain for all 0 ≤ t ≤ T 1 < T the inequality We have to show equality in (5.66), i. e. that there is no loss of mass as Λ → ∞ (more precisely no loss of H 1 -mass of the eigenfunctions of γ (Λ) t ). Equality in (5.66) is ensured by the time-reversal symmetry of the equation. Indeed for any 0 < T 1 < T , the path t ∈ [0, T 1 ] → Γ T 1 −t satisfies the same equation as (Γ t ) 0≤t<T , hence the same arguments as above give the reverse inequality of (5.66). We emphasize that the argument uses the obvious equality E(Γ) = E(Γ) and the local uniqueness (due to the Banach fixed-point argument) of the solution.
(i) Consistency of the regularization By Grönwall's method we show that where A was defined in (5.60) and B .
We thus get the inequality Similarly Γ (Λ) t (1 − P Λ ) B satisfies the same inequality and we obtain For any 0 < T 1 < T , we have Thus f t = 0 on [0, T 1 ] for all 0 < T 1 < T , that is f t identically vanishes on [0, T ).
(ii) Conservation laws for the regularized problem Point (i) ensures us that ω (Λ) t is Y-differentiable on [0, T ), and that We then observe that the energy is invariant under complex conjugation on Y, and can write the energy functional as where E HF (γ) = tr(−∆γ) + 1 2 tr (V(γ)γ). Notice furthermore that due to the assumption V (x) = V (−x) we have the identity tr Π V (α)β = tr αΠ V (β).
Taking explicitly the time derivative of E(Γ (Λ) t ) and using the last identity, it is a straightforward calculation to see that Conservation of the number of particles is proven similarly (and with less calculations).
(iii) Convergence of the regularization Let 0 < T 1 < T . We consider the mild form of the equation on α which gives the integral part in the integral inequality (5.67). All we have to show is that the lines (5.69) through (5.71) also converge to 0 as Λ → 0. By a similar approach we can deal with the terms of the decomposition of γ (Λ) t − γ t , and both estimates will give (5.67). We only estimate α (Λ) t − α t S 2 and leave γ (Λ) t − γ t S 1 to the reader.
We first describe two technical results, which will then be useful in dealing with the remaining lines. We emphasize that e −ish and e −ish Λ , s ∈ R are unitary operators which leave the Hilbert-Schmidt norm invariant. For γ (Λ) t − γ t , the conjugation by the unitary operators e is∆ leaves the trace-norm invariant.
• The first technical issue is to deal with the convergence of e −ish Λ , for s ∈ R. The key observation is that h Λ converges to h in the strong-resolvent sense. Indeed the resolvent identity gives For α ∈ S 2 , the integral kernel of K := (h + i) −1 α is in H 2 (R 3 × R 3 ). By compactness of K, we have Similarly, as the operator Π V (K) is Hilbert-Schmidt we have Then [52,Theorem VIII.20] ensures that for any bounded Borelian function f : R → C, the operator f (h Λ ) converges to f (h) in the strong operator topology. In particular for all s ∈ R, e −ish Λ converges to e −ish in the strong operator topology.
• For the second technical issue, we introduce a second level of cutoff Λ ′ > 0. For 0 ≤ s ≤ T 1 , the operator K 2 (ω s ) is compact (and its integral kernel is in H 1 ), hence we have point-wise in s: On [0, T 1 ], the norm K 2 (ω s ) S 2 is uniformly bounded by sup t∈[0,T 1 ] K 2 (ω t ) < ∞.
Hence by dominated convergence, we obtain Consider now the second line (5.69). Taking Λ ′ = Λ, we get that it converges to 0 as Λ → ∞.
Consider now the fourth line (5.71). The first term converges to 0 since we have lim Λ→0 P Λ α 0 P Λ − α 0 S 2 = 0. For its second term, splitting α 0 into two: The second summand vanishes as Λ ′ → ∞. At fixed Λ ′ , we then have This follows from the convergence in the strong operator topology of e −ith Λ and the fact that P Λ ′ α 0 P Λ ′ ∈ dom S 2 (h) = H 2 which gives by functional calculus the crude estimate where t 1 , t 2 ∈ R and where h denotes h or h Λ . Hence by an ε/2-argument, the fourth line converges to 0 as Λ → ∞.
As in Lemma 5.8, by using Lemma 5.7 we obtain the following estimate: where Γ i denotes the generalized density matrix corresponding to ω i . Point (ii) and Point (iii), namely the energy conservation of Γ
From the estimates above (and adjusting the choice of δ) we conclude that Finally, by using the symmetry α T = −α and going to Fourier space, we get