Stefano Marcantoni

Marcello Porta

Julien Sabin

Dynamics of mean-field Fermi systems with nonzero pairing

We study the dynamics of many-body Fermi systems, for a class of initial data which are close to quasi-free states exhibiting a nonvanishing pairing matrix. We focus on the mean-field scaling, which for fermionic systems is naturally coupled with a semiclassical scaling. Under the assumption that the initial datum enjoys a suitable semiclassical structure, we give a rigorous derivation of the time-dependent Hartree-Fock-Bogoliubov equation, a nonlinear effective evolution equation for the one-particle density matrix of the system, as the number of particles goes to infinity. Our result holds for all macroscopic times, and provides bounds for the rate of convergence.

Introduction

Let us consider a system of N ≫ 1 quantum particles in R 3 , with M internal degrees of freedom, or spins. At zero temperature, the state of the system is denoted by the wave function ψ N ∈ L 2 ((R 3 × {1, . . . , M }) N ), under the normalization condition

∥ψ N ∥ 2 2 := σ 1 ,...,σ N dx 1 . . . dx N |ψ N (x 1 , . . . , x N )| 2 = 1 , (1.1)
where x = (x, σ) with x ∈ R 3 the position label and σ ∈ {1, . . . , M } the spin label. We shall restrict the attention to fermionic particles, which correspond to antisymmetric wave functions: ψ N (x 1 , . . . x N ) = sgn(π)ψ N (x π(1) , . . . x π(N ) ) ,

where π is a permutation of the particle labels {1, . . . , N }. We shall denote by L 2 a ((R 3 × {1, . . . , M }) N ) the set of all antisymmetric wave functions in L 2 ((R 3 × {1, . . . , M }) N ). A simple example of fermionic wave functions is provided by Slater determinants, defined as antisymmetric tensor products of orthonormal functions in L 2 (R 3 × {1, . . . , M }):

ψ Slater (x 1 , . . . x N ) = 1 √ N ! π∈S N sgn(π)f 1 (x π(1) ) • • • f N (x π(N ) ) , (1.3)
with ⟨f i , f j ⟩ = δ ij , and where S N is the set of permutations of {1, 2, . . . , N }. These wave functions encode the least amount of correlations between the particles, besides the requirement of antisymmetry.

Let us now introduce the many-body Hamiltonian, which generates the dynamics of the system. A typical class of non-relativistic Hamiltonians is given by:

H trap N = N j=1 (-∆ j + V trap (x j )) + λ N i<j V (x i -x j ) , (1.4)
where V trap is a confining potential, and V is the two-body interaction, with coupling constant λ. The ground state energy, which physically corresponds to the energy of the system at zero temperature, is given by:

E N = inf ψ N ∈L 2 a ((R 3 ×{1,...,M }) N ) ⟨ψ N , H N ψ N ⟩ ∥ψ N ∥ 2 2
.

(1.5)

Let ψ N be a minimizer of H trap N . Clearly, ψ N is a stationary state for the quantum dynamics generated by H trap N . In order to observe a non-trivial evolution, we turn off the external potential, and evolve the state according to the time-dependent Schrödinger equation:

i∂ τ ψ N,τ = - N j=1 ∆ j + λ N i<j V (x i -x j ) ψ N,τ , ψ N,0 = ψ N . (1.6)
In physical applications, it is extremely hard to obtain quantitative understanding on the ground state of the system or on the solution of the Schrödinger equation, for realistic values of the number of particles N . One is therefore interested in approximate descriptions, which depend on much fewer degrees of freedom, and which allow to capture the collective behavior of the system on a macroscopic scale. The price to pay for the reduction in dimensionality is that the laws of motion for the effective theories are typically nonlinear. The paradigmatic scaling regime in which such effective nonlinear description is expected to arise is the meanfield scaling. In this setting, the potential V trap in (1.4) confines the system in a region whose volume is of order 1 in N , and the coupling constant λ is chosen to guarantee the balance between the kinetic and the interaction energy. Due to the growth as N 5/3 of the kinetic energy of the Fermi gas (which could be proved using the Lieb-Thirring inequality, see e.g. [START_REF] Lieb | The Stability of Matter in Quantum Mechanics[END_REF]), the natural choice for the coupling constant is λ = O(N -1/3 ). In this scaling regime, the equilibrium and the dynamical properties of the system are expected to be well captured by Hartree-Fock theory, defined as the restriction of the space of fermionic wave functions to the subsets of Slater determinants (1.3). The advantage in considering this class of states is that they behave as Gaussian states, in the following sense. Given ψ N ∈ L 2 a ((R 3 × {1, . . . , M }) N ) we define its k-particle reduced density matrix as the linear operator on L 2 ((R 3 × {1, . . . , M }) k ) with integral kernel:

γ (k) N (x 1 , . . . , x k ; y 1 , . . . , y k ) = N k dz k+1 • • • dz N ψ N (x 1 , .
. . , x k , z k+1 , . . . , z N )ψ N (y 1 , . . . , y k , z k+1 , . . . , z N ) ,

with the notation:

dz f (z) = M σ=1 R 3
dz f (z, σ) .

(1.8)

The reduced k-particle density matrix is important because it allows to compute expectation values of k-particle observables; it plays the same role as the k-particle marginal in classical statistical mechanics. For Slater determinants, it turns out that all k-particle reduced density matrices are determined by the one-particle density matrix. Let ω N be the one-particle density matrix of a Slater determinant (1.3). It is not difficult to see that:

ω N = N i=1 |f i ⟩⟨f i | ; (1.9)
in particular, observe that ω N is a rank-N orthogonal projector. The k-particle density matrix of a Slater determinant can be computed starting from the one-particle density matrix (1.9) via the application of the fermionic Wick's rule: thus, all the information about the state described by a Slater determinant is encoded in the reduced one-particle density matrix ω N . For instance, the energy of a Slater determinant is a functional of the one-particle density matrix only:

⟨ψ Slater , H trap N ψ Slater ⟩ = tr(-∆ + V ext )ω N + 1 2N 1/3 dxdy V (x -y) ρ(x)ρ(y) -|ω N (x; y)| 2 =: E HF N (ω N ) , (1.11)
with ρ(x) = ω N (x; x) the density of particles at x; E HF N (ω N ) is called the Hartree-Fock energy functional.

Hartree-Fock theory can also be used to investigate the dynamical properties of the system. Observe preliminarily that, in this confined setting, the relevant time scale of the system is τ = O(N -1/3 ), due to the fact that the typical kinetic energy per particle grows as N 2/3 . From now on, we will write τ = N -1/3 t, and we will be interested in studying the dynamics of the system for times t of order 1 in N . Suppose that at time t = 0 the external potential is turned off. The Schrödinger equation in the rescaled time variable reads:

iε∂ t ψ N,t = - N j=1 ε 2 ∆ j + 1 N N i<j V (x i -x j ) ψ N,t , ψ N,0 = ψ N , (1.12)
where ε = N -1/3 . Eq. (1.12) shows that the fermionic mean-field can be understood as the coupling of the usual mean-field with a semiclassical scaling. If one assumes that the state of the system is a Slater determinant for all times, then it is not difficult to see that the oneparticle density matrix of the system evolves according to the time-dependent Hartree-Fock equation:

iε∂ t ω N,t = [h HF (t), ω N,t ] , ω N,0 = ω N , (1.13)
where h HF (t) = -ε 2 ∆+N -1 V * ρ t -X t , with ρ t (x) = ω N,t (x; x) and X t is the linear operator with integral kernel given by X t (x; y) = N -1 V (x -y)ω N,t (x; y).

Of course, the restriction to the set of Slater determinants is a nontrivial assumption, and it is a natural mathematical question to understand when such approximation is justified. Concerning the equilibrium problem, the validity of the Hartree-Fock approximation for the ground state of many-body fermionic systems has been proved in [START_REF] Bach | Error bound for the Hartree-Fock energy of atoms and molecules[END_REF] for large atoms and molecules, and it has been extended to high density Coulomb systems (jellium) in [START_REF] Graf | A correlation estimate with applications to quantum systems with Coulomb interactions[END_REF]. More recently, there has been substantial progress in the development of nonperturbative methods to go beyond the Hartree-Fock approximation, and in particular to determine the large N asymptotics of the correlation energy, defined as the difference between the many-body ground state energy minus the Hartree-Fock ground state energy; see [START_REF] Hainzl | On the correlation energy of interacting fermionic systems in the mean-field regime[END_REF][START_REF] Benedikter | Optimal upper bound for the correlation energy of a Fermi gas in the mean-field regime[END_REF][START_REF] Benedikter | Correlation Energy of a Weakly Interacting Fermi Gas[END_REF][START_REF] Christiansen | The Random Phase Approximation for Interacting Fermi Gases in the Mean-Field Regime[END_REF][START_REF] Benedikter | Correlation Energy of a Weakly Interacting Fermi Gas with Large Interaction Potential[END_REF][START_REF] Christiansen | The Gell-Mann-Brueckner Formula for the Correlation Energy of the Electron Gas: A Rigorous Upper Bound in the Mean-Field Regime[END_REF].

Concerning quantum dynamics, the first derivation of the time-dependent Hartree-Fock equation (1.13) in the mean-field/semiclassical scaling has been obtained in [START_REF] Elgart | Nonlinear Hartree equation as the meanfield limit of weakly coupled fermions[END_REF], for analytic potentials and for short times. More recently, in [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF] (and later in [START_REF] Petrat | A new method and a new scaling for deriving fermionic mean-field dynamics[END_REF] with a different method), this result has been generalized to a much larger class of potentials and to all times. The result of [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF] holds in the sense of convergence of density matrices. Convergence in the L 2 -sense for a special class of initial data has been obtained in [START_REF] Benedikter | Bosonization of Fermionic Many-Body Dynamics[END_REF], using the rigorous bosonization methods of [START_REF] Benedikter | Optimal upper bound for the correlation energy of a Fermi gas in the mean-field regime[END_REF][START_REF] Benedikter | Correlation Energy of a Weakly Interacting Fermi Gas[END_REF][START_REF] Benedikter | Correlation Energy of a Weakly Interacting Fermi Gas with Large Interaction Potential[END_REF]. The result of [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF] has been extended to pseudo-relativistic fermions in [START_REF] Benedikter | Mean-field dynamics of fermions with relativistic dispersion[END_REF] and to mixed states in [START_REF] Benedikter | Mean-field Evolution of Fermionic Mixed States[END_REF]. In the mean-field regime without semiclassical scaling, the derivation of the time-dependent Hartree-Fock equation has been obtained in [START_REF] Bardos | Mean-field dynamics of fermions and the time-dependent Hartree-Fock equation[END_REF], for bounded potentials, and in [START_REF] Fröhlich | A microscopic derivation of the time-dependent Hartree-Fock equation with Coulomb two-body interaction[END_REF], for Coulomb potentials (see also [START_REF] Bach | Kinetic energy estimates for the accuracy of the time-dependent Hartree-Fock approximation with Coulomb interaction[END_REF]). See [START_REF] Benedikter | Effective evolution equations from quantum dynamics[END_REF] for a review. Concerning unbounded potentials in the mean-field/semiclassical scaling, the timedependent Hartree-Fock equation for particles interacting through a Coulomb potential has been derived in [START_REF] Porta | Mean Field Evolution of Fermions with Coulomb Interaction[END_REF], under the assumption that a suitable semiclassical structure of the initial datum propagates along the flow of the Hartree-Fock equation. The propagation of the semiclassical structure for a class singular potentials has been proven in [START_REF] Chong | From many-body quantum dynamics to the Hartree-Fock and Vlasov equations with singular potentials[END_REF] for mixed states; this work also allows to consider the Coulomb potential, for short times. Finally, we mention the recent work [START_REF] Fresta | Effective Dynamics of Extended Fermi Gases in the High-Density Regime[END_REF], where the derivation of [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF] has been generalized to the case of extended Fermi gases at high density, for bounded potentials.

Slater determiants are a special example of quasi-free states. In general, a quasi-free state is a state for which all correlation functions can be computed according to the Wick rule. If one works in the grand-canonical ensemble and relaxes the constraint that the number of particles is fixed, there exists other quantum states at zero temperature besides Slater determinants which are quasi-free. In order to introduce them, we will use the notations of second quantization and of the Fock space, introduced in detail in Section 2.

A pure state is described by a vector in the Fock space, ψ ∈ F = C ⊕ n≥1 L 2 a ((R 3 × {1, . . . , M }) n ), with ∥ψ∥ F = 1. In general, the vector ψ corresponds to a superposition of states with different particle numbers. The reduced one-particle density matrix of ψ can be expressed in Fock space as: γ

(1)

ψ (x; y) = ⟨ψ, a * y a x ψ⟩ , (1.14)
with a x and a * y the operator-valued distributions associated with the usual fermionic creation and annihilation operators. Since the state ψ does not have a definite particle number, we can also introduce a nontrivial pairing matrix:

π (1)
ψ (x; y) = ⟨ψ, a y a x ψ⟩ .

(1.15)

We shall denote by N the average number of particles of the state, defined as N = tr γ

ψ . In this more general setting, quasi-free states are completely determined by γ (1)

ψ ≡ ω N and π (1)
ψ ≡ α N , via a generalization of the Wick's rule (1.10); see Section 2.3. Furthermore, as discussed in Section 2.3, it turns out that the reduced one-particle density matrix and the pairing matrix of a pure quasi-free state satisfy the constraints:

0 ≤ ω N ≤ 1 , α T N = -α N , ω N (1 -ω N ) = α N α * N , ω N α N -α N ω N = 0 . (1.16)
Next, let us introduce the second quantization of the Hamiltonian H trap N ,

H trap N = 0 ⊕ n≥1 n j=1 (-∆ j + V ext (x j )) + 1 N 1/3 n i<j V (x i -x j ) .
(1.17)

The energy of a quasi-free state ψ is completely specified by ω N and α N . We have:

⟨ψ, H trap N ψ⟩ = tr (-∆ + V trap )ω N + 1 2N 1/3 dxdy V (x -y) ρ(x)ρ(y) -|ω N (x; y)| 2 + |α N (x; y)| 2 =: E HFB N (ω N , α N ) , (1.18)
where E HFB N (ω N , α N ) is called the Hartree-Fock-Bogoliubov (HFB) functional. This functional plays an important role in physics, for instance in the Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity and in superfluidity. It differs from the Hartree-Fock functional (1.11) because of the presence of the last term, taking into account the energy of Cooper pairs. In general, it is a difficult problem to determine whether the ground state of the HFB functional is attained in correspondence with a nonzero pairing matrix; a necessary condition for this to happen is that the potential is non-positive. See [START_REF] Bach | Generalized Hartree-Fock theory and the Hubbard model[END_REF][START_REF] Hainzl | The BCS functional of superconductivity and its mathematical properties[END_REF] for mathematical reviews of BCS theory and of the Hartree-Fock-Bogoliubov approximation, as well as [START_REF] Hainzl | The BCS Functional for General Pair Interaction[END_REF][START_REF] Bach | Bogolubov-Hartree-Fock mean field theory for neutron stars and other systems with attractive interactions[END_REF][START_REF] Lenzmann | Minimizers for the Hartree-Fock-Bogoliubov theory of neutron stars and white dwarfs[END_REF][START_REF] Lewin | A numerical perspective on Hartree-Fock-Bogoliubov theory[END_REF][START_REF] Bräunlich | Bogolubov-Hartree-Fock theory for strongly interacting fermions in the low density limit[END_REF][START_REF] Lauritsen | The BCS Energy Gap at Low Density[END_REF][START_REF] Henheik | The BCS Energy Gap at High Density[END_REF] for recent results.

Let us now consider the many-body evolution of an initial quasi-free state, after removing the external trap. The quantum dynamics of the system in Fock space, parametrized by the macroscopic time variable t, is defined as

iε∂ t ψ t = H N ψ t , ψ 0 = ψ , (1.19)
with H N given by:

H N = 0 ⊕ n≥1 - n j=1 ε 2 ∆ j + 1 2N n i<j V (x i -x j ) . (1.20)
Similarly to the case of Hartree-Fock theory, if one assumes that the state remains quasi-free for all times, one obtains the following coupled nonlinear evolution equations for the reduced one-particle density matrix and for the pairing matrix, which generalize (1.13) to the case in which the number of particles is not fixed:

iε∂ t ω N,t = h HF (t), ω N,t + Π t α * N,t -α N,t Π * t iε∂ t α N,t = h HF (t)α N,t + α N,t h HF (t) + Π t 1 -ω N,t -ω N,t Π t , (1.21)
where Π t a linear operator with integral kernel given by Π t (x, y) := 1 N V (x -y)α N,t (x, y). The equation (1.21) for the pair (ω N,t , α N,t ) is called the Hartree-Fock-Bogoliubov (HFB) equation. HFB theory is expected to be a good approximation of many-body quantum mechanics in the mean-field regime, for many-body systems with attractive interactions. However, in contrast to the Hartree-Fock case, until now there has been no rigorous proof of the validity of the HFB approximation from the many-body Schrödinger equation, neither at equilibrium nor for quantum dynamics. See [START_REF] Benedikter | The Dirac-Frenkel Principle for Reduced Density Matrices, and the Bogoliubov-de-Gennes Equations[END_REF] for the proof of well-posedness of the fermionic HFB equations, and for the derivation of the HFB equation (1.21) via the Dirac-Frenkel principle. See also [START_REF] Hainzl | On blowup for time-dependent generalized Hartree-Fock equations[END_REF] for the proof of blow-up at finite time of the relativistic fermionic HFB equations.

In this paper we prove the validity of the Hartree-Fock-Bogoliubov equation (1.21) for bounded and regular enough interaction potentials, in the mean-field regime. We shall consider the many-body dynamics of initial data ψ ∈ F that are approximated by a suitable class of quasi-free states with one-particle density matrix ω N and with pairing matrix α N , with N = tr ω N . We shall suppose that:

sup p∈R 3 1 1 + |p| ∥[ω N , e ip•x ]∥ HS ≤ CN 1/3 , ∥[ω N , ε∇]∥ HS ≤ CN 1/3 ∥[α N , ε∇]∥ HS ≤ CN 1/3 , ∥α N ∥ HS ≤ CN 1/3 . (1.22)
We shall suppose that the initial datum ψ is close, in a suitable sense to be specified below, to a quasi-free state parametrized by ω N , α N , satisfying the bounds (1.22). In particular:

∥γ (1) ψ -ω N ∥ HS ≪ N 1/2 , ∥π (1)
ψ -α N ∥ HS ≪ N 1/3 , (1.23)
that is, ∥γ

ψ -ω N ∥ HS /N 1/2 and ∥π

ψ -α N ∥ HS /N 1/3 converge to zero as N → ∞. The estimates (1.22) encode the semiclassical structure of the initial datum, and are expected to hold for quasi-free states at low energy. For α N = 0, similar estimates play a key role in the derivation of the time-dependent Hartree-Fock equation in the meanfield/semiclassical scaling [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF], and can be proved to hold for the ground state of confined, non-interacting systems [START_REF] Fournais | An optimal semiclassical bound on commutators of spectral projections with position and momentum operators[END_REF]. In Appendix A we prove the validity of the bounds (1.22) for translation-invariant quasi-free states on the torus of side 2π, at low enough energy. If translation-invariance is assumed, the lowest energy state of the HFB functional on the torus turns out to be the free Fermi gas, which in particular has α N = 0. This is ultimately a consequence of the existence of a spectral gap for the relative kinetic energy at the Fermi level, due to the finite volume of the torus. Thus, for confined translation-invariant states, the condition α N ̸ = 0 implies that the state is an excited state. We do not know whether the ground state of the HFB functional is attained in correspondence with a nonzero pairing matrix. The persistence of translation-invariance for the ground state of the BCS functional has been discussed in [START_REF] Deuchert | Persistence of translational symmetry in the BCS model with radial pair interaction[END_REF], in the thermodynamic limit. There it is shown that translationinvariance persists in two-dimensions, and also in three-dimensions in the case of vanishing angular momentum.

Let us denote by γ

N,t and π

N,t the many-body evolution of the reduced one-particle density matrix and of the pairing matrix associated with the solution of (1.19), and by ω N,t , α N,t the solution of Eqs. (1.21). We prove that:

∥γ (1) N,t -ω N,t ∥ HS ≪ C(t)N 1/2 , ∥π (1)
N,t -α N,t ∥ HS ≪ C(t)N 1/3 ∀t ∈ R, (1.24)
where the function C(t) grows as a double exponential in time, and it is independent of N . The bounds (1.24) prove convergence from many-body dynamics to the HFB dynamics at the level of the one-particle density matrices. More generally, our result also extends to higher order density matrices. The method of proof is based on an extension of the technique of [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF]. As in [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF], a key role is played by the fluctuation dynamics, a unitary evolution in Fock space that does not preserve the number of particles, and that allows to understand the distance between the many-body state and the reference quasi-free states in term of particle excitations over a suitable time-dependent Fock space vacuum. With respect to the Hartree-Fock case, however, in the present setting the control of the fluctuation dynamics is considerably more involved, due to the more complex structure of the effective evolution equation (1.21), and to the fact that ω N is not an orthogonal projection. Also, the propagation of the semiclassical structure (1.22) along the HFB flow is more involved than in the HF case. Let us point out that the proof could be simplified, at the cost of assuming a stronger notion of semiclassical structure (1.22), replacing in (1.22) the square of the Hilbert-Schmidt norm by the trace norm, similarly to [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF]. Furthermore, with this stronger notion of semiclassical structure we could also prove that the Hartree-Fock-Bogoliubov dynamics cannot be distinguished from the Hartree-Bogoliubov (HB) dynamics, given by Eqs. (1.21) after neglecting all terms involving the operators X t and Π t : the distance between the solutions of the HFB and HB evolutions can be proved to be smaller than the estimate quantifying the distance between many-body and HFB dynamics (1.24). Since we have not been able to verify the trace-norm semiclassical structure for interesting quasi-free states with non-zero pairing, we prefer to work under the weaker semiclassical structure given by (1.22).

In order to deal with the weaker form of semiclassical structure (1.22), we introduce a suitable auxiliary fluctuation dynamics, as already done in the case of fermionic mixed states (without pairing) [START_REF] Benedikter | Mean-field Evolution of Fermionic Mixed States[END_REF]. The growth of fluctuations over this new dynamics can be controlled relying only on (1.22); furthermore, the auxiliary dynamics and the fluctuation dynamics can be proved to be close in norm. Then, a bootstrap argument allows to control the growth of fluctuations on the original fluctuation dynamics, and ultimately to prove convergence of the many-body dynamics to the HFB evolution.

The paper is organized as follows. In Section 2 we review the formalism of second quantization, we recall useful estimates in Fock space, we introduce the Bogoliubov transformations and the quasi-free states. In Section 3 we state our main result, Theorem 3.3, about the derivation of the time-dependent HFB equation from the many-body dynamics in the semiclassical/mean-field regime. The rest of the paper is devoted to the proof of Theorem 3.3. In Section 4 we introduce the fluctuation dynamics, and prove a key identity for the growth of fluctuations, Proposition 4.1; in Proposition 4.3 we prove the propagation of the semiclassical structure (1.22), and we use it to control the growth of fluctuations on the auxiliary fluctuation dynamics, Proposition 4.5. Then, in Proposition 4.6 we prove norm-closeness of the auxiliary fluctuation dynamics and of the full fluctuation dynamics, and we use this result to prove a bound on the growth of fluctuations on the full fluctuation dynamics, Proposition 4.7. In Section 5 we put everything together and we prove Theorem 3.3. Finally, in Appendix A we prove the validity of the semiclassical structure on translation-invariant quasi-free states with low energy.

2 Fock space, Bogoliubov transformations and quasi-free states

Fock space

Let h = L 2 (R 3 × {1, . . . , M }), with M ∈ N, be the single-particle Hilbert space. Let F be the fermionic Fock space over h,

F = C ⊕ n≥1 h ∧n .
Vectors in the Fock space correspond to infinite sequences of functions (ψ (n) ) n with ψ (n) ∈ h ∧n . We shall use the notation

ψ (n) ≡ ψ (n) (x 1 , . . . , x n ) , (2.1)
where x = (x, σ) with x ∈ R 3 and σ ∈ {1, . . . , M }. A simple example of vector in F is the vacuum state, Ω = (1, 0, . . . , 0, . . .). Given ψ ∈ F, ψ = (ψ (0) , ψ (1) , . . . , ψ (n) , . . .), we introduce the Fock space scalar product as:

⟨ψ 1 , ψ 2 ⟩ = n≥0 ⟨ψ (n) 1 , ψ (n) 2 ⟩ h ⊗n = n≥0 dx 1 . . . dx n ψ (n) 1 (x 1 , . . . , x n )ψ (n) 2 (x 1 , . . . , x n ) ,
where we recall the notation dx f (x) := M σ=1 R 3 dx f (x, σ). Equipped with this natural inner product, the Fock space F is a Hilbert space. We denote by ∥ • ∥ the norm induced by this inner product, and with a slight abuse of notation we shall use the same notation to denote the operator norm of linear operators acting on F.

It is convenient to introduce creation and annihilation operators on F. Let f ∈ h. We define the creation operator a * (f ) and the annihilation operator a(f ) as

(a * (f )ψ) (n) (x 1 , . . . , x n ) := 1 √ n n j=1 (-1) j-1 f (x j )ψ (n-1) (x 1 , . . . , x j-1 , x j+1 , . . . , x n ) (a(f )ψ) (n) (x 1 , . . . , x n ) := √ n + 1 dx f (x) ψ (n+1) (x, x 1 , . . . , x n ) ,
for any ψ ∈ F. These definitions are supplemented by the requirement a(f )Ω = 0. It turns out that a * (f ) = a(f ) * . Also, one can check that the canonical anticommutation relations hold true:

{a(f ), a(g)} = {a * (f ), a * (g)} = 0 , {a(f ), a * (g)} = ⟨f, g⟩ h .
These relations imply that ∥a(f

)∥ ≤ ∥f ∥ 2 , ∥a * (f )∥ ≤ ∥f ∥ 2 .
We will also make use of operator-valued distributions a * x and a x , in terms of which:

a * (f ) = dx f (x) a * x , a(f ) = dx f (x) a x .
The creation and annihilation operators can be used to define the second quantization of observables. To begin, consider the number operator N , acting on a given Fock space vector as (N ψ) (n) = nψ (n) . In terms of the operator-valued distributions, it can be written as:

N = dx a * x a x .
The average number of particles over a state ψ is then:

N = ⟨ψ, N ψ⟩ . (2.2)
More generally, for a given one-particle operator O on h we define its second quantization dΓ(O) as the operator on the Fock space acting as follows:

dΓ(O) ↾ F (n) = n j=1 O (j)
where 1) . If O has the integral kernel O(x; y), we can represent dΓ(O) as:

O (j) = 1 ⊗(n-j) ⊗ O ⊗ 1 ⊗(j-
dΓ(O) = dxdy O(x; y)a * x a y .
In the next lemma we collect some bounds for the second quantization of one-particle operators, that will play an important role in the rest of the paper. We refer the reader to e.g.

[11, Lemma 3.1] for the proof.

Lemma 2.1 (Estimates for Fock space operators.). Let O be a bounded operator on h. We have, for any ψ ∈ F:

|⟨ψ, dΓ(O)ψ⟩| ≤ ∥O∥ op ⟨ψ, N ψ⟩ , ∥dΓ(O)ψ∥ ≤ ∥O∥ op ∥N ψ∥ .
Let O be a Hilbert-Schmidt operator. We have, for any ψ ∈ F:

∥dΓ(O)ψ∥ ≤ ∥O∥ HS N 1/2 ψ dxdx ′ O(x; x ′ )a x a x ′ ψ ≤ ∥O∥ HS N 1/2 ψ dxdx ′ O(x; x ′ )a * x a * x ′ ψ ≤ 2∥O∥ HS (N + 1) 1/2 ψ .
Let O be a trace class operator. We have, for any ψ ∈ F:

∥dΓ(O)ψ∥ ≤ 2 ∥O∥ tr ∥ψ∥ dxdy O(x; x ′ )a x a x ′ ψ ≤ 2 ∥O∥ tr ∥ψ∥ dxdy O(x; x ′ )a * x a * x ′ ψ ≤ 2 ∥O∥ tr ∥ψ∥ .
Given a state ψ ∈ F, we define its reduced one-particle density matrix as the non-negative trace class operator γ

N on h with integral kernel:

γ (1) N (x; y) = ⟨ψ, a * y a x ψ⟩ . (2.3) Thus, by (2.2), N = tr γ (1)
N . Furthermore, given a one-particle observable O, we have:

⟨ψ, dΓ(O)ψ⟩ = dxdy O(x; y) ⟨ψ, a * x a y ψ⟩ = tr Oγ (1)
N .

(2.4) Also, we define the pairing matrix as the operator π

N on h with integral kernel: π

N (x; y) = ⟨ψ, a y a x ψ⟩ .

(2.5)

In the special case in which the vector ψ in an N -particle state, N ψ = N ψ, the pairing operator π

N is trivially vanishing. In general, if the state ψ does not have a definite number of particles, as it will be the case in the present work, the pairing operator will be nonzero. It is not difficult to see that the reduced one-particle density matrix and the pairing operator satisfy:

0 ≤ γ (1) N ≤ 1 , ∥π (1)
N ∥ ≤ 1 , π (1)
N (x; y) = -π (1)
N (y; x) .

(2.6)

Next, we lift the many-body Hamiltonian to the Fock space, as follows. We define the Hamiltonian in Fock space as

H N = 0 ⊕ n≥1 H (n)
N , where

H (n) N = n i=1 -ε 2 ∆ i + 1 N n i<j V (x i -x j ) ,
with V a real-valued and even function (further regularity assumptions will follow). In terms of the operator-valued distributions, we can rewrite the Fock space Hamiltonian as:

H N = ε 2 dx ∇a * x ∇a x + 1 2N dxdy V (x -y)a * x a * y a y a x . (2.7)
Finally, the time evolution of a state in the Fock space is defined as ψ t = e -iH N t/ε ψ. That is, each sector evolves according to the n-particle Schrödinger equation,

iε∂ t ψ (n) t = H (n) N ψ (n) t , ψ (n) 0 = ψ (n) . (2.8)

Bogoliubov transformations

Here we shall give a brief introduction to the (fermionic) Bogoliubov transformations, which will play a crucial role in our analysis. We refer the reader to e.g. [START_REF] Solovej | Many Body Quantum Mechanics[END_REF][START_REF] Dereziński | Mathematics of quantization and quantum fields[END_REF] for a more detailed introduction to the topic. Given f, g ∈ h, let us define the operator A(f, g) as:

A(f, g) := a(f ) + a * (g) , A * (f, g) := (A(f, g)) * = a * (f ) + a(g) .
(2.9)

That is:

A * (f, g) = A(g, f ) . (2.10)
These operators can be viewed as generalizing the notion of creation and annihilation operators, to which they reduce for g = 0. They satisfy the anticommutation relations:

{A(f 1 , g 1 ), A * (f 2 , g 2 )} = ⟨(f 1 , g 1 ), (f 2 , g 2 )⟩ h⊕h . (2.11) Notice that {A(f 1 , g 1 ), A(f 2 , g 2 )}, {A * (f 1 , g 1 ), A * (f 2 , g 2 )} do not vanish, in general. A fermionic Bogoliubov transformation is a linear map ν : h ⊕ h → h ⊕ h, such that, for all f 1 , f 2 , g 1 , g 2 in h: A(ν(f 1 , g 1 )), A * (ν(f 2 , g 2 )) = A(f 1 , g 1 ), A * (f 2 , g 2 ) , (2.12)
and such that, for all f, g in h:

A * (ν(f, g)) = A(ν(g, f )) . (2.13)
In other words, the operators B(f, g) = A(ν(f, g)) satisfy the same anticommutation relations of A(f, g). The condition (2.12) is equivalent to ν * ν = 1, while the condition (2.13) is equivalent to, calling C the complex conjugation operator:

νC = Cν . (2.14)
Then, it is not difficult to check that ν is a Bogoliubov transformation if and only if it has the form:

ν = u v v u , (2.15)
where u, v are linear operators on h such that:

u * u + v * v = 1, u * v + v * u = 0 uu * + vv T = 1, uv * + vu T = 0 , (2.16)
with O = COC. We say that a Bogoliubov transformation is implementable if there exists a unitary operator R on F such that:

R * A(f, g)R = A(ν(f, g)) ,
for all f, g ∈ h.

(2.17)

It is well-known that a Bogoliubov transformation is implementable if and only if v is a Hilbert-Schmidt operator; this fact goes under the name of Shale-Stinespring condition. The operator R is called the implementor of the Bogoliubov transformation; with a slight abuse of terminology, in the following we will refer to R as the Bogoliubov transformation. We shall extensively use the following transformation property of the operator-valued distributions:

R * a x R = a(u x ) + a * (v x ) , R * a * x R = a * (u x ) + a(v x ) , (2.18)
where u x := u(•; x). Furthermore, it is not difficult to see that R * is the implementor of the Bogoliubov transformation associated with ν * ≡ ν -1 . In particular, the following relations hold:

Ra x R * = a(u * x ) + a * (v * x ), Ra * x R * = a * (u * x ) + a(v * x ) . (2.19)

Quasi-free states and Hartree-Fock-Bogoliubov equation

Given a normalized state ψ ∈ F, we define its one-particle density matrix Γ ψ as the operator on h ⊕ h such that:

⟨(f 1 , g 1 ), Γ ψ (f 2 , g 2 )⟩ h⊕h = ⟨ψ, A * (f 2 , g 2 )A(f 1 , g 1 )ψ⟩ F . (2.20)
A simple computation shows that Γ ψ can be expressed in terms of the one-particle reduced density matrix and on the pairing matrix as:

Γ ψ = γ (1) N π (1) N -Cπ (1)
N C 1 -Cγ (1) N C , (2.21)
where we used that π

N is antisymmetric, recall Eq. (2.6). Pure quasi-free states correspond to vectors in the Fock space of the form RΩ, with R a Bogoliubov transformation. Suppose that ⟨RΩ, N RΩ⟩ = N . We shall denote by ω N and by α N the one-particle density matrix and the pairing matrix of RΩ, respectively:

Γ RΩ = ω N α N -α N 1 -ω N . (2.22)
The density matrix of RΩ can be easily represented in terms of the operators u and v entering in the definition of ν, recall Eq. (2.15). We have:

ω N = v * v , α N = v * u . (2.23)
As a sanity check, notice that ω N ≤ 1 as a consequence of the first of (2.16), while α T N = -α N :

α T N = α * N = u * v = u * v = -v * u = -α N , (2.24)
where the fourth equality follows from the second of (2.16). One can also check that the density matrix of RΩ is a projection:

Γ RΩ = Γ 2 RΩ . (2.25)
In particular, the following relations hold:

ω N (1 -ω N ) = α N α * N , ω N α N -α N ω N = 0 . (2.26)
Quasi-free states are completely determined by their one-particle density matrix: all correlation functions can be computed starting from the one-particle density matrix via the Wick rule [START_REF] Solovej | Many Body Quantum Mechanics[END_REF]. Let us use the notation a ♯ to denote the operators a and a * . It turns out that:

RΩ, a ♯ 1 x 1 a ♯ 2 x 2 • • • a ♯ 2j x 2j RΩ = π∈P 2j sgn(π) RΩ, a ♯ σ(1)
x σ(1) a

♯ σ(2) x σ(2) RΩ • • • RΩ, a ♯ σ(2j-1)
x σ(2j-1) a

♯ σ(2j)
x σ(2j) RΩ ,

(2.27) where P 2j is the set of pairings:

P 2j = σ ∈ S 2j σ(2ℓ -1) < σ(2ℓ + 1) , ℓ = 1, . . . , j -1 , σ(2ℓ -1) < σ(2ℓ) , ℓ = 1, . . . , j , (2.28
) with S 2j the set of permutations of 2j elements. The identity (2.27) can be checked starting from the relation (2.18), combined with the canonical anticommutation relations and the property a x Ω = 0.

We are interested in the many-body evolution of initial data close to quasi-free states, generated by the Hamiltonian (2.7). It turns out that if one assumes that the many-body evolution preserves quasi-freeness for all times, one finds the following coupled evolution equations for the one-particle density matrix:

iε∂ t ω N,t = h HF (t), ω N,t + Π t α * N,t -α N,t Π * t , iε∂ t α N,t = h HF (t)α N,t + α N,t h HF (t) + Π t 1 -ω N,t -ω N,t Π t , (2.29)
where h HF (t) is the Hartree-Fock Hamiltonian,

h HF (t) = -ε 2 ∆ + ρ t * V -X t , (2.30)
with ρ t the normalized density and X t , Π t linear operators on h,

ρ t (x) = 1 N ω N,t (x; x) , X t (x; y) = 1 N V (x -y)ω N,t (x; y) , Π t (x, y) := 1 N V (x -y)α N,t (x, y) .
(2.31) Thus, Eqs. (2.29) define the evolution of the one-particle density matrix:

Γ N,t = ω N,t α N,t α * N,t 1 -ω N,t . (2.32)
It is not difficult to check that Eqs. (2.29) are equivalent to:

iε∂ t Γ N,t = [H(t), Γ N,t ] , (2.33)
with:

H(t) = h HF (t) Π t Π * t -h HF (t) . (2.34)
That is, H(t) generates a unitary dynamics for Γ N,t . Hence, since Γ N = Γ 2 N , recall Eq. (2.25), we also have Γ N,t = Γ 2 N,t . Eq. (2.33) is called the Hartree-Fock-Bogoliubov equation, and it describes interesting physical phenomena, such as the time evolution of superconductors. The well-posedness of the HFB equation (2.29) for a class of potentials that includes those considered in the present work, is well-known; see e.g. [START_REF] Benedikter | The Dirac-Frenkel Principle for Reduced Density Matrices, and the Bogoliubov-de-Gennes Equations[END_REF]. The main goal of the present work will be to give a rigorous derivation of the HFB equation starting from the full manybody quantum dynamics.

Main result

Before stating our main result, let us introduce the assumptions on the initial data that we will consider in this work. Assumption 3.1 (Assumptions on the initial datum). Let Γ N : h ⊕ h → h ⊕ h be the one-particle density matrix of a pure, quasi-free state:

Γ N = ω N α N α * N 1 -ω N , 0 ≤ Γ N ≤ 1 , Γ N = Γ 2 N . (3.1)
We shall suppose that:

tr ω N = N , ∥α N ∥ HS ≤ CN 1/3 . (3.2)
Furthermore, we shall suppose that ω N , α N satisfy the following commutator estimates:

sup p∈R 3 1 1 + |p| ∥[ω N , e ip•x ]∥ HS ≤ CN 1/3 , ∥[ω N , ε∇]∥ HS ≤ CN 1/3 , ∥[α N , ε∇]∥ HS ≤ CN 1/3 . (3.3)
Remark 3.2. (i) As we will discuss in Appendix A, all the above assumptions are satisfied for translation-invariant quasi-free states on a finite torus, at low energy.

(ii) We shall refer to the commutator estimates in (3.3) as the semiclassical structure of the initial datum.

The next theorem is the main result of the paper.

Theorem 3.3 (Main result). Let V ∈ L 1 (R 3 ) such that V (x) ∈ R , V (x) = V (-x) , dp | V (p)|(1 + |p| 2 ) ≤ C . (3.4)
Let ω N , α N be the reduced one-particle density matrix and the pairing matrix of a quasifree state satisfying Assumption 3.1. Let ψ = R 0 ξ with R 0 the Bogoliubov transformation associated to (ω N , α N ) and ξ ∈ F, ∥ξ∥ = 1, such that, for some k ∈ N:

⟨ψ, N ψ⟩ = N , ⟨ξ, N α(k) ξ⟩ ≤ C k (3.5)
with α(k) = (13/2)k if k is even and α(k) = (13/2)k + 3/2 if k is odd. Let ψ t = e -iH N t/ε ψ be the solution of the many-body Schrödinger equation, and let R t be the Bogoliubov transformation associated with the solution ω N,t , α N,t of the Hartree-Fock-Bogoliubov equation (2.29).

Then, for all t ≥ 0, and for all j ≤ k:

⟨ψ t , a ♯ 1 x 1 . . . a ♯ 2j x 2j ψ t ⟩ = ⟨R t Ω, a ♯ 1 x 1 . . . a ♯ 2j x 2j R t Ω⟩ + E ♯ 1 ,...,♯ 2j 2j;t (x 1 , . . . , x 2j ) (3.6)
where a ♯ x can be either a x or a * x , a ♯ 1 x 1 . . . a

♯ 2j
x 2j is a normal-ordered monomial, and

dx 1 . . . dx 2j E ♯ 1 ,...,♯ 2j 2j;t (x 1 , . . . , x 2j ) 2 ≤ C(j) exp(d(j) exp c|t|)N j-2/3-max(1/3, p/6) , (3.7)
where the constants C(j), d(j) depend on j and are independent of N and of t, while p is equal to the absolute value of the difference of the number of creation and annihilation operators in the left-hand side of (3.6).

Remark 3.4. Theorem 3.3 proves convergence of the many-body Schrödinger equation to the Hartree-Fock-Bogoliubov equation, in the sense of density matrices, in the Hilbert-Schmidt topology. In fact, the left-hand side of (3. ). The Hilbert-Schmidt norm squared of the first term in the right-hand side of (3.6) is much larger than the right-hand side of (3.7) for all times of order 1 in N . This can be seen writing it in terms of the one-particle density matrix using the Wick rule. For instance:

⟨R t Ω, a * x 1 a * x 2 a x 3 a x 4 R t Ω⟩ = ω N,t (x 3 ; x 2 )ω N,t (x 4 ; x 1 ) -ω N,t (x 3 ; x 1 )ω N,t (x 4 ; x 2 ) + α N,t (x 1 ; x 2 )α N,t (x 4 ; x 3 ) . (3.8)
Assuming that the bounds on ∥ω N,t ∥ HS and ∥α N,t ∥ HS are tight, namely ∥α N,t ∥ 2 HS ∼ N 2/3 and ∥ω N,t ∥ 2 HS ∼ N , the largest contribution to the two-particle density matrix of the quasi-free state comes from the term with the largest number of ω N,t . For gauge-invariant correlation functions, i.e. with an equal number of creation and annihilation operators, the Hilbert-Schmidt norm of the largest contribution is of order N j . In general, if m is the number of annihilation operators, say smaller than the number of creation operators m + p, the largest term has m = j -p/2 factors ω N,t and p/2 factors α N,t . Thus, its Hilbert-Schmidt norm squared is of order N j-p/6 , which is much larger than the estimate for the error term in (3.7).

Remark 3.5 (Resolution of the pairing operators).

It is interesting to compare the error term with the smallest contribution to the quasi-free approximation. If the error is subleading, this means that the quasi-free approximation can be fully resolved. Let us focus on gaugeinvariant correlation functions (p = 0). The smallest Wick pairing in Hilbert-Schmidt norm is the one that has the largest number of contractions of "aa" or "a * a * " monomials, which correspond to kernels of α N,t operators (or its complex conjugate). The number of such contractions is less or equal than j -1 if the number of "a" operators is odd, and it is less or equal than j if the number of "a" operators is even. In the p = 0 case, the number of "a" operators is exactly j, therefore the Hilbert-Schmidt norm squared of the Wick pairing with the largest number of α-contractions is:

(N 2/3 ) j-1 × N if j is odd (N 2/3 ) j if j is even. (3.9)
Comparing these estimates with the bound (3.7) for p = 0, we have that, for N ≫ 1:

(N 2/3 ) j-1 × N ≫ ∥E ♯ 1 ,...,♯ 2j 2j;t ∥ 2 HS if j = 1, 3 (N 2/3 ) j ≫ ∥E ♯ 1 ,...,♯ 2j 2j;t ∥ 2 HS if j = 2.
(3.10)

Therefore, all Wick pairings are individually much bigger than the error term, and hence they can be resolved, for the j-body reduced density matrices up to j = 3. More generally, it is clear that if p > 0 the bound (3.7) also allows to fully resolve correlations of higher order.

Remark 3.6 (Comparison with the Hartree-Bogoliubov dynamics.). As mentioned in the Introduction, the proof of Theorem 3.3 could be simplified if one replaces in the assumptions (3.2), (3.3) the Hilbert-Schmidt norm with the square root of the trace norm. Furthermore, as for the case of N -particle states [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF], under these stronger assumptions it is possible to prove that neglecting the exchange operators X t , Π t in the Hartree-Fock-Bogoliubov equations (2.29) does not deteriorate the error estimate (3.7) in the main result. Thus, for this class of initial data, we cannot distinguish the Hartree-Fock-Bogoliubov dynamics from the Hartree-Bogoliubov dynamics. We prefer to work under weaker assumptions that we can verify in concrete cases, see Appendix A. For this class of initial data, we do not know how to prove the indistinguishability of the Hartree-Fock-Bogoliubov and Hartree-Bogoliubov dynamics.

The rest of the paper is devoted to the proof of Theorem 3.3.

The fluctuation dynamics

Suppose that (ω N,t , α N,t ) solves the time-dependent Hartree-Fock-Bogoliubov equation (2.29).

Let us introduce the fluctuation dynamics as:

U N (t; s) := R * t e -iH(t-s)/ε R s . (4.1)
The main goal of this section will be to prove an estimate for the evolution of powers of the number operators N under the fluctuation dynamics. This will be the key technical tool behind the proof of Theorem 3.3.

Growth of fluctuations

The starting points is the following proposition.

iε d dt U N (t; 0)ξ, (N + 1) k U N (t; 0)ξ = - 4i N k j=1 Im dxdy V (x -y) × U N (t; 0)ξ, (N + 1) j-1 a * (u t,x )a(v t,y )a(u t,y )a(u t,x )(N + 1) k-j U N (t; 0)ξ + U N (t; 0)ξ, (N + 1) j-1 a(v t,x )a(v t,y )a(u t,y )a(u t,x )(N + 1) k-j U N (t; 0)ξ + U N (t; 0)ξ, (N + 1) j-1 a * (u t,y )a * (v t,y )a * (v t,x )a(v t,x )(N + 1) k-j U N (t; 0)ξ . (4.2)
Remark 4.2. The same identity holds true for the Hartree-Fock case [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF], see Proposition 3.3 there. Here, however, the operators u t and v t are not orthogonal.

Proof. To begin, recall the identity:

R t N R * t = N -2dΓ(ω N,t ) + N + dxdy a * x a * y α N,t (y, x) + a x a y α N,t (x, y) , (4.3)
where we used that tr ω N,t = N . Therefore, since U N (t; s) = R * t e -iH(t-s)/ε R s , we have:

iε d dt U N (t; 0) * N U N (t; 0) = iε d dt R * 0 e iHt/ε R t N R * t e -iHt/ε R 0 = iε d dt R * 0 e iHt/ε -2dΓ(ω N,t ) + dxdya * x a * y α N,t (y, x) + h.c. e -iHt/ε R 0 = R * 0 e iHt/ε -2dΓ(iε∂ t ω N,t ) + dxdya * x a * y (iε∂ t α N,t )(y, x) -h.c. e -iHt/ε R 0 + R * 0 e iHt/ε 2 H, dΓ(ω N,t ) - dxdy H, a * x a * y α N,t (y, x) -h.c. e -iHt/ε R 0 (4.4)
where we used that the first and the third term in (4.3) do not contribute to the time derivative. Then, by unitarity of the Bogoliubov transformation.

iε d dt U N (t; 0) * N U N (t; 0) = 2U N (t; 0) * R * t H, dΓ(ω N,t ) -dΓ(iε∂ t ω N,t ) R t U N (t; 0) + U N (t; 0) * R * t dxdy a * x a * y (iε∂ t α N,t )(y, x) -H, a * x a * y α N,t (y, x) -h.c. R t U N (t; 0) . (4.5)
In this expression, it is not difficult to see that all terms involving the Laplacian cancel out exactly. Next, let us consider the commutators involving the many-body interaction. We compute:

2 dxdy V, a * x a y ω N,t (x, y) = 2 N dxdydz V (x -y)ω N,t (z, y) a * z a * x a y a x -h.c., (4.6)
and

dxdy V, a * x a * y α N,t (y, x) = 1 N dxdy V (x -y)α N,t (y, x) a * x a * y + 2 N dxdydz V (x -y)α N,t (z, y) a * x a * y a * z a x . (4.7)
Therefore, focusing on the following term,

A := R * t 2 V, dΓ(ω N,t ) - dxdy V, a * x a * y α N,t (y, x) -h.c. R t , (4.8)
by Eqs. (4.6), (4.7), we have:

A = 2 N dxdy V (x -y)R * t dz ω N,t (z, y)a * z a * x a y a x -α N,t (z, y)a * x a * y a * z a x R t - 1 N dxdy V (x -y)R * t α N,t (y, x)a * x a * y R t -h.c.. (4.9)
Let us denote by A 2 and A 4 the quadratic and the quartic terms in Eq. (4.9),

A 2 := - 1 N dxdy V (x -y)α N,t (y, x)R * t a * x a * y R t -h.c. A 4 := 2 N dxdy V (x -y)R * t dz ω N,t (z, y)a * z a * x a y a x -α N,t (z, y)a * x a * y a * z a x R t -h.c.. ( 4

.10)

Consider A 4 . We rewrite it as:

A 4 = 2 N dxdy V (x -y)R * t a * x dz -ω N,t (z, y)a * z + α N,t (z, y)a z a y a x R t -h.c. = 2 N dxdy V (x -y)R * t a * x R t • dz -ω N,t (z, y) a * (u t,z ) + a(v t,z ) + α N,t (z, y) a(u t,z ) + a * (v t,z ) R * t a y a x R t -h.c..
Performing the z integration we obtain:

dz -ω N,t (z, y) a * (u t,z ) + a(v t,z ) + α N,t (z, y) a(u t,z ) + a * (v t,z ) = dz dr a * r -u t (r, z)ω N,t (z, y) + v t (r, z)α N,t (z, y) + a r -v t (r, z)ω N,t (z, y) + u t (r, z)α N,t (z, y) = dr a * r v t α N,t -u t ω N,t (r, y) + a r u t α N,t -v t ω N,t (r, y) = -a(v t,y ) . (4.11)
To prove the last identity we used that:

v t α N,t -u t ω N,t = v t v T t u t -u t v * t v t = v t v T t u t + v t u T t v t = v t (v T t u t + u T t v t ) = 0 , (4.12)
where: the first step follows from the representation of ω N,t and α N,t in terms of u t and v t , Eq. (2.23); the second step follows from the second property in the second line of (2.16); and the third step follows from the second property in the first line of (2.16), after taking a complex conjugate. Eq. (4.12) proves the vanishing of the creation contribution to the right-hand side of (4.11). To obtain the last identity in (4.11), we used:

u t α N,t -v t ω N,t = u t v T t u t -v t v * t v t = -v t u * t u t -v t v * t v t = -v t (u * t u t + v * t v t ) = -v t . (4.
13) The first step follows from the representation of ω N,t and α N,t in terms of u t and v t ; the second step follows from the second property in the second equation in (2.16), after taking a complex conjugate; the last step follows from the first property in the first equation of (2.16). This concludes the check of (4.11).

Let us now go back to A 4 in (4.11). Using (4.11), we have:

A 4 = - 2 N dxdy V (x -y)R * t a * x R t a(v t,y )R * t a y a x R t -h.c. = - 2 N dxdy V (x -y)(a * (u t,x ) + a(v t,x ))a(v t,y )(a(u t,y ) + a * (v t,y ))(a(u t,x ) + a * (v t,x )) -h.c., (4.14)
and expanding the products: The second and the fifth term in the right-hand side of (4.15) are self-adjoint, and hence cancel with minus the hermitian conjugate. The sum of the first and the sixth term is also self-adjoint (each one is the adjoint of the other), so they cancel as well when subtracting the hermitian conjugate. In order to proceed, we put the remaining terms in the right-hand side into normal order. We get: Using that V (x -y) = V (y -x), the second and the fourth term are self-adjoint, and hence they cancel with minus their hermitian conjugates. Then, we rewrite (4.16) as:

A 4 = - 2 N dxdy V (x -y) a * (u t,x )a(v t,
A 4 = - 2 N dxdy V (x -y) a * (u t,x )a(v t,
A 4 = - 2 N dxdy V (x -y) a * (u t,x )a(v t,y )a(u t,y )a(u t,x ) + a(v t,x )a(v t,y )a(u t,y )a(u t,x ) -a * (v t,x )a(v t,x )a(v t,y )a(u t,y ) -2 dxdy (a * (u t,x ) + a(v t,x ))Π t (x; y)a(v t,y ) + 2 dxdy a(v t,x )X t (x, y)a(u t,y ) -2 dy a(v t,y )(ρ t * V )(y)a(u t,y ) -h.c.. (4.17)
We claim that the quadratic terms appearing in (4.17) cancel out with A 2 and with the quadratic terms in Eq. (4.5). Consider the quadratic terms in Eq. (4.5), minus the kinetic term. The sum of these terms, plus the term A 2 in Eq. (4.10), gives:

-R * t dΓ(C t ) + dxdy a * x a * y D t (x; y) R t -h.c. (4.18)
where the operators C t and D t are:

C t = ρ t * V -X t , ω N,t + Π t α * N,t -α N,t Π * t D t = ρ t * V -X t α N,t + α N,t ρ t * V -X t -Π t ω N,t -ω N,t Π t ; (4.19)
these operators take into account the interaction part of the equations (2.29), after subtracting the linear term in α N,t (thanks to A 2 ). Let us now plug these expressions in (4.18). Consider the contribution to (4.18) due to ρ t * V . We have:

-R * t dxdy a * x (ρ t * V )(x) -(ρ t * V )(y) ω N,t (x; y)a y + dxdy a * x (ρ t * V )(x) + (ρ t * V )(y) α N,t (x; y)a * y R t -h.c. = -2R * t dxdy a * x (ρ t * V )(x)ω N,t (x; y)a y + dxdy a * x (ρ t * V )(x)α N,t (x; y)a * y R t -h.c. (4.20)
Acting with the Bogoliubov transformation, we get:

(4.20) = -2 dxdy (a * (u t,x ) + a(v t,x ))(ρ t * V )(x)ω N,t (x; y)(a(u t,y ) + a * (v t,y )) -2 dxdy (a * (u t,x ) + a(v t,x ))(ρ t * V )(x)α N,t (x; y)(a * (u t,y ) + a(v t,y )) -h.c. = -2 dx (a * (u t,x ) + a(v t,x ))(ρ t * V )(x) a((u t ω t ) x ) + a * ((v t ω t ) x ) + 2 dx (a * (u t,x ) + a(v t,x ))(ρ t * V )(x) a * ((u t α t ) x ) + a((v t α t ) x ) -h.c.. (4.21)
Using the relations (4.12), (4.13) we obtain:

(4.20) = -2 dx (a * (u t,x ) + a(v t,x ))(ρ t * V )(x)a * (v t,x ) -h.c. = -2 dx a * (u t,x )(ρ t * V )(x)a * (v t,x ) -h.c. . ( 4
R * t dΓ(X t ω t ) + dxdy a * x a * y (X t α t )(x; y) R t -R * t dΓ(ω t X t ) -dxdy a * x a * y (α t X t )(x; y) R t -h.c. ≡ 2R * t dΓ(X t ω t ) + dxdy a * x a * y (X t α t )(x; y) R t -h.c. , (4.23)
where we used the antisymmetry of α N,t . We have:

2R * t dΓ(X t ω t ) + dxdy a * x a * y (X t α t )(x; y) R t = 2 dxdy (a * (u t,x ) + a(v t,x ))(X t ω N,t )(x; y)(a(u t,y ) + a * (v t,y )) + 2 dxdy (a * (u t,x ) + a(v t,x ))(X t α N,t )(x; y)(a * (u t,y ) + a(v t,y )) = 2 dxdz (a * (u t,x ) + a(v t,x ))X t (x; z)(a((u t ω t ) z ) + a * ((v t ω t ) z ) -2 dxdy (a * (u t,x ) + a(v t,x ))X t (x; z)(a * ((u t α t ) z ) + a((v t α t ) z ) . (4.24)
Using the relations (4.12), (4.13), we get

2R * t dΓ(X t ω t ) + dxdy a * x a * y (X t α t )(x; y) R t -h.c. = 2 dxdz (a * (u t,x ) + a(v t,x ))X t (x; z)a * (v t,z ) -h.c. = 2 dxdz a * (u t,x )X t (x; z)a * (v t,z ) -h.c. . (4.25)
Finally, consider the terms corresponding to Π t in (4.19). Their contribution to (4.18) is:

-

R * t dΓ(Π t α * t -α t Π * t )R t + R * t dxdy a * x a * y Π t ω t + ω t Π t (x; y)R t -h.c. ≡ 2R * t dΓ(Π t α t )R t + 2R * t dxdy a * x a * y Π t ω t (x; y)R t -h.c. , (4.26)
where we used the antisymmetry of α N,t and of Π t . We have:

2R * t dΓ(Π t α t ) + dxdy a * x a * y Π t ω t (x; y) R t -h.c. = 2 dxdy (a * (u t,x ) + a(v t,x )) Π t α t (x; y)(a(u t,y ) + a * (v t,y )) + 2 dxdy (a * (u t,x ) + a(v t,x )) Π t ω t (x; y)(a * (u t,y ) + a(v t,y )) = -2 dxdz (a * (u t,x ) + a(v t,x ))Π t (x; z)(a((u t α t ) z ) + a * ((v t α t ) z )) + 2 dxdz (a * (u t,x ) + a(v t,x ))Π t (x; z)(a * ((u t ω t ) z ) + a((v t ω t ) z )) . (4.27)
Using again the relations (4.12), (4.13), we get:

2R * t dΓ(Π t α t ) + dxdy a * x a * y Π t ω t (x; y) R t -h.c. = 2 dxdz (a * (u t,x ) + a(v t,x ))Π t (x; z)a(v t,z ) -h.c. (4.28)
Putting together (4.22), (4.25), (4.28), we obtained:

-R * t dΓ(C t ) + dxdy a * x a * y D t (x; y) R t -h.c. = -2 dx a * (u t,x )(ρ t * V )(x)a * (v t,x ) + 2 dxdz a * (u t,x )X t (x; z)a * (v t,z ) + 2 dxdz (a * (u t,x ) + a(v t,x ))Π t (x; z)a(v t,z ) -h.c. (4.29)
Thus, as claimed, we found that the sum of the quadratic terms in (4.17) and in (4.29) is vanishing. This implies:

iε d dt U N (t; 0) * N U N (t; 0) = - 1 N dxdy V (x -y)U N (t; 0) * a * (u t,x )a(v t,y )a(u t,y )a(u t,x ) + a(v t,x )a(v t,y )a(u t,y )a(u t,x ) -a * (v t,x )a(v t,x )a(v t,y )a(u t,y ) U N (t; 0) -h.c. (4.30)
The final statement, Eq. (4.2), follows from the identity:

iε d dt U N (t; 0) * (N + 1) k U N (t; 0) = k j=1 U N (t; 0) * (N + 1) j-1 U N (t; 0) × U N (t; 0) * iε d dt U N (t; 0) * N U N (t; 0) U N (t; 0) * (N + 1) k-j U N (t; 0) . (4.31)
This concludes the proof of Proposition 4.1.

Propagation of the semiclassical structure

In order to bound the growth of fluctuations, we will make crucial use of the propagation of the semiclassical structure of Assumption 3.1 along the flow of the HFB equation. Let (ω N,t , α N,t ) be the solution of the Hartree-Fock-Bogoliubov equation (2.29). Then, for all times t ≥ 0: Proof. The proof is based on an extension of [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF] from the Hartree-Fock equation to the Hartree-Fock-Bogoliubov equation. To begin, let us prove the invariance of the trace-norm of ω N,t . Since 0 ≤ ω N ≤ 1, we have:

sup p∈R 3 1 1 + |p| ∥[ω N,t , e ip•x ]∥ HS ≤ Ce c|t| (N 1/3 + |t|∥α N ∥ HS ) ∥[ω N,t , ε∇]∥ HS ≤ Ce c|t| (N 1/3 + |t|∥α N ∥ HS ) ∥[α N,t , ε∇]∥ HS ≤ Ce c|t| (N 1/3 + |t|∥α N ∥ HS ) . ( 4
∥ω N,t ∥ 2 HS ≤ tr ω N,t . (4.34)
Then, from (2.29):

iε∂ t tr ω N,t = tr Π t α * N,t -α N,t Π * t (4.35)
where we used that the trace of the commutator is zero. Next, we rewrite:

tr Π t α * N,t = dxdy Π t (x; y)α * N,t (y; x) = - 1 N dxdy V (x -y)|α N,t (x; y)| 2 (4.36)
thanks to the antisymmetry of α N,t . Being the right-hand side of Eq. (4.35) equal to the imaginary part of tr Π t α * N,t , and since the right-hand side of (4.36) is real, we immediately conclude that tr ω N,t is constant in time. This proves the first of (4.33).

Next, let us prove the propagation of the commutator estimates. Let us start by computing, using the HFB equation for ω N,t :

iε d dt e ip•x , ω N,t = e ip•x , h HF (t), ω N,t + e ip•x , Π t α * N,t -α N,t Π * t = h HF (t), e ip•x , ω N,t + ω N,t , h HF (t), e ip•x + e ip•x , Π t α * N,t -α N,t Π * t = h HF (t), e ip•x , ω N,t -ω N,t , ε 2 ∆, e ip•x -ω N,t , X t , e ip•x + e ip•x , Π t α * N,t -α N,t Π * t , (4.37)
where in the second step we used the Jacobi identity, and in the third step we used that the direct term ρ t * V commutes with e ip•x . Consider the term involving the commutator with the Laplacian. Noticing that:

ε 2 ∆, e ip•x = iε∇ • εp e ip•x + e ip•x εp • iε∇ , (4.38)
we obtain:

ω N,t , ε 2 ∆, e ip•x = ω N,t , iε∇ • εp e ip•x + iε∇ • εp ω N,t , e ip•x + εp e ip•x • ω N,t , iε∇ + ω N,t , e ip•x εp • iε∇ . (4.39)
The terms involving the gradient operator outside of the commutator, which are unbounded, can be reabsorbed thanks to a suitable unitary conjugation, together with the first term in the right-hand side of the last line of (4.37). Consider the unitary groups U 1 (t; s) and U 2 (t; s), generated by the following differential equations:

iε∂ t U 1 (t, s) = A(t)U 1 (t, s) , iε∂ t U 2 (t, s) = B(t)U 2 (t, s) , U 1 (s; s) = U 2 (s; s) = 1 (4.40) where: A(t) = h HF (t) + iε 2 ∇ • p , B(t) = h HF (t) -iε 2 ∇ • p . (4.41)
We then get:

iε d dt U * 1 (t, 0) e ip•x , ω N,t U 2 (t, 0) = U * 1 (t, 0) -A(t) e ip•x , ω N,t + e ip•x , ω N,t B(t) + iε d dt e ip•x , ω N,t U 2 (t, 0) = U * 1 (t, 0) -ω N,t , iε∇ • εp e ip•x -εp e ip•x • ω N,t , iε∇ -ω N,t , X t , e ip•x + e ip•x , Π t α * N,t -α N,t Π * t U 2 (t, 0). (4.42)
Integrating in time this identity, and using ∥U * 1 (t, 0) e ip•x , ω N,t U 2 (t, 0)∥ HS = ∥ e ip•x , ω N,t ∥ HS , we find:

e ip•x , ω N,t HS ≤ e ip•x , ω N,0 HS + 2|p| t 0 ds ε∇, ω N,s HS + 1 ε t 0 ds ω N,s , X s , e ip•x HS + 1 ε t 0 ds e ip•x , Π s α * N,s -α N,s Π * s HS . (4.43)
Next, it is useful to rewrite the direct term V * ρ s , the exchange term X s and the term Π s as, for x = (x, σ), acting as (xf )(x, σ) = (x, σ)f (x, σ):

(V * ρ s )(x) = 1 N dq V (q)ρ s (q, σ)e iq•x X s = 1 N dq V (q)e iq•x ω N,s e -iq•x Π s = 1 N dq V (q)e iq•x α N,s e -iq•x , (4.44)
where ρs (q, σ) = dx e iq•x ω N,s (x, σ; x, σ). Using these representations, one can bound the third and the fourth terms in the right-hand side of Eq. (4.43) as, recalling that ∥ω N,s ∥ ≤ 1,

∥α N,s ∥ ≤ 1: ω N,s , X s , e ip•x HS ≤ 2 N dq | V (q)| e ip•x , ω N,s HS ≤ C N e ip•x , ω N,s HS e ip•x , Π s α * N,s -α N,s Π * s HS ≤ 4 N dq | V (q)| α N,s HS ≤ C N α N,s HS , (4.45)
where we used the assumptions (3.4) on the potential. Therefore, plugging these bounds in (4.43) we obtain:

e ip•x , ω N,t HS ≤ e ip•x , ω N,0 HS + C t 0 ds |p| ε∇, ω N,s HS + 1 N ε e ip•x , ω N,s HS + 1 N ε ∥α N,s ∥ HS . (4.46)
Consider now the commutator with the gradient in Eq. (4.46). We compute:

iε d dt ε∇, ω N,t = ε∇, h HF (t), ω N,t + Π t α * N,t -α N,t Π * t = h HF (t), ε∇, ω N,t + ω N,t , h HF (t), ε∇ + ε∇, Π t α * N,t -α N,t Π * t = h HF (t), ε∇, ω N,t + ω N,t , V * ρ t , ε∇ -ω N,t , X t , ε∇ + ε∇, Π t α * N,t -α N,t Π * t . (4.47)
In order to get rid of the first term, we conjugate the commutator with the unitary group U 3 (t; s), generated by the differential equation:

iε∂ t U 3 (t, s) = h HF (t)U 3 (t, s) , U 3 (s; s) = 1 . (4.48)
Then:

iε d dt U * 3 (t, 0) ε∇, ω N,t U 3 (t, 0) = U * 3 (t, 0) ω N,t , V * ρ t , ε∇ -ω N,t , X t , ε∇ + ε∇, Π t α * N,t -α N,t Π * t U 3 (t, 0) . (4.49)
Integrating in time, and using the invariance of the Hilbert-Schmidt norm under unitary conjugations, we obtain:

ε∇, ω N,t HS ≤ ε∇, ω N HS + 1 ε t 0 ds ω N,s , V * ρ s , ε∇ HS + 1 ε t 0 ds ω N,s , X s , ε∇ HS + 1 ε t 0 ds ε∇, Π s α * N,s -α N,s Π * s HS . (4.50)
Consider the second term in Eq. (4.50). We have:

ω N,s , V * ρ s , ε∇ HS = ε ω N,s , ∇(V * ρ s ) HS ≤ ε σ dp | V (p)||ρ s (p, σ)||p| e ip•x , ω N,s HS ≤ ε dq | V (q)|(1 + |q|) 2 sup p 1 1 + |p| e ip•x , ω N,s HS , (4.51)
where we used Eq. (4.44) and

|ρ s (p, σ)| ≤ 1 N σ dx ω N,t (x, σ; x, σ) = 1 . (4.52)
Consider now the third term in Eq. (4.50). We have, thanks to the representation (4.44):

ω N,s , X s , ε∇ HS ≤ 1 N dq | V (q)| ω N,s , e iq•x ω N,s e -iq•x , ε∇ HS ≤ 2 N dq | V (q)| e iq•x ω N,s , ε∇ e -iq•x HS ≤ 2 N dq | V (q)| ε∇, ω N,s HS , (4.53)
where we used that e -iqx ∇e iq•x = ∇ + iq in the second line. Next, consider the last term in Eq. (4.50). In order to estimate it, we shall use the following identity:

ε∇, e iq•x α N,s e -iq•x α * N,s -ε∇, α N,s e iq•x α * N,s e -iq•x = e iq•x α N,s e -iq•x ε∇, α * N,s + ε∇, e iq•x α N,s e -iq•x α * N,s -α N,s ε∇, e iq•x α * N,s e -iq•x -ε∇, α N,s e iq•x α * N,s e -iq•x = e iq•x α N,s e -iq•x ε∇, α * N,s + e iq•x ε∇, α N,s e -iq•x α * N,s -α N,s e iq•x ε∇, α * N,s e -iq•x -ε∇, α N,s e iq•x α * N,s e -iq•x . (4.54)
Plugging this identity into the expression for the last term in Eq. ( 4.50), we get:

ε∇, Π s α * N,s -α N,s Π * s HS ≤ 4 N dq | V (q)| ε∇, α N,s HS . (4.55)
Using the bounds (4.51), (4.53), (4.55) to estimate the right-hand side of (4.50), we find:

ε∇, ω N,t HS ≤ ε∇, ω N HS + C t 0 ds sup p 1 1 + |p| e ip•x , ω N,s HS + 1 N ε ε∇, ω N,s HS + 1 N ε ε∇, α N,s HS . (4.56)
In order to close the Gronwall argument, we need to estimate the time derivative of the commutator of α N,s with the gradient. We compute:

iε d dt ε∇, α N,t = ε∇, h HF (t)α N,t + α N,t h HF (t) + Π t (1 -ω N,t ) -ω N,t Π t = h HF (t) ε∇, α N,t + ε∇, h HF (t) α N,t + α N,t ε∇, h HF (t) + ε∇, α N,t h HF (t) + Π t ε∇, 1 -ω N,t + ε∇, Π t (1 -ω N,t ) -ω N,t ε∇, Π t -ε∇, ω N,t Π t . (4.57)
The first and the fourth term can be reabsorbed thanks to a unitary conjugation of the commutator with U 3 (t; 0), recall the definition (4.48). Notice that this time the relative sign between these terms is plus, and that the fourth term comes with a complex conjugate. For this reason, we shall consider U * 3 (t; 0) ε∇, α N,t U 3 (t; 0), with U 3 (t; 0) := CU 3 (t; 0)C and C the complex conjugation operator. We then have:

iε d dt U * 3 (t; 0) ε∇, α N,t U 3 (t; 0) = U * 3 (t; 0) -h HF (t) ε∇, α N,t -ε∇, α N,t h HF (t) + iε d dt ε∇, α N,t U 3 (t; 0) , (4.58)
and the first two terms cancel the first and the fourth term in Eq. (4.57). Thus, integrating in time and using the invariance of the trace norm under unitary conjugation, we get:

ε∇, α N,t HS ≤ ε∇, α N HS + 1 ε t 0 ds ε∇, h HF (s) α N,s HS + 1 ε t 0 ds α N,s ε∇, h HF (s) HS + 2 N ε dq | V (q)|
(ε∇ρ s * V )α N,s HS ≤ ∥(ερ s * ∇V )∥ ∞ α N,s HS ≤ ε dp |ρ s (p)|| V (p)| α N,s HS ≤ Cε α N,s HS , (4.61)
while the second term is bounded as:

ε∇, X s α N,s HS ≤ 1 N dq | V (q)| ε∇, e iq•x ω N,s e -iq•x α N,s HS ≤ 1 N dq | V (q)| e iq•x ε∇, ω N,s e -iq•x HS ≤ 1 N dq | V (q)| ε∇, ω N,s HS , (4.62)
where in the second step we used that ∥α N,s ∥ ≤

∥Π s ∥ HS = 1 N dp V (p)e ip•x α N,s e -ip•x HS ≤ 1 N dp | V (p)|∥α N,s ∥ HS ≤ C N ∥α N,

The auxiliary dynamics

In this section we shall introduce a new dynamics, called the auxiliary dynamics and denoted by U N (t; s), on which we will be able to control powers of the number operator via a Gronwalltype inequality. Later, we will prove that U N (t; s) and U N (t; s) are close, in a strong sense. In order to define the auxiliary dynamics, let us denote by L N (t) the generator of the fluctuation dynamics U N (t; s):

iε∂ t U N (t; s) = L N (t)U N (t; s) , U N (s; s) = 1 . (4.73)
The generator L N (t) could be explicitly computed (proceeding as in [START_REF] Benedikter | Mean-field Evolution of Fermionic Mixed States[END_REF]), but this will not be needed in the present setting. Let us define the self-adjoint operator:

L N (t) := L N (t) - 2 N
Re dxdy V (x -y) a * (u t,x )a(v t,y )a(u t,y )a(u t,x ) + a * (u t,y )a * (v t,y )a * (v t,x )a(v t,x ) .

(4.74)

We define the auxiliary dynamics as the solution of the following differential equation:

iε∂ t U N (t; s) = L N (t) U N (t; s) , U N (s; s) = 1 . (4.75)
Comparing with Eq. ( 4.2), it is not difficult to see that:

iε d dt U N (t; 0)ξ, (N + 1) k U N (t; 0))ξ = - 4i N k j=1
Im dxdy V (x -y) × U N (t; 0)ξ, (N + 1) j-1 a(v t,x )a(v t,y )a(u t,y )a(u t,x )(N + 1) k-j U N (t; 0)ξ .

(4.76)

The next proposition allows to bound all powers of the number operator on the auxiliary dynamics. • U N (t; 0)ξ, (N + 1) j-1 a(v t,x )a(v t,y )a(u t,y )a(u t,x )(N + 1) k-j U N (t; 0)ξ .

(4.78)

It is convenient to introduce the notation:

B t,p := dy e -ip•y a(v t,y )a(u t,y ) = dr 1 dr 2 a r 1 a r 2 v t e -ip•x u * t (r 1 ; r 2 ) .
(4.79)

The operator v t e -ip•x u * t can be conveniently expressed in terms of ω N,t and α N,t via the following key rewriting:

v t e -ip•x u * t = v t [ω t , e -ip•x ]u * t + u t α * t e -ip•x u * t -v t e -ip•x α t v T t ; (4.80)
this identity can be checked opening the commutator in the right-hand side, and using the properties (4.12), (4.13). In particular, using that ∥u t ∥ ≤ 1, ∥v t ∥ ≤ 1, we obtain the following estimate:

v t e -ip•x u * t HS ≤ [ω t , e -ip•x ] HS + 2∥α t ∥ HS ≤ Ce c|t| N 1/3 (1 + |p|) , (4.81)
where the last step follows from Proposition 4.3, combined with Assumption 3.1. Thus, by Lemma 2.1 we have:

∥B t,p ξ∥ ≤ Ce c|t| N 1/3 (1 + |p|)∥N 1/2 ξ∥ ∥B * t,p ξ∥ ≤ Ce c|t| N 1/3 (1 + |p|)∥(N + 1) 1/2 ξ∥ ∀ξ ∈ F . (4.82)
Let us estimate I k,j in Eq. (4.78). Inserting 1 = (N + 5) k/2+1-j (N + 5) -k/2-1+j , and using the Cauchy-Schwarz inequality together with the fact that [B t,p , N ] = 2B t,p , we get:

|I k,j | ≤ 1 N dp | V (p)| • (N + 5) k/2+1-j (N + 1) j-1 U N (t; 0)ξ (N + 5) -k/2-1+j B t,p B t,-p (N + 1) k-j U N (t; 0)ξ ≤ C N dp | V (p)| (N + 1) k/2 U N (t; 0)ξ B t,p B t,-p (N + 1) k/2-1 U N (t; 0)ξ . (4.83)
Using twice the estimate (4.82), we have:

B t,p B t,-p (N + 1) k/2-1 U N (t; 0)ξ ≤ Ce c|t| N ε(1 + |p|) 2 (N + 1) k/2 U N (t; 0)ξ , (4.84)
where in the last step we used the bound (4.81). Thus, we obtained: The final claim, Eq. (4.77), follows from the application of Gronwall's lemma.

|I k,j | ≤ Ce c|t| ε dp | V (p)|(1 + |p|) 2 ⟨ U N (t; 0)ξ, (N + 1) k U N (t; 0)ξ⟩ ≤ Ke c|t| ε⟨ U N (t; 0)ξ, (N + 1) k U N (t; 0)ξ⟩ .
We now prove that the fluctuation dynamics and the auxiliary dynamics are close in norm.

Proposition 4.6 (Norm approximation for the fluctuation dynamics.). Under the same assumptions of Theorem 3.3, the following holds. For all ξ ∈ F and for all times t ∈ R:

∥(U N (t; 0) -U N (t; 0))ξ∥ ≤ CN -1/3 exp(4 exp(c|t|)) (N + 1) 3/2 ξ . (4.87)
Proof. The argument is similar to the proof of Proposition 4.5 in [START_REF] Benedikter | Mean-field Evolution of Fermionic Mixed States[END_REF]. By unitarity:

∥(U N (t; 0) -U N (t; 0))ξ∥ = ∥(1 -U N (t; 0) * U N (t; 0))ξ∥ , (4.88)
where:

1 -U N (t; 0) * U N (t; 0) = - t 0 ds ∂ s U N (s; 0) * U N (s; 0) = i ε t 0 ds U N (s; 0) * (-L N (s) + L N (s)) U N (s; 0) = - 2i εN t 0 ds U N (s; 0) * Re dxdy V (x -y) a * (u s,x )a(v s,y )a(u s,y )a(u s,x ) + a * (u s,y )a * (v s,y )a * (v s,x )a(v s,x ) U N (s; 0) , (4.89)
where the last step follows from Eq. (4.74). Plugging this identity in (4.88) we obtain:

∥(U N (t; 0) -U N (t; 0))ξ∥ ≤ 2 εN t 0
ds Re dxdy V (x -y) a * (u s,x )a(v s,y )a(u s,y )a(u s,x ) + a * (u s,y )a * (v s,y )a * (v s,x )a(v s,x ) U N (s; 0)ξ .

(4.90)

Consider the contribution due to the operator a * (u s,x )a(v s,y )a(u s,y )a(u s,x ); the other term is estimated in a similar way. We have:

I := 1 εN t 0 ds dxdy V (x -y) a * (u s,x )a(v s,y )a(u s,y )a(u s,x ) U N (s; 0)ξ ≤ 1 εN t 0 ds dp | V (p)| dΓ(u s e -ip•x u * s )B s,p U N (s; 0)ξ , (4.91)
recall the notation (4.79). By Lemma 2.1, using that ∥u s e -ip•x u * s ∥ ≤ 1, we obtain:

I ≤ 1 εN t 0 ds dp | V (p)| N B s,p U N (s; 0)ξ ≤ 1 εN t 0 ds dp | V (p)| B s,p N U N (s; 0)ξ ≤ Ce c|t| N 1/3 t 0 ds ∥(N + 1) 3/2 U N (s; 0)ξ∥ , (4.92)
where the last step follows from the bound (4.82). Thus, using Proposition 4.5 we find: The other term contributing to the right-hand side of (4.90) can be estimated in the same way. This concludes the proof of Proposition 4.6.

I ≤ Ke c|t| N 1/3

Bound on the growth of fluctuations

We shall now use Proposition 4.5, 4.6, to estimate the growth of fluctuations on the fluctuation dynamics U N (t; s). The next is the main result of Section 4. Proof. Without loss of generality we suppose that k > 0, otherwise the statement is trivial. The proof follows the argument of Section 4.3 of [START_REF] Benedikter | Mean-field Evolution of Fermionic Mixed States[END_REF]. We claim that for every k, n ∈ N there exists c(k) > 0 such that, for all ξ ∈ F:

⟨U N (t; s)ξ, (N + 1) k U N (t; s)ξ⟩ ≤ N (k-(2/3)n) + e 2(n+1
U N (t; s)ξ⟩ ≤ C k ⟨R s ξ, (N + N ) k R s ξ⟩ ≤ C 2k ⟨ξ, (N + 2N ) k ξ⟩ (4.97)
where the first inequality follows from (4.96) and from the fact that the Hamiltonian commutes with the number operator, while the second follows from R * s N R s ≤ C(N + N ), which can be proved similarly to (4.96). This establishes (4.95) for n = 0.

Suppose now that (4.95) holds for n ≥ 0, and let us prove it for n replaced by n + 1, assuming that k -(2/3)(n + 1) ≥ -1/3. We start by writing: Consider now the term II. We have:

⟨U N (t; s)ξ, (N + 1) k U N (t; s)ξ⟩ = ⟨ U N (t; s)ξ, (N + 1) k U N (t; s)ξ⟩ + 2Re ⟨(U N (t; s) -U N (t; s))ξ, (N + 1) k U N (t; s)ξ⟩ + ⟨(U N (t; s) -U N (t; s))ξ, (N + 1) k (U N (t; s) -U N (t; s))ξ⟩ ≡ I + II + III .
|II| ≤ 2 U N (t; s) -U N (t; s) ξ (N + 1) k U N (t; s)ξ ≤ CN -1 3 exp((4 + k) exp(c|t|)) (N + 1) 3 2 ξ (N + 1) k ξ ≤ CN -1 3 exp((4 + k) exp(c|t|)) (N + 1) k+ 3 2 (n+1) ξ 2 (4.100)
where the second step follows from Propositions 4.6, 4.5. Finally, consider the term III. To control it, we use that, setting ∆L N (t) = L N (t) -L N (t):

U N (t; s) -U N (t; s) = 1 iε t s ds ′ U N (t; s ′ )∆L N (s ′ ) U N (s ′ ; s) . ( 4

.101)

Plugging this identity in the expression for III, we have:

III = 1 ε 2 t s ds ′ t s ds ′′ ⟨U N (t; s ′ )∆L N (s ′ ) U N (s ′ ; s)ξ, (N + 1) k U N (t; s ′′ )∆L N (s ′′ ) U N (s ′′ ; s)ξ⟩ , (4
.102) which we estimate as, by Cauchy-Schwarz inequality:

|III| ≤ 2|t -s| ε 2 t s ds ′ ⟨U N (t; s ′ )∆L N (s ′ ) U N (s ′ ; s)ξ, (N + 1) k U N (t; s ′ )∆L N (s ′ ) U N (s ′ ; s)ξ⟩ . (4.103)
Now, by the inductive assumption (4.95), we have:

⟨U N (t; s ′ )∆L N (s ′ ) U N (s ′ ; s)ξ, (N + 1) k U N (t; s ′ )∆L N (s ′ ) U N (s ′ ; s)ξ⟩ ≤ N (k-(2/3)n) + e 2(n+1)c|t| exp(2(n + 1)c(k) exp(c|t -s ′ |)) × ⟨∆L N (s ′ ) U N (s ′ ; s)ξ, (N + 1) 2k+3n ∆L N (s ′ ) U N (s ′ ; s)ξ⟩ . (4.104)
In order to proceed, we rewrite the difference of the generators as ∆L N (s

′ ) = A N (s ′ )+A * N (s ′ ), with A N (s ′ ) = 1 N dxdy V (x -y) a * (u s ′ ,x )a(v s ′ ,y )a(u s ′ ,y )a(u s ′ ,x ) + a * (v s ′ ,x )a(v s ′ ,x )a(v s ′ ,y )a(u s ′ ,y ) . ( 4

.105)

By the Cauchy-Schwarz inequality:

⟨∆L N (s ′ ) U N (s ′ ; s)ξ, (N + 1) 2k+3n ∆L N (s ′ ) U N (s ′ ; s)ξ⟩ ≤ C (N + 1) k+ 3n 2 A N (s ′ ) U N (s ′ ; s)ξ 2 + C (N + 1) k+ 3n 2 A * N (s ′ ) U N (s ′ ; s)ξ 2 ≤ K A N (s ′ )(N + 1) k+ 3n 2 U N (s ′ ; s)ξ 2 + K A * N (s ′ )(N + 1) k+ 3n 2 U N (s ′ ; s)ξ 2 . ( 4

.106)

Both terms can be estimated as in the proof of Proposition 4.6. We find:

⟨∆L N (s ′ ) U N (s ′ ; s)ξ, (N + 1) 2k+3n ∆L N (s ′ ) U N (s ′ ; s)ξ⟩ ≤ Ce 2c|s ′ | N 4/3 (N + 1) k+ 3(n+1) 2 U N (s ′ ; s)ξ 2 ≤ Ce 2c|s ′ | N 4/3 exp(c(k) exp(c|s ′ -s|)) (N + 1) k+ 3(n+1) 2 ξ 2 , (4.107)
for a k-dependent constant c(k) (recall the condition k -(2/3)(n + 1) ≥ -1/3). The last step follows from the control of the growth of powers of the number operator on the auxiliary dynamics, Proposition 4.5. Plugging this estimate in (4.104), and recalling (4.103), we obtain:

|III| ≤ C|t -s|N (k-(2/3)n) + N -2/3 e 2(n+2)c|t| exp((2(n + 1) + 1)c(k) exp(c|t -s|)) × (N + 1) k+ 3(n+1) 2 ξ 2 ≤ C|t -s|N (k-(2/3)(n+1)) + e 2(n+2

Proof of the main result

We are now ready to prove our main result, Theorem 3.3. We start by writing, by unitarity of R t , for j ≤ k:

⟨ψ t , a ♯ 1 x 1 • • • a ♯ 2j x 2j ψ t ⟩ = ⟨R * t ψ t , R * t a ♯ 1 x 1 R t • • • R * t a ♯ 2j
x 2j R t R * t ψ t ⟩ = ⟨U N (t; 0)ξ, (a ♯ 1 (u t,x 1 ) + a ♯ 1 (v t,x 1 ) * ) • • • (a ♯ 2j (u t,x 2j ) + a ♯ 2j (v t,x 2j ) * )U N (t; 0)ξ⟩ .

(5.1)

We then expand the products, and we put the result into normal order, using the canonical anticommutation relations. We call "Wick contraction" any nontrivial anticommutator appearing in this process. In doing this, we produce a constant term, given by:

⟨Ω, R * t a ♯ 1 x 1 R t • • • R * t a ♯ 2j x 2j R t Ω⟩ (5.2)
which is equal to the first term in the right-hand side of (3.6). We are left with estimating the remainder, E ♯ 1 ,...,♯ 2j 2j;t (x 1 , . . . , x 2j ), which is associated to a sum of normal-ordered monomials in the fermionic creation and annihilation operators, each multiplying an integral kernel produced by ℓ contractions, with 0 ≤ ℓ ≤ j -1. The general term contributing to the remainder is bounded as: We observe that 0 ≤ ω N ≤ 1, tr ω N = N , α T N = -α N . Furthermore, it is not difficult to check that:

ω N α N α * N 1 -ω N 2 = ω N α N α * N 1 -ω N . (A.16)
Thus, (ω N , α N ) parametrizes a pure, quasi-free state, with nonzero pairing. Let us estimate the HFB energy of such state. From (A.8), we have:

E HFB N (ω N , α N ) ≤ σ k ε 2 |k| 2 ωN (k, σ) + N V (0) 2 + C∥V ∥ ∞ , (A.17)
where the last term is an upper bound for the exchange energy and for the pairing energy. Next, we estimate:

σ k ε 2 |k| 2 ωN (k, σ) = σ k ε 2 |k| 2 ωFFG N (k, σ) + σ k ε 2 |k| 2 (ω N (k, σ) -ωFFG N (k, σ)) = σ k ε 2 |k| 2 ωFFG N (k, σ) + σ k (ε 2 |k| 2 -µ)(ω N (k, σ) -ωFFG N (k, σ)) , (A.18)
where in the last step we used that ω N and ω FFG N have the same number of particles with given spin. The first term in the right-hand side of (A.18) is the kinetic energy of the free Fermi gas. The second term is bounded above as: Thus, recalling that the direct energy is the same for all translation-invariant states with equal number of particles with a given spin, Eqs. (A.17), (A. [START_REF] Dereziński | Mathematics of quantization and quantum fields[END_REF]) imply the bound

E HFB N (ω N , α N ) ≤ E HFB N (ω FFG N , 0) + KN 1/3 + 2∥V ∥ ∞ , (A.20)
where the factor 2 takes into account the bound for the exchange energy of the free Fermi gas. This concludes the check of the assumptions of Proposition A.1.

Proof. Let us start by proving the first of (A.11). By translation invariance and boundedness of the potential:

E HFB N (ω N , α N ) ≥ k,σ ε 2 |k| 2 ωN (k, σ) + N V (0) 2 -C∥V ∥ ∞ . (A.21)
The last term is a lower bound of the last line in (A.8). Consider the kinetic energy. We rewrite it as:

Nω

  (x 1 , . . . , x k ; y 1 , . . . , y k ) N (x i ; y π(k) ) ;(1.10)

Proposition 4 . 1 (

 41 Growth of fluctuations). Under the same assumptions of Theorem 3.3, the following is true. For any ξ ∈ F, ∥ξ∥ = 1, and k ∈ N,

Proposition 4 . 3 (

 43 Propagation of the semiclassical structure). Under the same assumptions of Theorem 3.3, the following is true. Suppose that ω N and α N satisfy Assumption 3.1.

  N,s HS + ε∇, α N,s HS . (4.59) Consider the second term. We have: ε∇, h HF (s) α N,s HS ≤ (ε∇ρ s * V )α N,s HS + ε∇, X s α N,s HS ; (4.60) the first term is bounded as, setting ρs (p) = σ ρs (p, σ):

Proposition 4 . 5 (

 45 Bound on the growth of fluctuations on the auxiliary dynamics). Under the same assumptions of Theorem 3.3, the following is true. There exists constants C, c > 0 independent of N such that, for all k ∈ N, t ∈ R, ξ ∈ F: ⟨ U N (t; 0)ξ, (N + 1) k U N (t; 0)ξ⟩ ≤ C(exp(k exp(c|t|)))⟨ξ, (N + 1) k ξ⟩ . (4.77) Proof. Let us define: I k,j := 1 N dxdy V (x -y)

( 4 .

 4 85) In conclusion, thanks to the identity (4.76): iε d dt U N (t; 0)ξ, (N + 1) k U N (t; 0)ξ ≤ Cke c|t| ε⟨ U N (t; 0)ξ, (N + 1) k U N (t; 0)ξ⟩ . (4.86)

t 0 ds

 0 exp(3 exp(c|s|))∥(N + 1) 3/2 ξ∥ ≤ K N 1/3 exp(4 exp(c|t|))∥(N + 1) 3/2 ξ∥ . (4.93)

Proposition 4 . 7 (

 47 Bounds for the growth of fluctuations). Under the same assumptions of Theorem 3.3 the following holds true. For all k ∈ N there exists d(k) > 0 such that, for all ξ ∈ F, t ∈ R :⟨U N (t; 0)ξ, (N + 1) k U N (t; 0)ξ⟩ ≤ exp(d(k) exp(c|t|))⟨ξ, (N + 1) α(k) ξ⟩ ,(4.94)with α(k) = (13/2)k for k even and α(k) = (13/2)k + 3/2 for k odd.

( 4 .

 4 98)The term I is estimated using Proposition 4.5:|I| ≤ C(exp(k exp(c|t|)))⟨ξ, (N + 1) k ξ⟩ .(4.99)

2 ξ 2 .( 4 2 ξ 2 ,( 4

 224224 )c|t| exp((2(n + 1) + 1)c(k) exp(c|t -s|))× (N + 1) k+ 3(n+1) .108)Putting together (4.99), (4.100), (4.108) we obtain, possibly for a larger c(k):⟨U N (t; s)ξ, (N + 1) k U N (t; s)ξ⟩ ≤ N (k-(2/3)(n+1)) + e 2(n+2)c|t| exp(2(n + 2)c(k) exp(c|t -s|))× (N + 1) k+3(n+1) .109) which reproduces (4.95) with n replaced by n + 1. This concludes the proof of Proposition 4.7.

1 ) 2 × 1 xr 1 ) 2 =( 5 . 6 ) 2 ≤ 3 λ 3 λ

 1211256233 n∈A |α N,t (x n ; x π(n) )| r∈B |ω N,t (x r ; x π(r) )| × ⟨U N (t; 0)ξ, a ♯r 1 (η r 1 xr • • • a ♯r 2j-2ℓ (η r 2j-2ℓ xr 2j-2ℓ )U N (t; 0)ξ⟩ ,(5.3)where: A, B are disjoint subsets of (1, 2, . . . , 2j) such that |A| + |B| = ℓ; π is a permutation of (1, . . . , 2j) such thatA ∩ π(A) = B ∩ π(B) = A ∩ π(B) = B ∩ π(A) = ∅ ;(5.4) the labels r 1 , . . . , r 2j-2ℓ are the labels of the operators that are not touched by the pairings, that is:(A ∪ π(A) ∪ B ∪ π(B)) ∪ (r 1 , . . . , r 2j-2ℓ ) = (1, 2, . . . , 2j) (A ∪ π(A) ∪ B ∪ π(B)) ∩ (r 1 , . . . , r 2j-2ℓ ) = ∅ ; (5.5)the operators η r i are such that ∥η r i ∥ ≤ 1 (they are either u t or v t ). Let us now estimate the Hilbert-Schmidt norm of the generic contribution. We have:dx 1 . . . dx 2j n∈A |α N,t (x n ; x π(n) )| 2 r∈B |ω N,t (x r ; x π(r) )| ⟨U N (t; 0)ξ, a ♯r 1 (η r • • • a ♯r 2j-2ℓ (η r 2j-2ℓ xr 2j-2ℓ )U N (t; 0)ξ⟩ ∥α N,t ∥ 2|A| HS ∥ω N,t ∥ 2|B| HS dz 1 . . . dz 2j-2ℓ ⟨U N (t; 0)ξ, a ♯r 1 (η r 1 z 1 ) • • • a ♯r 2j-2ℓ (η r 2j-2ℓ z 2j-2ℓ )U N (t; 0)ξ⟩ 2 .Let us denote by m the number of creation operators in the scalar product, 0 ≤ m ≤ 2j -2ℓ.Being the monomial normal ordered, we easily get:dz 1 . . . dz 2j-2ℓ ⟨U N (t; 0)ξ, a ♯r 1 (η r 1 z 1 ) • • • a ♯r 2j-2ℓ (η r 2j-2ℓ z 2j-2ℓ )U N (t; 0)ξ⟩ ∥N m/2 U N (t; 0)ξ∥ 2 ∥N j-ℓ-m/2 U N (t; 0)ξ∥ 2 . (k)|f k ⟩⟨f k | , α N,σ-σ = σ k∈Z (k)(1 -λ(k))|f k ⟩⟨f k | . (A.15)

(ε 2

 2 |k| 2 -µ)(ω N (k, σ) -ωFFG N (k, σ)) = ε 2 σ k (|k| + k F )(|k| -k F )(ω N (k, σ) -ωFFG N (k, σ)) = ε 2 σ k:||k|-k F |≤C (|k| + k F )(|k| -k F )(ω N (k, σ) -ωFFG N (k, σ)) ≤ Kε k:||k|-k F |≤C 1 ≤ KN 1/3 .(A.[START_REF] Dereziński | Mathematics of quantization and quantum fields[END_REF]

k,σ ε 2 F(ε 2 F(ε 2 F|ε 2 |ε 2 |ε 2 ( 1 -F( 1 -F( 1 - 3 σ

 2222221113 |k| 2 ωN (k, σ) = k,σ ε 2 |k| 2 ωFFG N (k, σ) + k,σ ε 2 |k| 2 (ω N (k, σ) -ωFFG N (k, σ)) , (A.22)where ωFFGN (k, σ) = χ(k ∈ B σ F ) with B σF the Fermi ball of an N σ -particle state. Since both states have the same number of particles with given spin,k,σ ε 2 |k| 2 ωN (k, σ) = k,σ ε 2 |k| 2 ωFFG N (k, σ)+ k,σ (ε 2 |k| 2 -µ σ )(ω N (k, σ)-ωFFG N (k, σ)) . (A.23)Consider the last term in the right-hand side of Eq. (A.[START_REF] Fresta | Effective Dynamics of Extended Fermi Gases in the High-Density Regime[END_REF]). We rewrite it as:k,σ (ε 2 |k| 2 -µ σ )(ω N (k, σ) -ωFFG N (k, σ)) = σ,k∈B σ |k| 2 -µ σ )(ω N (k, σ) -1) + σ,k / ∈B σ |k| 2 -µ σ )ω N (k, σ) = σ,k∈B σ F |ε 2 |k| 2 -µ σ |(1 -ωN (k, σ)) + σ,k / ∈B σ |k| 2 -µ σ |ω N (k, σ) . (A.24)We bound below the right-hand side of (A.24) as:k,σ (ε 2 |k| 2 -µ σ )(ω N (k, σ)-ω FFG N (k, σ)) |k| 2 -µ σ |(1 -ωN (k, σ)) |k| 2 -µ σ |ω N (k, σ) ωN (k, σ)) + Cε σ,k / ∈B σ ωN (k, σ)) ≥ σ,k∈B σ ωN (k, σ)) -CN 2/3 = tr ω FFG N (1 -ω N ) -CN 2/σ) -CN 2/3 = tr (1 -ω FFG N )ω N -CN 2/3 , (A.26)we obtain the lower bound: k,σ (ε 2 |k| 2 -µ σ )(ω N (k, σ) -ωFFG N (k, σ)) ≥ Cεtr (1 -ω FFG N )ω N -CN 1/3 . (A.27)

  [START_REF] Benedikter | Optimal upper bound for the correlation energy of a Fermi gas in the mean-field regime[END_REF] is the Hilbert-Schmidt norm squared of the operator E

	♯ 1 ,...,♯ 2j 2j;t	with integral kernel E 2j;t ♯ 1 ,...,♯ 2j	(x 1 , . . . , x 2j

  (v t,x )a(v t,x )a(v t,y )a(u t,y ) + a(v t,x )a(v t,y )α N,t (x, y) -a(v t,x )a(u t,y )ω N,t (x, y)

	+ a(v t,y )a(u t,y )ω N,t (x, x) -h.c..
	(4.16)

y )a(u t,y )a(u t,x ) + a * (u t,x )a * (v t,x )a(v t,y )a(u t,y ) + a * (u t,x )a(v t,y )α N,t (x, y) -a * (u t,x )a(u t,y )ω N,t (x, y) + a(v t,x )a(v t,y )a(u t,y )a(u t,x ) -a *

  .32) HS ≤ ∥ω N ∥ tr , ∥α N,t ∥ HS ≤ e C|t|/(N ε) ∥α N ∥ HS . (4.33) Remark 4.4. The above estimates hold regardless of the size of ∥α N ∥ HS . If ∥α N ∥ HS ≤ CN 1/3 , Eqs. (4.32) allow to propagate the estimates of Assumption 3.1 for all times.

	Furthermore,
	∥ω N,t ∥ 2

  To conclude, let us consider ∥α N,t ∥ HS . We have: N,t U 3 (t; 0) = U * 3 (t; 0) Π t (1 -ω N,t ) -ω N,t Π t U 3 (t; 0) .

	iε	d dt	U * 3 (t; 0)α (4.66)
	Using the invariance of the norm under unitary conjugation:
	∥α N,t ∥ HS ≤ ∥α N ∥ HS + ≤ ∥α N ∥ HS +	1 ε ε 2	0 0	t t	ds Π t (1 -ω N,t ) HS + ω N,t Π t HS ds ∥Π s ∥ HS .	(4.67)
	From the representation (4.44), we get:
									1. Therefore, plugging the bounds (4.61),
	(4.62) in the right-hand side of (4.60):
				ε∇, h HF (s) α N,s HS ≤ Cε α N,s HS +	C N	ε∇, ω N,s HS ,	(4.63)
	and the same estimate holds true for the third term in Eq. (4.59). Thus,
	ε∇, α N,t HS ≤ ε∇, α N HS
				+ C		0	t	ds α N,s HS +	C N ε	ε∇, ω N,s HS +	C N ε	ε∇, α N,s HS .
									(4.64)
	Putting together Eqs. (4.46), (4.56), (4.64), we have:
	sup	1 1 + |p|		e ip•x , ω N,t HS + ε∇, ω N,t HS + ε∇, α N,t HS
		≤ sup	1 1 + |p| t	e ip•x , ω N HS + ε∇, ω N HS + ε∇, α N HS	(4.65)
				+ C	0	ds α N,s HS + ε∇, ω N,s HS
				+ sup	1 1 + |p|	e ip•x , ω N,s HS +	1 N ε	ε∇, α N,s HS .

  s ∥ HS . , ω N HS + ε∇, ω N HS + ε∇, α N HS + |t|e Ct/(N ε) ∥α N ∥ HS , ω N,t HS + ε∇, ω N,t tr + ε∇, α N,t HS , ω N tr + ε∇, ω N HS + ε∇, α N HS + |t|e Ct/(N ε) ∥α N ∥ HS e C|t| ,

							(4.68)
	Plugging this estimate in (4.67), we get:
				∥α N,t ∥ HS ≤ ∥α N ∥ HS +	C N ε	0	t	ds ∥α N,s ∥ HS ;	(4.69)
	by Gronwall's lemma,
				∥α N,t ∥ HS ≤ e Ct/(N ε) ∥α N ∥ HS ,	(4.70)
	which proves the last of (4.33). Coming back to (4.65), we obtained:
	sup p	1 1 + |p|	e ip•x , ω N,t HS + ε∇, ω N,t HS + ε∇, α N,t HS
	≤ sup p e ip•x + C 1 1 + |p| t 0 ds ε∇, ω N,s HS + sup p	1 1 + |p|	e ip•x , ω N,s HS +	1 N ε	ε∇, α N,s HS .
							(4.71)
	By Gronwall's lemma:
	sup p e ip•x ≤ sup 1 1 + |p| p 1 1 + |p|	e ip•x (4.72)
	which implies the estimates (4.32). This concludes the proof of Proposition 4.3.

  + the positive part. The claim (4.94) follows after choosing the smallest n for which(k -(2/3)n) + = 0, that is n = (3/2)k if k is even, or n = (3/2)(k + 1/3) if k is odd.We are left with proving (4.95). The proof is by induction. Without loss of generality, we can suppose that k -(2/3)n ≥ -1/3. Let us check Eq. (4.95) for n = 0. By Eq. (4.3),

	we easily get:
	R t N R * t ≤ C(N + N ) .	(4.96)
	Therefore,
	⟨U N (t; s)ξ, (N + 1) k

)c|t| exp(2(n + 1)c(k) exp(c|t -s|))⟨ξ, (N + 1) 2k+3n ξ⟩ , (4.95) with (•)

Acknowledgements. We thank Benjamin Schlein for useful discussions. S.M. and M.P. acknowledge financial support by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program ERC StG MaMBoQ, n. 802901. S.M. acknowledges financial support from the MSCA project ConNEqtions, n. 101056638. The work of S.M. and M.P. has been carried out under the auspices of the GNFM of INdAM.

The right-hand side of (5.7) can be bounded thanks to Proposition 4.7. We get:

2j exp(d(j) exp(c|t|)) .

(5.8)

Since ∥ω N,t ∥ 2 HS ≤ N and ∥α N,t ∥ 2 HS ≤ CN 2/3 , we obtain:

recall that |A| + |B| = ℓ ≤ j -1. Let p be the absolute value of the difference between the number of a's and a * 's in the left-hand side of (5.1). If p = 0, it is possible to have |A| = 0 and |B| = j -1. Using that the number of terms contributing to E ♯ 1 ,...,♯ 2j 2j;t (x 1 , . . . , x 2h ) is bounded by C j (2j)!, we easily get:

∥ 2 HS ≤ C(j)N j-1 exp(d(j) exp(c|t|))

if p = 0.

(5.10) Suppose now that p ≥ 2. Then, |B| ≤ j -p/2, and |A| ≤ (j -1) -|B|. Summing over the allowed |A| and |B|, it is not difficult to see that:

∥ 2 HS ≤ C(j)N j-2/3-p/6 exp(d(j) exp(c|t|)) if p ≥ 2.

(5.11)

Putting together (5.10) and (5.11), the final claim (3.7) follows. This concludes the proof of Theorem 3.3.

A Check Assumption 3.1 for quasi-free states of the torus

Here we shall check the validity of Assumption 3.1 for translation-invariant fermionic states on the torus T 3 of side 2π. We shall consider states with N σ particles with spin σ ∈ {↑, ↓}, N = N ↑ + N ↓ . Let ω N , α N be the reduced one-particle density matrix and the pairing matrix of a pure quasi-free state. Translation invariance means that:

Therefore, the density ρ ω (x) = ω N (x; x) is spatially constant and equal to N σ /(2π) 3 . We introduce the Fourier symbol associated to ω N (x, σ; y, σ) as:

It follows that 0 ≤ ωN (k, σ) ≤ 1 and that k∈Z 3 ωN (k, σ) = N σ . In terms of this object, the kinetic energy can be written in momentum space as:

A simple example of translation invariant state is the free Fermi gas:

where µ has to be thought as a diagonal matrix in spin space, (µf )(x, σ) = µ σ f (x, σ). In this appendix we shall suppose that the Fermi ball is completely filled:

Given a general quasi-free state (ω N , α N ), its Hartree-Fock-Bogoliubov energy is:

If the state is translation-invariant, the HFB functional can be written as:

Notice that the direct energy is constant, for all translation-invariant states with N particles.

In what follows, besides translation invariance and a given number N σ of particles with spin σ, we shall also assume that the state is pure:

We shall be interested in the properties of quasi-free states which are energetically close to the free Fermi gas.

Proposition A.1 (Bounds for α N ). Let V be bounded. Suppose that (ω N , α N ) is a pure, quasi-free and translation invariant state, and that:

Then, there exists C > 0 such that:

with α N ̸ = 0 satisfying (A.10). For simplicity, suppose that N is even and

and λ(k) smoothly interpolates between 1 and 0 for ||k| -k F | ≤ C. We shall suppose that:

Then, we consider (ω N , α N ) such that:

Coming back to (A.21), we have:

where in the last step we used that the sum of the first two terms in the right-hand side of the first line in (A.28) reproduces the Hartree energy of the free Fermi gas, and we added and subtracted the exchange energy of the free Fermi gas, which is bounded by a constant. Hence, assumption (A.10) and Eq. (A.28) immediately imply the estimate:

To conclude, we observe that:

combined with (A.29), this gives the first estimate in (A.11). Let us now prove the second estimate in (A.11). We have:

where the last step follows from Cauchy-Schwarz inequality for traces. Thus, by Eq. (A.9):

The right-hand side of (A.32) can be bounded proceeding as for the first estimate in (A.11).

Consider the difference of the kinetic energies of ω N and ω FFG N . Proceeding as in (A.24):

which we bound below by dropping the absolute values:

Next, we use that, similarly to (A.30):

-tr

where we used that [ω FFG N , ε∇] = 0. Thus, we obtained:

where the right-hand side of (A. [START_REF] Porta | Mean Field Evolution of Fermions with Coulomb Interaction[END_REF]) is what appears in the right-hand side of (A.34). Hence, plugging this bound in (A.34), we get:

where in the last step we used (A.29). Proceeding as in (A.27), (A.28), the bound (A.37) implies:

combined with assumption (A.10), we obtain:

Recalling (A.32), the second bound in (A.11) follows. This concludes the proof of Proposition A.1.

The next proposition establishes the validity of the commutator estimates for ω N in (3.3), in Hilbert-Schmidt norm, for translation-invariant states energetically close to the free Fermi gas. Notice that for a translation invariant state the second estimate in (3.3) is trivially true. The next proposition proves the first estimate in (3.3).

Proposition A.3 (Bounds for ω N ). Let V be bounded. Suppose that (ω N , α N ) is a pure, quasi-free and translation invariant state, and that:

Then, there exists C > 0 such that:

Proof. We start by observing:

where we used that 0 ≤ ω N ≤ 1, and the cyclicity of the trace. We will use the energetic constraint (A.40) to show that:

which combined with (A.42) gives the final claim (for p = 0, this has been obtained in the proof of Proposition A.1). Let us denote by ωN,p (k, σ) the Fourier symbol of (

) and 0 ≤ ωN,p (k) ≤ 1, we have, from Eq. (A.25):

This implies:

(A.45) Therefore, we obtain the following lower bound:

(A.46)

Using that the direct energy of (ω N , α N ) and of (ω FFG N , 0) are equal, and the fact that the exchange energy and the pairing energy are bounded by a constant, Eq. (A.46) gives:

Thus, by the assumption (A.40), we get:

where the last inequality follows from an explicit computation of the trace norm, see e.g. [START_REF] Benedikter | Mean-field evolution of fermionic systems[END_REF]. Combining this estimate with (A.42), we immediately get:

which concludes the proof of (A.41).

Finally, we conclude the section by showing that the only translation invariant ground state of the HFB functional on the torus is the free Fermi gas.

Proposition A.4 (Translation invariant ground state). Let V be bounded. Suppose that (ω N , α N ) is a pure, quasi-free and translation invariant state. Then, for N large enough:

Equality holds if and only if (ω N , α N ) = (ω FFG N , 0). Proof. We introduce the notation:

Then, is it not difficult to see that:

, and where we recall that µ has to be thought as a diagonal matrix in spin state: (µf )(x, σ) = µ σ f (x, σ). We will use the first term in the right-hand side to control the other two. We write:

Since by assumption dist(σ(-∆), k σ2 F ) = 1/2, and using the bound ∥X ω ∥ ≤ C/N together with the fact that h ω is invariant by translation and thus commutes with -∆, we have for N large enough:

ω FFG N = 1(-ε 2 ∆ -µ ≤ 0) ≡ 1(h ω -µ ≤ 0) .