On the effective dynamics
Bose-Fermi mixtures ]

On the effective dynamics
Bose-Fermi mixtures

Esteban Cárdenas
[Esteban Cárdenas]Department of Mathematics,
University of Texas at Austin,
2515 Speedway,
Austin TX, 78712, USA

Joseph K. Miller
[Joseph K. Miller]Department of Mathematics,
University of Texas at Austin,
2515 Speedway,
Austin TX, 78712, USA

Nataša Pavlović
[Nataša Pavlović]Department of Mathematics,
University of Texas at Austin,
2515 Speedway,
Austin TX, 78712, USA

In this work, we describe the dynamics of
a Bose-Einstein condensate interacting with a degenerate Fermi gas, at zero temperature.
we analyze the mean-field approximation
of the
many-body Schrödinger dynamics
and prove
emergence   of     a    coupled Hartree-type system of equations.
We obtain rigorous error control
that yields a non-trivial scaling window in which the approximation
is meaningful.
Second, starting from this Hartree system,
we identify a  novel scaling regime
in which the fermion  distribution
behaves semi-clasically, but the boson field remains quantum-mechanical;
this is  one of the main contributions of the present article.
In this regime, the  bosons are much lighter and more numerous than the fermions.
We then prove convergence to a coupled Vlasov-Hartee system of equations
with an explicit convergence rate.

§ INTRODUCTION

In this work, we study the dynamics of a gas composed of $M$ identical fermions and $N$ identical bosons
moving in   Euclidean space $\R^d$,
for spatial dimensions $d\geq2.$
The Hilbert space for the system is the tensor product

\begin{equation}
\mathscr{H}
\equiv
L^2_a( \R^{dM}) \otimes L^2_s( \R^{dN}) \  ,
\end{equation}

where  $L_a^2$ and $L_s^2$ correspond to the subspaces of antisymmetric and symmetric functions, respectively;
we neglect any internal degrees of freedom the particles may have.
We assume  that the two systems are non-relativistic,
and interact by means of a   two-body potential $V : \R^d \rightarrow \R$.
Thus, we  study the   Hamiltonian in
dimensionless variables

\begin{equation}
\label{hamiltonian HNM}
\equiv
\frac{\hbar^2}{2 m_F} \sum_{i =1}^M  (- \Delta_{x_i}) \otimes \1
\frac{\hbar^2}{2 m_B} \sum_{j =1 }^N \1 \otimes   (- \Delta_{y_j})
\lambda
\sum_{i,j =1 }^{N,M}
V(x_i - y_j) \    .
\end{equation}

Here,  all the physical parameters are effective (dimensionless) quantities,
and for simplicity we shall make no distinction in the notation.
$\hbar>0$  corresponds to  Planck's   constant,
$\lambda>0$ is the strength   of the interaction, and $m_F$ and $m_B$ are the respective masses of fermions and bosons.
The first two terms in (<ref>)
correspond to the kinetic energies of the two subsystems,
whereas the third one corresponds to the potential energy associated to their interaction.
In particular,
$V$ will always be assumed to be regular enough so that the Hamiltonian
$H$ is self-adjoint in its natural domain, and
its time evolution is well-defined.

Since Bose-Fermi mixtures   have not been extensively studied in the mathematical literature,
to better focus on the interplay between bosons and fermions through the  potential $ V $,
this work intentionally omits two-body boson-boson and fermion-fermion interactions:

\begin{equation}
\qquad
\mathbb W _{BB}
\mu_{BB}
\sum_{j_1 < j_2}
W(y_{j_1}  - y_{j_2} )
\qquad
\t{and}
\qquad
\mathbb W _{FF} =
\mu_{FF}
\sum_{i_1 < i_2}
W (x_{i_1}- x_{i_2})
\end{equation}

in the Hamiltonian (<ref>).
On the other hand, systems of interacting bosons, and interacting fermions have been broadly studied in the literature in the last few decades and we now give a brief historical review of the corresponding mathematical results.

§.§ Historical background on single species systems

The main goal of this paper is to study the   dynamics generated by the Hamiltonian
As is often the case in many-body systems,   the associated   Schrödinger equation

\begin{eqnarray} \label{schrodinger intro}
i \hbar \partial_t \psi(t) = H \psi(t) \ , \qquad \psi (t) \in \mathscr H
\end{eqnarray}

is impossible to solve exactly
and one must rely on effective approximations
to understand the  physical   behaviour of the system.
In this context, one of the best understood approximations is the mean-field approximation.
In broad terms,
one assumes here that the two-body interaction potential    between particles is weak,
but the number of particles is large–one then   replaces
the total   interaction  by its average over position densities.
Let us briefly describe what this idea leads to   in the context of single species gases.

For cold systems of interacting bosons,
assuming that the  potential $W(y_i - y_j)$
that mediates their interaction
varies over length scales
comparable to the size of the gas,
the mean-field approximation leads
to the derivation of the Hartree equation

\begin{eqnarray}
\label{hartree}
&  i \partial_t \vp(t) = -  \frac{1}{2}
\Delta \vp(t) +
\big(W* |\vp (t)|^2 \big) \,  \vp(t)  \ ,
\end{eqnarray}

where the solution $ \vp(t) \in L^2(\R^d)$  corresponds to the wave-function
of a Bose-Einstein condensate;
if $ W(x)$ varies on much smaller scales, one obtains
the Gross-Pitaevskii or nonlinear Schrödinger (NLS) equation.
The  Hartree equation (<ref>)
was first rigorously derived in [58] for bounded potentials,
and later in [27]
for Coulomb systems.
More recently,
the second quantization formalism
was employed  in [55, 16]–inspired
by studies on the fluctuation dynamics originated in
[43, 34, 35]–to obtain a quantitative convergence rate;
see also [19] where uniform-in-time bounds for error estimates are proven.
We also refer the reader to the following
non-exhaustive list of references
on related works [32, 2, 45, 53, 40, 41, 42]
on the derivation of the Hartree equation,
and to
[23, 26, 24, 25, 1, 44, 14] on the derivation of the Gross-Pitaevskii equation.

For  cold gases of $M$ interacting fermions,  one obtains the Hartree-Fock equation

\begin{equation}
\label{hartree fock}
i \hbar  \partial_t \omega(t)
\Big[
-\frac{\hbar^2}{2} \Delta + W*\rho (t) - X(t)  , \omega(t)
\Big] \ ,
\end{equation}

$ \rho(t,x)= M^{-1} \omega(t;x,x)$
is the density of particles,
and $X(t) $ is the so-called
exchange term.
Here, the solution
$\omega(t) $
is a positive,  trace-class operator on $ L^2(\R^d)$
whose trace is equal to $M$;
it ought to describe an interacting    Fermi gas of $M$ particles.
The Hartree-Fock equation
has been historically studied in two scaling regimes.
The first derivation in the “microscopic regime" (namely, physical scales for which $\hbar=1$)
was carried out in
[4] for regular interactions,
and later improved in
[31]
for Coulomb systems.
In the “macroscopic regime"
(namely, physical scales for which
$ \hbar = M^{-1/d}$),
the first derivation was carried out in [22] for real analytic potentials, yielding an optimal convergence rate for short macroscopic times.
More recently, the derivation was revisited in
[9]
using second quantization methods,
relaxing significantly the regularity of the potentials,
and extending the time-validity of the derivation–as a tradeoff,
one here requires additional semi-classical structure on the initial data.
This inspired substantial work in the literature;
see for instance
[10, 7, 54, 30].

On the other hand, the $\hbar\downarrow 0 $ limit of the Hartree-Fock equation (<ref>)
leads to the Vlasov equation

\begin{align}
\label{vlasov}
& (\partial_t + p \cdot \nabla_x + F_f(t) \cdot \nabla_p) f (t,x,p)= 0
\end{align}

where $ F_f (t,x) = - \int \nabla W(x-y ) f(t,y,p) d y dp$
is a mean-field force term,
and $f(t )\in L^1_+ (\R_{x,p}^{2d} )$ is a  macroscopic phase-space distribution function.
In particular, the Pauli exclusion principle, $ 0 \leq \omega \leq 1$ viewed in the sense of quadratic forms, still holds in the macroscopic limit   $ 0 \leq f \leq 1$, in the pointwise sense.
One can therefore understand
the solution of the Vlasov equation (<ref>) as the description
of a macroscopic gas with quantum features.
As for  the  derivation of the Vlasov equation
from interacting quantum systems,
the first works on the subject are
[52, 59].
Here, the derivation is carried out
in the   macroscopic regime,
by studying directly
the BBGKY hierarchy
associated to
the many-body Schrödinger dynamics.
The convergence from the
Hartree/Hartree-Fock equation
to the Vlasov equation
was later analyzed in
[48, 49, 33],
although providing no convergence rate.
The first work
to provide a convergence rate for regular potentials
was [3].
Later, the derivation of a convergence rate
from the Hartree to the Vlasov dynamics was revisited
and established in
[8, 17]
for a larger class of potentials.

§.§ Bose-Fermi mixtures

Investigating degenerate mixtures of bosons and fermions is an extremely active area of research in experimental physics for constructing and understanding novel quantum bound states such as those in superconductors, superfluids, and supersolids [28, 56, 21]. These ultra-cold Bose-Fermi mixtures are fundamentally different from degenerate gases with only bosons or fermions. They not only show an enriched phase diagram, but also a fundamental instability due to energetic considerations coming from the Pauli exclusion principle [51]. In particular, the fermionic particles maintain a higher energy than the bosonic particles in the ground state, causing a physical instability due to the energetic difference. This difference bounds from above the number of fermions allowed to exist in these doubly degenerate mixtures [18]. On the other hand, by varying the ratio of masses of bosons and fermions in these mixtures, experimentalists have studied Bose-Einstein condensates with bose-bose interactions mediated by fermions [20].

Inspired by this activity in the physics community, in this paper we start exploring
the mathematical theory of Bose-Fermi mixtures,
by studying the mean-field dynamics of the Hamiltonian introduced in (<ref>).
Here, one of the main challenges is understanding the physical scales of the system
that allow for suitable analysis.
Indeed, the Pauli Exclusion Principle implies that for  confined gases of fermions at low temperatures,
fermions have a characteristic energy
that varies in
a scale
\hbar^2 M^{\frac{2}{d}}/ m_F,
for bosons
this is only of order
$\hbar^2 / m_B$.
Thus, finding a   scaling regime in which one can capture
the effective dynamics of the system presents a challenge in itself that we have to address.

Let us informally describe the main results of this paper, stated rigorously in Section <ref>. The first result, formulated in Theorem <ref>,
a quantum mean-field approximation
of the many-body Schrödinger dynamics.
Here, we prove that
the one-particle reduced density matrices for the corresponding fermionic and
bosonic subsystems (see (<ref>) for the definition) are
effectively described by a pair of  interacting variables

\begin{equation}
(\omega, \vp )
: \R
\rightarrow
\mathscr{L}^1
( L^2 ( \R^d)   )
\times L^2 ( \R^d  )   \ ,
\end{equation}

satisfying the system of self-consistent equations,
which we  shall refer to
as the  Hartree-Hartree equation.

\begin{align}
\label{hartree hartree intro}
\begin{cases}
& i \hbar \partial_t \omega
[   - (\hbar^2 / 2m_F)\Delta
\lambda N   \,
(V  * \rho_B )  ,
\omega
]			\\
&   i \hbar \partial_t  \vp
- (\hbar^2 / 2m_B)\Delta  \vp
\lambda  M    \,
(V  *  \rho_F)
\vp
\end{cases}   \   .
\end{align}

Here $\rho_F (t , x) = \frac{1}{M } \omega( t; x , x )$
and $\rho_B(t,x) = |\vp(t,x)|^2$
correspond to the
fermionic and bosonic position densities, respectively.
A few comments are in order.

* The derivation of the above equation can be heuristically justified as follows.
We assume that at time $t=0$,
the system
has been externally confined by means of an harmonic trap,
at zero temperature.
The many-body wave function is   then
expected to be of the form
$\psi (0) =     \psi_F(0) \otimes \psi_B(0) \in \mathscr H ,$
where the fermionic component $\psi_F(0)$ is a
Slater determinant,
and the bosonic component $\psi_B(0)$ is a fully factorized state–for a precise meaning, see Assumption <ref>.
If the interactions between particles are weak enough,
then the above structure for the wave function is approximately
valid also for later times.
In particular, a direct calculation shows that plugging the formal ansatz
$\psi(t)  =   \psi_F(t) \otimes \psi_B(t)$–preserving the initial data  structure–in the Schrödinger dynamics (<ref>)
leads to  the Hartree-Hartree equation  (<ref>)
as a self-consistent approximation.

* While the informal justification of the emergence of  (<ref>)
may not present a   challenge,
obtaining a scaling regime in which
the above system
describes the leading order dynamics
is non-trivial.
More specifically, using Second Quantization methods,
we provide rigorous control
of error terms, which
become small only in a particular non-empty parameter window–including
both macroscopic and microscopic scaling regimes–see Theorem <ref> for more details.

In our   second result, stated in Theorem <ref>,
we study a new scaling regime (contained in the aforementioned parameter window)
for the combined system.
In this   scaling regime,
the fermion component $\omega$
is described semi-classically,
but the bosonic component $\vp$
remains quantum-mechanical.
More precisely, the regime that we focus on  is given by

\begin{equation}
\label{scaling regime intro}
\lambda = \frac{1}{N}
\ ,
\qquad \hbar =  \frac{1}{M^{\frac{1}{d}}}  \ ,
\qquad  m_B = \hbar \ ,
\qquad m_F =1 \,
\quad \t{and} \quad
N = M^{	 1 +		\frac{1}{d}	} \ .
\end{equation}

Let us note here that under this scaling regime,
the Hamiltonian that drives the boson field  $\vp$
is proportional to $\hbar$.
In other words,  it has the form
$	\hbar ( -  (1/2) \Delta +   V*\rho_F^\hbar ) 	$.
Thus, it follows that factors of  $\hbar$ cancel out in the second equation of (<ref>),
which enable us to analyze the semi-classical limit.
To this end, we consider $f^\hbar = W^\hbar[\omega^\hbar]$,
the Wigner transform of the fermionic component (see (<ref>) for the definition of the Wigner transform)
and prove that in the $\hbar \downarrow 0$ limit,
there is convergence
$ (f^\hbar , \vp^\hbar) \rightarrow (f, \vp)$,
where the latter variables
satisfy a coupled system of equations.
This system has the following form,
and we shall refer to it  as the Vlasov-Hartree equation

\begin{align}
\label{vlasov hartree intro}
\begin{cases}
& 		(\partial_t  + p \cdot  \nabla_x + F_B(t,x)  \cdot  \nabla_p   ) f  = 0 	\\
&		  i \partial_t \vp
- \frac{1 }{2} \Delta   \vp +   ( V*\rho_F   )   \vp
\end{cases}  \ .
\end{align}

Here,  $F_B(t,x) \equiv
- \int\nabla V(x-y ) | \vp(t,y)|^2 \d y
is a mean-field force that the bosons exert over the fermions, and
$\rho_F 	(t,x) \equiv  \int_{	\R^d	} f (t,x,p ) \d p 	 $
is the fermionic position density.
Our proof of convergence is quantitative,
and implements   for the problem at hand  recently developed techniques of Quantum Optimal Transportation (QOT) [36, 37, 38, 39, 46].

In addition to
recognizing a mean-field scaling regime
that allows us to  rigorously derive the Hartree-Hartree system (<ref>),
one of the main contributions of this article
is the identification of a novel mean-field
semi-classical scaling regime
in which the limiting dynamics of  (<ref>) is non-trivial.
To the authors best knowledge, this regime
had not been identified previously in the literature.
We would like to note that
such a  regime can be related
to    systems of particles which  have been studied (individually) in experiments,
including    Ytterbium-171 for fermions,  and Helium-4 for bosons.
Here, the ratio of the masses is of order $m_B/ m_F \simeq 10^{-2}$;
comparable to the value of the effective Planck's constant $\hbar$
at particular length, time and energy scales.
The reader is referred to  Subsection
for more details
on the physical interpretation of our mathematical results.

§.§ Organization of this paper

In Section <ref>
we formulate our main results in Theorem <ref>
In Section <ref>
we give preliminaries on the Second Quantization formalism, that we will extensively use.
In Section <ref> we study the dynamics of the fluctuations
around a combined Bose-Einstein condensate and degenerate Fermi gas,
which then we use to prove Theorem <ref> in Section <ref>.
Next, in Section <ref>
we adapt  the formalism of   Quantum Optimal Transportation
and utilize it
to prove Theorem <ref>.
Finally, we include  Appendix <ref> where we  state some basic well-posedness results
regarding the PDEs introduced in this paper,
and  Appendix <ref>
where we give details of the calculation of the infinitesimal generator of the fluctuation dynamics.

§.§ Acknowledgements

E.C is very thankful to  F. Golse for an enlightening conversation regarding QOT,
and to N. Benedikter, M. Porta and  C. Saffirio for  helpful discussions
regarding the mean-field dynamics of Fermi systems.
E.C. gratefully acknowledges support from the Provost’s Graduate Excellence Fellowship at The University of Texas at Austin and from the NSF grant DMS-2009549,
and the NSF grant DMS-2009800 through Thomas Chen.
J.M. gratefully acknowledges support from the Provost’s Graduate Excellence Fellowship at The University of Texas at Austin and from the NSF grants No. DMS-1840314.
N.P. gratefully acknowledges support from the NSF under grants No. DMS-1840314, DMS-2009549
and DMS-2052789.

§ MAIN RESULTS

In this section, we describe the main results of this article,
that have already been announced in the introductory section.
In particular, in subsection <ref> we present Theorem <ref>,
the quantum mean-field approximation
of the many-body Schrödinger dynamics.
Here, we prove
an upper bound on the error term
that comes from the approximation of
the one-particle reduced density matrices for the corresponding fermionic and
bosonic subsystems,
and the solution of the Hartree-Hartree equation (<ref>).
In subsection <ref> we present Theorem <ref>, in which
we study the scaling regime (<ref>) for the Bose-Fermi system.
We prove that in the $\hbar \downarrow 0$ limit,
there is convergence towards
the Vlasov-Hartree equation (<ref>).
As stated in the Introduction,
one of the main contributions of this article
is the identification of a semi-classical scaling regime
in which the limiting dynamics of the coupled system is non-trivial–to the authors best knowledge, this regime
had not been identified previously in the literature.
In subsection <ref> we briefly discuss
the strategy of our proofs and the methods that we employ.
Finally, in subsection <ref>
we discuss the physical interpretation of our main mathematical results.

Before we move on to the
main results of this section, let us introduce some notation that we will be using in the rest of the article.

* $L^p (\R^n )$   denotes the   Lebesgue spaces of $p$-th integrable functions,
for $p \in [1,\infty]$.
The subset of non-negative functions is denoted by $L^p_+(\R^n )$.

* $\mathscr P_m (\R^n)$ is the space of probability measures on $\R^n$
that have finite $m\in \mathbb N $ moments.

* $\mathscr S(\R^n ) $ denotes the space of Schwartz functions of rapid decay.

* $W^{k,p} (\R^n )$
for $k\in \mathbb N$ and $p\in[1,\infty]$,
denotes the Sobolev space of functions with derivatives of order $k$, that are $p$-th integrable.

* $H^s (\R^n ) = W^{s,2} (\R^n )$ for $s \geq 1 $
and $\dot H ^s (\R^n) $ is the usual homogeneous Sobolev space.

* $\mathscr L ^1 (X)$ stands for the Banach space of trace-class operators over $X$, endowed
with the  norm  $\|   A \|_{\tr} \equiv \tr |A|$.
Similarly, $\mathscr L ^2 (X)$ is the space of Hilber-Schmidt
operators with norm $\| A\|_{HS} \equiv   \|   A^*A  \|_{\tr}^{1/2}$.

* We say that $C>0$ is a
constant if it is a positive number,
independent of the physical parameters
$N,M , \hbar , \lambda , m_F, m_B$ and $t$.

* $\<\xi \> =  (1 + \xi^2)^{1/2}$ denotes the standard angle bracket.

§.§ The mean-field approximation

As we have previously discussed, the main interest in this article is
to consider the mean-field  dynamics generated by the Hamiltonian $H$,
introduced in
To this end, we introduce the wave function of the system at time $t \in \R$

\begin{eqnarray}
\label{wave function}
\psi(t)
\equiv
\exp \Big(
- i t H / \hbar
\Big)  \psi
\end{eqnarray}

where $\psi \in \mathscr H $ is the initial data of the system.
Since our gas corresponds of two subsystems–each composed of identical particles–it will be crucial
to introduce
the corresponding fermionic and bosonic one-particle reduced density matrices.
are the time-dependent
trace-class operators
$\gamma_F(t), \,  \gamma_B(t) \in \mathscr{L}^1(L^2(\R^d))$
whose kernels are defined as  the partial traces

\begin{align}
\label{definition kernels}
\begin{cases}
&  \gamma_F (t ; x ,x')
\equiv
M \,
\int_{  \R^{d(M-1)} \times \R^{dN }	}
\psi (  t ;  x ,  \bx_{M-1};  \by_N  )
\overline{ \, \psi \,  }   ( t; x' ,   \bx_{M-1};  \by_N )
\d   \bx_{M-1} \d    \by_N  	\\
& 	\gamma_B (t ; y ,y')
\equiv
N  \,
\int_{  \R^{dM } \times \R^{d(N-1)  }	}
\psi   (   t;    {\bx}_M ;   y  ,   {\by}_{N-1}  )
\overline{ \,  \psi \, }  (   t ;   {\bx}_M ;   y'  ,   {\by}_{N-1}  )
\d  {\bx}_M  \d  {\by}_{N-1}    \\
\end{cases}
\end{align}

for $t \in \R$ and $x,x',y, y' \in \R^d$.
we denote by $\bx_{M-1} = (x_1, \ldots, x_{M-1})$,
$\by_N = (y_1, \ldots, y_N)$
and similarly
and $\by_{N-1}$, the variables that are being traced out.
In particular, we note here that
the normalizations are chosen so that
for all times $t\in \R $
there holds

\begin{equation}
\tr\gamma_F(t) = M
\qquad
\t{and}
\qquad
\tr \gamma_B(t) =N  \ .
\end{equation}

We describe now
the conditions that we shall impose
in the initial data $\psi \in \mathscr H $
associated to the solution of the Schrödinger dynamics (<ref>).
Physically, the situation we consider
concerns the description of an initially prepared cold gas of fermions and bosons.
At zero temperature, we
expect the
fermion component to be described as a degenerate Fermi gas–parametrized by a   Slater determinant–whereas the boson gas
undergoes  Bose-Einstein condensation, described by a single-particle wave function. This will be made rigorous in Assumption <ref>.
In addition, we require additional assumptions on the scales
in which the Fermi gas varies–see Remark <ref> for more details.

Let us now discuss the effective dynamics of this system.
Indeed, if the interactions between particles are weak enough, we expect the zero temperature structure described above to approximately persist for times $t>0$.
Thus, we expect to find a mean-field approximation for the reduced density matrices
$\gamma_F$ and $\gamma_B$
in terms of  a pair of  interacting   variables

\begin{equation}
(\omega, \vp )
: \R
\rightarrow
\mathscr{L}^1
( L^2 ( \R^d)   )
\times L^2 ( \R^d  )
\end{equation}

that solve a self-consistent equation.
A formal calculation using a time-dependent Slater determinant/fully factorized
ansatz combined with replacing the full interaction $V(x-y)$
with an average over the position densities then yields

\begin{align}
\label{hartree hartree}
\begin{cases}
& i \hbar \partial_t \omega
[   - (\hbar^2 / 2m_F)\Delta
\lambda N   \,
(V  * \rho_B )  ,
\omega
]			\\
&   i \hbar \partial_t  \vp
- (\hbar^2 / 2m_B)\Delta  \vp
\lambda  M    \,
(V  *  \rho_F)
\vp  	\\
& (\omega ,\vp)(0)  = (\omega_0 , \vp_0) \in
\mathscr{L}^1
( L^2 ( \R^d)   )
\times L^2 ( \R^d  )
\end{cases}   \   ,
\end{align}

up to leading order.
Here $\rho_F (t , x) = \frac{1}{M } \omega( t; x , x )$
and $\rho_B(t,x) = |\vp(t,x)|^2$
correspond to the
fermionic and bosonic position densities, respectively.
We shall refer to (<ref>)
as the  Hartree-Hartree equation.

We are now ready to rigorously state our assumptions which we indicated above.

[Schrödinger initial data]

We assume that  the initial data $\psi\in\mathscr H$
satisfies the following conditions.

* (Zero temperature)
$\psi$ is  a factorized state of the form

\begin{equation}
\psi = \psi_F \otimes \psi_B \ .
\end{equation}

Additionally, each factor satisfies the following assumptions.

There exists a rank-$M$ orthogonal projection $\omega_0 = \sum_{i=1}^M \ket{\phi_i} \bra{\phi_i} $
on the one-particle space $L^2(\R^d)$
such that

\begin{equation}
\psi_F
(x_1 , \ldots, x_M )
\frac{1}{\sqrt{ M! }}
\det_{1 \leq i,j\leq  M }
\big[
\phi_i(x_j)
\big]  \ .
\end{equation}

(1.2) There exists a unit vector  in the one-particle space $\vp_0 \in L^2(\R^d)$,
such that

\begin{equation}
\psi_B (y_1, \ldots, y_N )
\vp_0 (y_1) \cdots \vp_0 (y_N ) \  .
\end{equation}

* (Semi-classical bounds)
Let $(\omega,\vp)$ be the solution of the Hartree-Hartree system
with initial data $(\omega_0,\vp_0 ) \in \mathscr L ^1 (L^2(\R^d)) \times L^2(\R^d)  $.
Then, we assume that there exists a constant $C>0$ such that
for all $t\in \R$
there holds

\begin{eqnarray}
\label{semiclassical structure}
\|     [ e^{ i \xi \cdot x } , \omega (t) ]	   \|_{\tr }
\leq C
\exp(C | t| )
\hbar \<\xi \> \  \ , \qquad \forall \xi \in \R^d  .
\end{eqnarray}

Let us observe that under the above assumptions,
one can calculate that the following relations hold at $t = 0$

\begin{equation}
\gamma_F(0) =
\omega_0
\qquad
\t{and}
\qquad
\gamma_B(0)
\ket{\vp_0}
\bra{\vp_0} \ .
\end{equation}

In other words, the initial data
is such that
the one-particle reduced density matrices are given by
\ket{\vp_0}
\bra{\vp_0}
$, respectively.

Let us now comment on the semi-classical bounds
that we present in Assumption <ref>.
Two comments are in order.

These estimates were first considered by the authors in
[9]
in the derivation of the Hartree-Fock equation
from interacting Fermi systems.
In comparison to the present work,
the scaling is chosen so that $m_F =1 $, $\lambda =1/M$ and
$\hbar = 1/M^{1/d}$.
In this situation, the bounds (<ref>)
follow from the propagation-in-time
along the solutions of (<ref>)
of the following trace estimates
on the initial data

\begin{equation}
\label{semiclassical 1}
\|       [ x , \omega_0 ]        \|_{\tr}
\leq 	 C M \hbar
\quad
\t{and}
\quad
\|       [     i\hbar \nabla  , \omega_0 ]      \|_{\tr}
\leq  C M \hbar
\end{equation}

From the physical point of view,
the commutator estimates
contained in
record macroscopic or semi-classical behaviour.
Namely, in macroscopic units   the value of the effective Planck constant $\hbar$ is small and
so are the quantum-mechanical commutators on the left hand side of (<ref>).
We refer the reader to the original reference for a more in-depth physical discussion.
Finally, let us note
that–to the authors best knowledge–the only     example in $\R^d$
in which the bounds
have been verified  consists
of states describing a non-interacting system of fermions, trapped by an harmonic potential

In this work, we   view   the  bounds (<ref>)
as a condition that we impose  on both the initial data, and the scaling regime under consideration.
In particular, one might ask whether
the bounds (<ref>)
can be guaranteed by properly choosing initial data
by restricting the physical parameters.
We show that this is exactly the case.

Let us further explain.
First of all, one may repeat the argument carried out in <cit.>
without any scaling restrictions
and prove that  the   trace estimates for the same (externally trapped) initial data take the form

\begin{equation}
\label{semiclassical initial data}
\| [ x ,\omega_0]   \|_{\tr}
\leq
m_F^{-1/4 }
M \sqrt{\hbar M^{-1/d} }
\qquad
\t{and}
\qquad
\| [ i\hbar \nabla  ,\omega_0]   \|_{\tr}
\leq
C m_F^{1/4 }
M \sqrt{\hbar M^{-1/d}}
\end{equation}

where $C$ is a constant that only depends on the strength of the trap $V_{\t{ext}}(x) = \kappa x^2/2$.
As a side remark, we note that
with  the scaling mentioned in (i)
recovers (<ref>).
Furthermore, one may adapt the proof of
to show that the following bounds are propagated along the solutions of
the Hartree-Hartree equation  (<ref>)

\begin{align}
\sup_{ \xi \in \R^d}
\<\xi \>^{-1 }
\|         [e^{ i \xi \cdot x } , \omega(t)]   \|_{\tr}
& \leq
\big(
m_F^{ -1/4 }
m_F^{ -3/4}
\big)
e^{\big[     \max\big\{ 1, \lambda N / m_F\big\}       |t| \big] }
M \sqrt{\hbar M^{-1/d}}
\end{align}

where we have kept track of all the physical parameters for comparison.
In particular,
one here obtains (<ref>)
provided  we assume that

\begin{align}
\label{scaling 2}
m_F   \geq 1  , \qquad
\lambda N \leq m_F
\quad
\t{and}
\quad
M^{ -\frac{1}{d}}
\leq \hbar \ .
\end{align}

In other words, by considering externally trapped initial data satisfying (<ref>),
and by constraining the parameter window as in (<ref>), it holds true   that the solution of the Hartree-Hartree equation (<ref>)
verifies the semi-classical bounds contained in Assumption <ref>.
let us mention here that–while the above discussion is important in the macroscopic regime–in the microscopic regime (i.e. when $\hbar=1$),
the bounds (<ref>)
are automatically satisfied,
independently of the initial data $\omega_0$,
the constraint (<ref>).

The natural  topology in which convergence is expected to hold corresponds to that of trace-class operators.
Our main theorem is the following result.

Assume that the interaction potential satisfies
\int_{ \R^d 	}
\<\xi\>^2	|\hat V (\xi)	| \d \xi < \infty .
Let us consider the following:

* Let $\psi(t)= \exp(-itH/\hbar)\psi $ be the wave function of the system,
with initial data verifying
Assumption <ref>.
Let $\gamma_F(t)$
and $\gamma_B(t)$
be the one-particle reduced density matrices, as defined in (<ref>).

* Let $(\omega(t), \vp(t))$
be the solution of the Hartree-Hartree equation  (<ref>).

Additionally, assume that
the scaling regime is chosen so that the following is satisfied:
for all $\ell \geq 1$ there exists
$ k_\ell  \geq  1$ such that
for all physical parameters $\lambda$, $\hbar$, $N $ and $M$
there holds

\begin{equation}
\label{k condition}
\frac{ \lambda \sqrt N }{\hbar}
\leq
( \hbar M )^{k_\ell }.
\end{equation}

Then, there exists
a constant $C   > 0  $
such that
for all $ t\in \R$
there holds

\begin{align}
\label{thm1 fermion estimate}
\frac{1}{M }
\|
\gamma_F  (t)  - \omega (t)
\|_{\tr}
& \leq
\frac{C}{\sqrt  M }
{  \exp   \bigg[  C   \lambda \sqrt{  \frac{NM}{\hbar }	 }
\bigg( 1 + \sqrt{\frac{\hbar M}{N}}    \bigg)  \exp |t|      \bigg]		}   \ ,  \\
\label{thm 1 boson estimate}
\frac{1}{N}
\|
\, \gamma_B (t)  -
\ket{\vp (t)}
\bra{\vp(t) }
\|_{\tr}
& \leq
\frac{C}{\sqrt  N  }
{  \exp   \bigg[  C    \lambda \sqrt{   \frac{NM}{\hbar }	  }
\bigg( 1 + \sqrt{\frac{\hbar M}{N}}    \bigg)  \exp |t|      \bigg]		}  \ .
\end{align}

The above Theorem provides an explicit convergence rate
from the many-body Schrödinger dynamics
to the solution of the Hartree-Hartree system.
Note that we have chosen not to fix the
parameter regime in the theory–this is in contrast to
most works in the literature.
The reason is that the scaling regime in which Theorem <ref>
provides a reasonable approximation was not known
by the authors in the onset of this investigation.
Our interest then was not to prove an optimal convergence rate,
but to actually find a meaningful scaling regime.

Regarding the previous remark, 	Theorem <ref> contains a meaningful approximation
whenever the argument in the time dependent function is $\O(1)$
with respect to the physical parameters.
Let us   describe two scaling regimes that we regard as interesting.

Microscopic regime.
If one is working in microscopic units, we may set $\hbar  =1 $.
One can then investigate the mean-field regime  in which the number of bosons and
fermions is the same. Namely

\begin{equation}
\lambda  =  \frac{1}{N }  \ , \qquad \hbar =1
\quad \t{and} \quad
N = M
\end{equation}

and we also set $m_F = m_B =2 $ for completeness.
Clearly, the condition (<ref>)
is verified with $k_\ell= \ell $.
In this case, one should regard Theorem <ref>
as capturing the emergence of the mean-field equations

\begin{align}
\begin{cases}
& i    \partial_t \omega
[   -   \Delta
(V   * \rho_B )  ,
\omega
]			\\
& i   \partial_t  \vp
-  \Delta  \vp
(V  *  \rho_F)
\vp   \  ,
\end{cases}
\end{align}

as the leading order term driving the dynamics of the Hamiltonian $H$, for our choice of initial data.
Note that in this case, the semiclassical condition
imposed in (<ref>)
is verified immediately, independently of the structure of the initial data.
However, since $ \mathrm{Tr}  \omega (t) = M $, the above equation  does not yield a non-trivial limit when  $M \rightarrow \infty$.

Macroscopic regime.  In macroscopic units, the value of $\hbar$ becomes small.
As is well-known,
for a system of  confined fermions,
the energy  scale of each particle is
$\hbar^2 M^{\frac{2}{d}}/m_F$.
One is then interested in the regime for which
this scale is  $\O(1)$ – this is the so-called semi-classical limit that has been studied extensively in the literature for systems of interacting fermions.
On the other hand, for bosons the energy per particle has
the scale
$\hbar^2/m_B $.
We can then tune the parameters so that the total energy of the system is balanced.
For instance, we may look at

\begin{equation}
\label{scaling regime}
\lambda = \frac{1}{N}
\ ,
\qquad \hbar =  \frac{1}{M^{\frac{1}{d}}}  \ ,
\qquad  m_B = \hbar \ ,
\qquad m_F =1 \,
\quad \t{and} \quad
N = M^{	 1 +		\frac{1}{d}	} \ .
\end{equation}

It is possible to   check that   condition (<ref>)
is verified
$ k_\ell
\frac{1 + d (2\ell-1 )}{2(d-1) }.
one may readily verify that the condition (<ref>) is satisfied.
This leads to a natural candidate on the initial data for the fermionic component $\omega_0$
that verifies Assumption <ref>;
see Remark <ref> for more details.

§.§ The semi-classical limit

In this subsection, we adopt the macroscopic
scaling regime given by the equations
Let us now motivate the upcoming semiclassical analysis of the
coupled Hartree system.
In what follows, we shall   incorporate
the $\hbar$ dependence  on the solution $( \omega^\hbar, \vp^\hbar)$
of the coupled Hartree-Hartree equation (<ref>).
We start by noting that one of the main consequences of the scaling regime (<ref>)  is that
$\lambda M = \hbar $.
Hence, the  Hamiltonian for the boson field  $\vp_t^\hbar$ is
proportional to $\hbar$, i.e. it has the form
$	\hbar ( -  (1/2) \Delta +   V*\rho_F^\hbar ) 	$–it follows that factors of  $\hbar$ cancel out in the equation.
the solution of the   coupled Hartee equation can now be analyzed semi-clasically, in the limit
$\hbar \downarrow 0$.
Indeed,  for $t\in\R $
we   consider the Wigner transform of the fermionic density matrix

\begin{equation}
\label{wigner transform}
f^\hbar  (t)
\equiv
W^\hbar[\omega^{\hbar} (t) ]
\quad
\t{where}
\quad
W^\hbar[\omega ] ( x, p )
\equiv
\frac{ 1}{(2 \pi)^{\frac{d}{2}}}
\int_{\R^d }
\omega
\Big( x +   \frac{y}{2} , x -  \frac{y}{2}
\Big)
e^{		 - i  \frac{y\cdot p }{ \hbar }	}		 \d y   \ ,
\end{equation}

where $(x,p  ) \in      \R^d  \times \R^d .$
Heuristically, the pair $( f^\hbar , \vp^\hbar)$   converges to a solution $(f, \vp) \in L^1_+(\R^{2d})\times L^2(\R^d)$

\begin{align}
\label{vlasov hartree}
\begin{cases}
& 		(\partial_t  + p \cdot  \nabla_x + F_B(t,x)  \cdot  \nabla_p   ) f  = 0 	\\
&		  i \partial_t \vp
- \frac{1 }{2} \Delta   \vp +   ( V*\rho_F   )   \vp 	\\
&  (f,\vp) ( 0 )  = (f_ 0 , \vp_0)
\in L^1_+(\R^{2d})\times L^2(\R^d)
\end{cases}
\end{align}

where $F_B(t,x) \equiv
- \int\nabla V(x-y ) | \vp(t,y)|^2 \d y
$\rho_F 	(t,x) \equiv  \int_{	\R^d	} f (t,x,p ) \d p  	. $
We shall refer to (<ref>) as the Vlasov-Hartee equation.

As we have mentioned previously,
we are interested in the case in which the initial system is
at zero temperature.
In this situation,
Thomas-Fermi theory suggests that
one cannot expect  the initial data of
the classical fermion subsystem to have plenty of regularity (see for reference
In what follows, we introduce distances
which we will use throughout this article.
In particular, they will be necessary in our analysis of convergence to the Vlasov-Hartee system
in this context of “low regularity”.

* Wasserstein distance.
Given $n \in \mathbb N$, we denote the  $n$-th Wasserstein distance between probability measures
$\mu,\nu \in \mathscr{P}_n (\R^{2d})$ by

\begin{equation}
W_n (\mu , \nu )
\equiv
\Big(
\inf_{ \pi }
\int_{\R^{2d} \times \R^{2d}}  | z - z ' |^n  \pi( \d z \otimes  \d z  ' )
\Big)^{\frac{1}{n}}
\end{equation}

where the inifimum is taken over all couplings of $\mu$ and $\nu$, i.e.
probability measures $\pi \in \mathscr{P}(\R^{2d} \times \R^{2d})$
with  first marginal $\mu$, and second marginal $\nu$.

* Fourier-based norms.
Given $s \in \R$,
and $ g : \R^{2d} \rightarrow \C $
we introduce the following
Fourier-based norm

\begin{equation}
\label{fourier based norm definition}
| g |_s
\equiv
\sup_{ \zeta \in \R^{2d}}
(1 +  |\zeta|)^{-s} \ |  \hat g ( \zeta  )  |   \  .
\end{equation}

In applications, we take $s \geq 0 $. Hence,  we also regard $ |\cdot|_s$
as a negative Sobolev norm.

Our assumption for the initial data now reads.
[Hartree initial data]

The pair $(\omega_0^\hbar ,\vp_0^\hbar) \in \mathscr{L}^1 (L^2(\R^d))\times L^2(\R^d)$ satisfies the following conditions.

* $\omega_0^\hbar \in \mathscr{L}^1 ( L^2(\R^d)  )$
satisfies $0  \leq \omega_0^\hbar  = (\omega_0^\hbar)^* \leq 1$,  $\t{Tr} \omega_0^\hbar = M $
$\tr \omega_0^\hbar( - \Delta + x^2  ) < \infty  $.
Further, 	we assume that there exists
a real-valued
$      f_0 \in L^1   (\R^{2d})$
such that :

(1.1) $ 0 \leq f_0(x,p) \leq 1$
$\int_{	\R^{2d}	} f_0 (x,p ) \d x \d p  = 1$.

There are finite second moments:
$f_0 \in \mathscr P _2 (\R^{2d})$.

\lim_{\hbar \downarrow 0}
| f _0 - f_0^\hbar |_{1}
= 0
where  $f_0^\hbar = W^\hbar [\omega_0^\hbar]$ .

* There exists $\vp_0 \in L^2  (\R^d)$
with $\| \vp_0	\|_{L^2}=1 $
such that
$  \lim_{\hbar \downarrow 0} \|   \vp_0 - \vp_0^\hbar		\|_{L^2}=0.$

Let us note that in Assumption <ref>
there are no regularity requirements
on the limits of the sequence of initial data.
Of course, this comes with a price.
First, we shall need  the interaction potential to be  at least $V \in C^{1,1} (\R^d;\R)$.
Second, the metric that we use to measure the distances between fermion densities is rather weak.
Namely, it involves testing over functions $h (x,p)$
for which both the integrals
$ \int_{	\R^{2d}	}  \< \zeta\>^2 |\hat h (\zeta)|^2 \d \zeta  $
$\int_{	\R^{2d}	} |\zeta | \,  |\hat h (\zeta)| \d \zeta $
are finite.
Third–compared to similar results in the literature–we need
two moments in phase space, rather than only one.

In Assumption <ref> we require  that  $f_0^\hbar \rightarrow f_0$
with respect to the negative Sobolev norm $|\cdot|_1$.
Two comments are in order.

This assumption  on the initial data
can be verified
for examples  of non-interacting Fermi gases,
that arise
as minimizers of variational problems in the presence of an external trap $V_{ \t{ext} }(x)$.
More precisely, let  $\omega_0^\hbar $
be  the one-particle reduced density matrix of an approximate ground state $\psi_F^\hbar $ of the   minimization problem

\begin{equation}
E (M) \, = \,
\t{inf}_{\bigwedge_{i=1}^M L^2(\R^d )}
\
\sigma
\bigg[
\sum_{i=1}^M  - \hbar^2  \Delta_{x_i}    + V_{\t{ext}}   (x_i)
\bigg]  \  , \qquad \hbar = M^{ -1 /d } \ .
\end{equation}

Without loss of generality,
one may assume that  $\omega_0^\hbar$ is an approximate minimizer of the associated
Hartree-Fock problem

\begin{align}
E_{HF}( M )
\inf
\Big\{
\tr  ( \omega [- \hbar^2 \Delta + V_{\t{ext}} (x) ] )
\ :  \
\omega = \omega^2 = \omega^*  , \ \tr \omega = M
\Big\}  \ .
\end{align}

can be assumed to be
an orthogonal projection, i.e.
$\omega_0^\hbar = (\omega_0^\hbar)^2$,
which  is equivalent to $\psi_F^\hbar$ being
a Slater determinant.
It has been proven in  <cit.>
that,   as $\hbar \downarrow0$ (and, up to extraction of a subsequence),
the Wigner transform $f_0^\hbar = W^\hbar[\omega_0^\hbar]$
converges   in a weak sense to  the function  $f_0 (x,p) = \1 (     p^2 \leq C_{TF}\rho(x)^{2/d}  )$.
Here,   $C_{TF}$ is the Thomas-Fermi constant,
$\rho$ is the minimizer of  the associated Thomas-Fermi problem

\begin{eqnarray}
\mathcal E (\rho) = \frac{d}{d+2} \int_{\R^d } \rho(x)^{1 + \frac{2}{d}} \d x
+ \int_{\R^d }V_{\t{ext}}  (x)    \rho (x) \d x
\end{eqnarray}

with constraints $ \rho (x) \geq 0 $, $\int \rho(x) dx =1 $ and $\rho \in L^1 \cap L^{1 + \frac{2}{d}}$.
In particular,
the $\hbar\downarrow0 $ convergence
can be shown to hold
with respect to the negative Sobolev norm  (<ref>) as well.
In other words,   there holds $\lim_{\hbar\downarrow 0 }   |        f_0^\hbar - f_0   |_1  =0 $.
This is proved by the first author in upcoming work [13].

The above discussion  should be compared with the
$L^1$-norm convergence   considered in
for the  initial data, in the context of interacting Fermi gases.
While their conclusion is strictly stronger–that is, stronger mode of convergence–to the authors best knowledge
the only examples in $\R^d$ for which the $L^1$ convergence  has been verified correspond to coherent states.
these are not examples  of  zero temperature states
(i.e. orthogonal projections).
We believe that there is value in our approach
since–as the examples considered in (i) arise form orthogonal projections–we are able to put together Theorem <ref>
and our next result.
In particular, with this approach   we obtain  a
quantitative convergence
from the Schrödinger to the Vlasov-Hartee dynamics,
in the situation of low regularity or–put differently–the zero temperature situation.

Our main result concerning the semi-classical limit of the coupled Hartree equations is the following theorem.

Assume that the interaction potential satisfies
\int_{ \R^d 	}
\<\xi\>^2 |\hat V (\xi)	| \d \xi < \infty .
Let us consider the following:

* Let $(\omega^\hbar , \vp^\hbar)$
be the solution of the Hartree-Hartree  system
with intial data $(	 \omega_0^\hbar , \vp_0^\hbar	)$ satisfying Assumption <ref>.
Denote by 		  $f^\hbar (t  ) = W^\hbar[\omega^\hbar (t)]$
its Wigner transform.

* Let $(f, \vp)$ be the solution to the Vlasov-Hartree
system     (<ref>),
with initial data $(f_0, \vp_0). $

Then,  there exists $C>0$ such that
for all times
$ t \in \R $
and test functions $h : \R^{2d} \rightarrow \C$,
the following inequalities hold true

\begin{align}
\nonumber
|	 \<   h , ( f_t   -  f^\hbar_t  )  \>		|
& \leq
\|	  \<\zeta\>  \hat h  	\|_{L^1}
\Big(
|	   f^\hbar_0  -   f_0  			|_1
\hbar
\Big)
\|	   |\zeta|  \hat h  	\|_{L^2}
\Big(
\|	  \vp_0 -   \vp_0^\hbar 	\|_{L^2}
+ 		\hbar^{1/2}
\Big)   \, \\
\| \vp_t    -  \vp^\hbar_t		\|_{L^2 }
& 	\leq
\Big(
|	   f^\hbar_0  -   f_0  			|_1
+ \hbar
\Big)
C_1 (t)
\Big(
\|	  \vp_0 -   \vp_0^\hbar 	\|_{L^2}
\hbar^{1/2}
\Big) \  .
\end{align}

Here, we are denoting
$C_1(t) =
C\exp(C t^2 )
$C_2(t) = C\exp ( C \exp  C  | t |  ) $.

The above result gives an explicit convergence rate from the
Hartree-Hartree to the Vlasov-Hartee dynamics.
Of course, this is  not the optimal convergence rate in $\hbar$, which
we believe should be $\mathcal O (\hbar)$
on the right hand side.
In this work, we have not tried to optimize this rate.
Indeed, our main goal was to identify
the leading order equations that drive the effective dynamics
of the Bose-Fermi mixture,
which Theorem <ref> appropriately does.
In a similar spirit, we have not tried to optimize the growth-in-time
of the constants  involved in our estimates.

Here we have formulated our theorem in terms of test functions.
Alternatively,  it can be formulated in terms of the norm

\begin{equation}
\|		f \| \equiv \sup
\big\{
\<  h , f \>
h\in \mathscr{S} (\R^{2d})
\|   \< \zeta \> \hat h		\|_{L^1} \leq 1
\	    \t{ and } \
\|    | \zeta| \,  \hat h		\|_{L^2} \leq 1
\big\}  ,
\end{equation}

which is
strictly weaker than the norms
\|		\cdot \|_{\dot H ^{-1}}$
$ | \cdot |_{1 }$.

§.§ Strategy of the proofs

Let us  outline the proofs of our main results, Theorem <ref> and <ref>.

The proof of Theorem <ref>
consists of the study of an appropriate fluctuation dynamics.
For  gases of interacting bosons
the approach was  first carried out in [55],
whereas for gases of interacting fermions
the approach  was employed  in [9].
The difficulty  of tackling the Bose-Fermi
mixture lies in how to properly combine these two approaches.
In the present paper, we adapt the approach of studying fluctuation dynamics for the problem at hand. Namely, we introduce in Section <ref>
the formalism of Second Quantization on Fock space $\F$.
In this formalism,
coherent states describing Bose-Einstein condensates
are parametrized by a Weyl operator $\W[\sqrt N \vp(t)]$,
whereas degenerate Fermi gases are implemented
by a particle-hole transformation $\calR[\omega(t)]$; see Sections <ref>  and <ref> for more details.
In Section <ref>
we then study the dynamics of fluctuations
around the tensor product of these states.
Roughly speaking,  the problem  is then reduced
to estimating
the “number of excitations" outside of $(\vp(t), \omega(t))$.
We implement this point of view by introducing a new  unitary transformation on $\F$,
denoted by $\U(t,s)$ and defined in (<ref>).
Its understanding is  fundamental in our analysis, and leads to the number estimates contained in Theorem <ref>.
The proof of these estimates
is based on the analysis of  its infinitesimal generator, which has the form (see Lemma <ref> for details)

\begin{equation}
\L(t)
\d \Gamma_F [h_F(t) ] \otimes \1
\1 \otimes \d \Gamma_B [h_B(t)]
+ \lambda \sqrt N \L_{2,1}(t) + \lambda \L_{2,2}(t)
\end{equation}

Here, the difficulty lies in controlling the terms $\L_{2,1} (t)$ and $\L_{2,2}(t)$,	which do not commute with
particle number operators and can potentially increment the number of fluctuations.

The proof of Theorem <ref>  is essentially divided in two steps.
First, we rely  on techniques developed in
[8]
to understand the stability of the Hartree-Hartree equation (<ref>)
with respect  to the metric in $\mathscr{L}^1( L^2( \R^d)) $
that is induced by the norm $|\cdot |_s$,
defined in (<ref>).
using some recently developed tools
from Quantum Optimal Transporation,
we are able to
show that
the convergence from the Hartree-Hartree  (<ref>)
to the Vlasov-Hartree (<ref>) dynamics
can be controlled in
the negative Sobolev space $\dot H ^{-1}$.
These tools include the introduction in
[36]
of the quantum analogue of the classical Wasserstein distance between two probability measures,
which have been later developed and applied to the analysis of single species many-particle systems in a series of papers–see for instance
[37, 38, 39, 46].
One of the main advantages
of   these techniques is    the fact that they require no regularity
on the initial data under consideration.
This is compatible with Assumption <ref>, in which we assume our initial data corresponds
to a zero temperature state–for fermionic systems,
the $\hbar\downarrow 0 $ limit of the  Wigner function of
these states is expected to be of Thomas-Fermi type,
which fails  for instance to be in $W^{1,1}$.

§.§ Physical discussion of our results

Let us recall that we have chosen to work with dimensionless variables.
Thus, the physical parameters  $\hbar$, $m_F$, $m_B$ and $\lambda$
that appear in the Hamiltonian (<ref>)
correspond to the numerical values
that the original parameters take,
when measured with respect to
some chosen scales of length $\ell_0$,  time $\tau_0$ and mass $m_0$.
In other words, denoting by $E_0 = m_0 \ell_0^2 \tau_0^{-2}$
the characteristic energy scale of the system, we have

\begin{equation}
\hbar = \frac{\hbar_0}{E_0 \tau_0},
\qquad
m_B = \frac{m_{B,0}}{m_0},
\qquad
m_F = \frac{m_{F,0}}{m_0},
\quad
\t{and}
\quad
\lambda = \frac{\lambda_0}{E_0}
\end{equation}

$\hbar_0 = 10^{-34}{\tt kg \, m^2 \, s^{-1}},$
$m_{F,0}$ and $m_{B,0}$ are the original masses
of fermions and bosons,
and $\lambda_0$ is the  original interaction strength.
Recall that in the scaling regime (<ref>)
we have considered $m_F =1$ and $m_B = \hbar$.
In particular, let us set
$m_0  =  m_{F,0}$ as our   mass scale.
Hence, $m_F=1 $
and we find that $m_B = \hbar $ is equivalent to
the dimensionless equation

\begin{eqnarray}
\frac{m_{B,0}}{m_{F,0}}
\frac{\hbar_0}{ m_{F,0} \ell_0^2 \tau_0^{-1} } \ .
\end{eqnarray}

Let us now comment on the physical validity of this regime.
Amongst the lightest bose atoms that have been
used to create Bose-Einstein condensates lies Helium-4, with a mass of $ \sim  4 \,  \tt{Da}$.
On the other hand, Ytterbium-171 has been cooled down to form degenerate Fermi gases,
and has a mass of
$ 171  \, \tt{Da} \sim 3 \times 10^{-25} \,  \tt{kg} $.
The ratio between these two masses then gives a value
$ m_B   \sim  2 \times 10^{-2}$.
In particular, in an experiment that has been confined
to length scales of $\ell_0 \sim  10 \tt{\mu m}$, and is observed for $ \tau_0 \sim  10^{-2} \tt{s}$,
gives a value for the effective Planck's constant of
$ \hbar \sim  3 \times 10^{-2} $–comparable to the ratio of  the masses of the Bose-Fermi mixture.

We also note that some small mass limits have been considered in theoretical physics in the search for evidence of ultralight bosons in cosmology.
For example, in [50], the authors consider the Schrödinger-Poisson system, and compare it to the Vlasov-Poisson equation in the classical limit with $\hbar/m_B \rightarrow 0$.
In their analysis, they consider bosonic masses $m_{B,0}$ in the range of $10^{-22}{ \tt  eV / c^2 } \sim 10^{-58} \tt{kg}$, corresponding to axion dark matter.

§ SECOND QUANTIZATION I: PRELIMINARIES

It is convenient to study the Hamiltonian $\eqref{hamiltonian HNM}$
in the second quantization formalism.
we allow the number of particles to fluctuate
and thus consider the
Hilbert space composed of the corresponding Fock spaces.
Namely, we let

\begin{equation}\label{fock space definition}
\F \equiv \F_F \otimes \F_B\ ,
\quad \t{where} \quad
\mathscr{F}_F \equiv \C \oplus \bigoplus_{n=1}^\infty
L^2_a( \R^{dn} )
\quad \t{and} \quad
\mathscr{F}_B \equiv \C \oplus \bigoplus_{n=1}^\infty
\end{equation}

are the fermionic and bosonic Fock spaces, respectively.
Here,  $L^2_a $
and  $ L^2_s  $
correspond to $L^2$-functions  that are antisymmetric and symmetric, respectively, with respect to permutations of  the particle   variables.

On $\F_F$ we introduce fermionic creation- and annihilation operators $a_x$ and $a_x^*$
as the operator-valued distributions
that satisfy the Canonical Anticommutation Relations (CAR)

\begin{equation}
\{  a_x , a_{x'}^*	\} = \delta(x- x '  )
\qquad
\t{and}
\qquad
\{ a_x ,a_{x' } \}
\{ a_x^* ,a_{x' }^*  \}  =0
\quad x,x ' \in \R^d  \,  ,
\end{equation}

where $ \{ \cdot, \cdot \} $ is the anticommutator bracket.
Similarly, on $\F_B$ we introduce bosonic creation- and annihilation operators
$b_y$ and $b_y^*$
as the operator-valued distributions
that satisfy the Canonical Commutation Relations (CCR)

\begin{equation}
[b_y,  b_{y ' }^*]= \delta( y - y'  )
\qquad
\t{and}
\qquad
[  b_y  ,b_{y' }   ]
[  b_y^* ,b_{y'}^*  ]   =0
\quad y,y ' \in \R^d  \,  ,
\end{equation}

where now $[ \cdot  , \cdot ] $ is the commutator bracket.
We shall denote by
$\Omega_F = (1, \0)$
$\Omega_B = (1,\0 )$
the vaccum vector in each space,
and by $\Omega \equiv \Omega_F \otimes \Omega_B$
the vaccum state of the combined system.

In this setting, the  many-particle Hamiltonian
introduced in the previous section
can be written in terms of
creation- and annihilation-
in the following form

\begin{equation}
\label{hamiltonian}
\H =
\frac{\hbar^2 }{2m_F}
\int_{  \R^d }
a_x^* (  - \Delta_x )
a_x \d x
\frac{\hbar^2 }{2m_B}
\int_{  \R^d }
b_y^* (  - \Delta_y )
b_y \d y
+ \lambda
\int_{ \R^{2d} }
V(x-y) \
a_x^*a_x   b_y^* b_y
\ \d x \d y  \ ,
\end{equation}

where we do not display explicitly the tensor product symbols.
As for the dynamics,
we introduce
the time evolution
of the second quantized system

\begin{equation}
\Psi(t)
\equiv
\exp\big(
- i t \H / \hbar
\big)
\Psi(0)  \ , \qquad \forall t \in \R \ .
\end{equation}

Since the Hamiltonian $\H$
is quadratic and diagonal in creation-
and annihilation operators,
it commutes
with the \textit{fermionic}
and \textit{bosonic}
number operators
\begin{equation}
\N_F
\equiv
\int_{	\R^d	} a_x^* a_x
\d x
\qquad
\t{and}
\qquad
\N_B
\equiv
\int_{	\R^d	}
b_y^* b_y \d y \ .
\end{equation}
Consequently, if
δ_n,N δ_m,M ψ(0)$,
then for all $ t ∈$ it holds that
\begin{equation}
\Psi(t)^{(n,m)}
\delta_{n,N} \delta_{m,M} \psi(t) \   ,
\end{equation}
where $ψ(t)$
is the state
corresponding to the
$(N+M)$-particle system, defined in \eqref{wave function}.
Most importantly,  one may verify that
the following relations hold true
for the one-particle reduced density matrices
\begin{equation}
\label{density matrices second quant}
\gamma_F(t ; x_1 ,x_2)
\<
\Psi (t) ,
a_{x_2}^* a_{x_1}\otimes
\Psi(t)
\>_\F
\quad
\t{and}
\quad
\gamma_B (t ; y_1 ,y_2)
\<
\Psi (t) ,
\1 \otimes
b_{y_2}^* b_{y_1}
\Psi(t)
\>_\F
\   .
\end{equation}
The equations given in \eqref{density matrices second quant}
are  the starting point in the proof of Theorem \ref{thm1}.

\vspace{2mm}

In the rest of this section,
we introduce   preliminaries
that we will need to prove Theorem \ref{thm1}.
In Subsections \ref{subsection fermions} and \ref{subsection fermions}
we give a more detailed account   of the second quantization formalism for both fermions and bosons.
Our goal here is not   to be thorough, but to collect basic results and fix the notation that will be used throughout the article.
The reader is refered to the book [11]
and the lecture notes
[57]
for more details.

\subsection{Fermions}
\label{subsection fermions}
Throughout this subsection, we will write the Fermionic Fock space
$_F$  as follows
\begin{equation}
\mathscr{F}_F
\bigoplus_{n=0}^\infty
\mathscr{F}_F^{(n)}
\quad
\t{where}
\quad
\F_F^{(0)} \equiv \C
\quad
\t{and}
\quad
\mathscr{F}_F^{(n)} \equiv 	 L^2_a( \R^{dn} ) , \ \forall n \geq 1  ,
\end{equation}
where  $	L^2_a( ^dn )  $
to the space of $L^2$ functions that are    antisymmetric
with respect to the permutation of the particles position variables.
%	That is
%	\begin{equation}
%		\label{L2 antisymmetric}
%		L^2_a( \R^{dn} )    =
%		\{  \psi \in L^2(\R^{dn}) :
%		\psi ( x_{\pi (1)}    , \ldots, x_{	 \pi ( n) 	}	   )
%		=
%		\t{sgn}(\pi)
%		\psi (	 x_1, \ldots, x_n	)  \forall \pi \in S_n \} \ .
%	\end{equation}
The space $_F$  becomes a Hilbert space when endowed  with the inner product
\begin{equation}
\< \Psi_1 ,  \Psi_2 \>_{\F_F}
\equiv
\sum_{n=0}^\infty
\langle
\,
\Psi_1^{(n) } , \Psi_2^{(n)}
\,
\rangle_{\F_F^{(n)}}  \  ,
\qquad \forall \Psi_1 ,\Psi_2 \in \F_F \ .
\end{equation}
On the Fock space $_F$ one introduces the smeared-out
creation- and annihilation operators
as follows.
Given $f ∈L^2(^d) $,
we let   $a^*(f)$ and $a(f)$
be defined   for  $Ψ∈_F$ as
\begin{align}
\big(		 a^*(f)  \Psi \big)^{(n)}
(x_1, \ldots, x_n)
& \equiv
\frac{1}{\sqrt n }
\sum_{i=1}^n
f( x_i) \Psi^{(n-1)}
(x_1 , \ldots, x_{i-1} , x_{i+1} , \ldots, x_n )	\ , \\
\big(	 a  (f)	 \Psi \big)^{ (n)}
(x_1, \ldots, x_n)
& \equiv
\sqrt{n+1 }
\int_{\R^d }  \overline{f(x)} \Psi^{(n+1)} (x , x_1 , \ldots, x_n) \d x  \ .
\end{align}
In particular, they satisfy
the following version of the CAR
\begin{equation}
\{ a(f)  , a^*(g) \}    = \< f,g\>_{L^2}
\quad
\t{and}
\quad
\{	  a^\hash (f) , a^\hash (g)		 \} = 0
\end{equation}
where we recall ${ ·, ·}$ stands for the anticommutator.
In particular,  it is easy to see that
the CAR   turns them into bounded operators, with norms
\begin{equation}
\|	  a^*(f)	\|_{B (\F_F)} =
\|	  a(f)	\|_{B (\F_F)}  =
\|	 f\|_{L^2 } \ .
\end{equation}
Let us finally mention that the connection with the operator-valued distributions $a_x^*$ and $a_x$ is by means of the formulae
\begin{equation}
a^*(f) = \int_{\R^d }  f(x)  a_x^* \d x
\qquad
\t{and}
\qquad
a(f) = \int_{	\R^d	} \overline{f(x)} a_x \d x   \  ,
\qquad f \in L^2(\R^dw)
\end{equation}

\subsubsection{Fermionic operators on $\F_F$}
Given a closed linear operator
$Ø: (Ø) ⊂L^2(^d ) →L^2(^d )$,
we consider its second quantization
$Γ̣_F [Ø]$
as the diagonal operator
on $_F$, defined as
\begin{equation}
( 	\d\Gamma_F [\O]  )^{ (n)}
\equiv
\sum_{ i = 1 }^n
\mathbb{1}^{ i -1 } \otimes \O \otimes \mathbb{1}^{ n - i }  ,
\ n \geq 1
\qquad
\t{and}
\qquad
( 	 	\d\Gamma[\O]  )^{ (0)} \equiv 0  \   ,
\end{equation}
initially   on tensor products of elements of $(Ø)$, and then closed.
In most cases of interest $Ø$ is bounded (or even trace-class)
and the domain of $ Γ̣_F [Ø]$ is contained in $(_F)$--the only exception will be the  Laplacian $Ø= - Δ$, in which case
$ħ^2/2m_F Γ̣_F[-Δ]$ is the
associated kinetic energy of the system.

\vspace{1mm}

The reader should be aware that, at least formally,
if   $Ø$ has an operator kernel
$Ø(x,x')$, then  one may write in terms of creation- and annihilation operators
\begin{equation}
\d\Gamma_F  [\O ] = \int_{	\R^{2d}	} \O (x,x')      a_{x'}^*a_x     \d x \d x' \ .
\end{equation}
In this context, one of the most important observables in second quantization corresponds to the
the {fermionic number operator}.
It is defined as   the second quantization of the
the identity operator  on $L^2(^d)$, which
has the distributional  kernel
$ (x,x') = δ(x-x')$.
\begin{equation}
\label{fermion identity N}
\N_F =  \bigoplus_{n=0}^\infty  n = \d \Gamma_F [\1 ] = \int_{	\R^d	} a_x^*a_x \d x  \ .
\end{equation}

\vspace{1mm }

Let us now  collect some basic results concerning estimates
for the second quantization of operators in fermionic Fock space
in the following lemma.
For a proof, we refer the reader to
\cite[Lemma 3.1]{BenedikterPortaSchlein2014}.

\begin{lemma}[Estimates for fermionic operators]
\label{lemma fermion estimates}
Let $\O: L^2(\R^d) \rightarrow L^2(\R^d)$ 	 be a bounded operator.
Then, the following holds true.

\begin{enumerate}[leftmargin=1cm]
\item
For all $\Psi , \Phi \in \calD(\N_F)$
\begin{align}
\|	 \d \Gamma_F (    \O	  	) \Psi	 \|_{\F_F}			&
\leq	\|	   \O  	\|_{B(  L^2 )}
\|	  \N_F  \Psi 	\|_{\F_F}		 \ , \\
| 		\< \Psi , \d  \Gamma_F  (	 \O 	)\Phi \>_{\F_F}	|
\leq	\|	   \O  	\|_{B(  L^2 )}
\| 			\N_F^{\frac{1}{2}}	 		\Psi 	\|_{\F_F}
\| 				\N_F^{\frac{1}{2}}				 \Phi 	\|_{\F_F} \ .
\end{align}
\item
If $\O$ is Hilbert-Schmidt   with kernel $\O(x,y)$, then for all $\Psi \in \calD (\N_F^{\frac{1}{2}})$
\begin{align}
\|		 \d \Gamma_F (\O)	 \Psi 	\|_{\F_F}		& \leq 		\|	 \O	\|_{HS}		\|		 \N_F^{\frac{1}{2}}	\Psi 	\|_{\F_F}		\\
\|		\int_{	 \R^d \times \R^d 	}	 \O(x_1 , x_2 )  a_{x_1} a_{x_2} \d x_1 \d x_2  		\|_{\F_F}				 & \leq 		 \|	 \O	\|_{HS}		\|		 \N_F^{\frac{1}{2}}	\Psi 	\|_{\F_F}		\\
\|		\int_{	 \R^d \times \R^d 	}	 \O(x_1 , x_2 )  a^*_{x_1} a^*_{x_2} \d x_1 \d x_2  		\|_{\F_F}				 & \leq 		 \|	 \O	\|_{HS}		\|	( 	 \N_F+\1)  ^{\frac{1}{2}}	\Psi 	\|_{\F_F}	 \ .
\end{align}
\item
If $\O$ is trace-class with kernel $\O(x,y)$, then for all $\Psi \in \F_F$
\begin{align}
\|		 \d \Gamma_F  (\O)	 \Psi 	\|_{\F_F}		& \leq 		\|	 \O	\|_{Tr}		\|		 	\Psi 	\|_{\F_F}			\\
\|		\int_{	 \R^d \times \R^d 	}	 \O(x_1 , x_2 )  a_{x_1} a_{x_2} \d x_1 \d x_2  		\|_{\F_F}			 & \leq 			\|	 \O	\|_{Tr}		\|		 	\Psi 	\|_{\F_F}			\\
\|		\int_{	 \R^d \times \R^d 	}	 \O(x_1 , x_2 )  a^*_{x_1} a^*_{x_2} \d x_1 \d x_2  		\|_{\F_F}			 & \leq 			\|	 \O	\|_{Tr}		\|		 	\Psi 	\|_{\F_F}		 \ .
\end{align}
\end{enumerate}
\end{lemma}

\subsubsection{Particle-hole transformation}
\label{subsubsection PH}

In this subsection we introduce a class of
Bogoliubov transformations
on Fock space
that we will use in the proof of Theorem \ref{thm1};
they will be useful in
the number of fluctuations
outside of a degenerate Fermi gas.

\vspace{2mm}
More precisely,
let us consider
to be a   rank-$M$ orthogonal projection
on $L^2(^d)$.
Thus, there exists an orthonormal basis
${  ϕ_i }_i=1^∞⊂L^2(^d)$ such that
\begin{equation}
\textstyle
\omega = \sum_{i=1}^M \ket{\phi_i}		\bra{\phi_i}   \ .
\end{equation}
We   introduce a       map on Fock space
$[ω] : _F →_F$,
which we shall refer to as a
\textit{particle-hole transformation} associated to  $ω$.
We define it
according to its action on creation-
and annihilation
operators as follows
\begin{equation}
\label{R omega}
\calR ^*[\omega] a^*(\phi_i) 	\calR [\omega]
\equiv
\begin{cases}
a^* (\phi_i) \, ,  & i \leq M \\
a(\phi_i)  \, , & i> M
\end{cases} \ ,
\end{equation}
and its action on the vaccum
$	[ω] Ω_F
≡a^*(ϕ_1) ⋯a^*(ϕ_M) Ω_F .$
Note that since
the span
of vectors of the form
$   a^*(ϕ_i_1)   ⋯a^*(ϕ_i_n    )Ω_F$
is dense in $_F$,
the above prescription
completely determines $[ω]$.

\vspace{2mm}

Let us now collect additional properties of the map $[ω]$
In order to state them, we need to introduce some important notation.
Indeed, we consider the operators
on the one-particle space
$u,v ∈B(L^2(^d))$
defined as
\begin{equation}
\label{definition u and v}
\textstyle
u \equiv \1 - \omega
\qquad
\t{and}
\qquad
v= \sum_{ i = 1 }^M \ket{\overline{\phi_i}} \bra{\phi_i} \ .
\end{equation}
The following properties are recorded in the following Lemma.
We refer the reader to [9] for more details.

\begin{lemma}
[Properties of $\calR$]
\label{lemma fermion properties}
Let $\omega$, $u$, $v$ and $\calR[\omega]$ be as above.
Then, the following statements hold true.

\begin{enumerate}
\item
$\calR[\omega]$ is a unitary transformation on $\F_F$, and $\calR^*[\omega] = \calR[\omega]$.

\item
We denote  $u_y (x) \equiv u(x,y)$   and $v_y(x) \equiv v(x,y)$.
Then, for all $x \in \R^d$
\begin{align}
\calR ^*[\omega]
\ a^*_x  \
\calR [\omega]
& 	=
a^* ( u_{x})
a( \overline{  v_{ x } } )	  \ ,   \\
\calR ^*[\omega]
\ a_x  \
\calR [\omega]
& 	=
a ( u_{ x })
a^*( \overline{  v_{ x} } )	 \ .
\end{align}

\item
For  all $x,y  \in \R^d$  there holds
\begin{equation}
\< 	\calR [\omega]  \Omega_F    ,  \, a_{y}^*  \, a_x \,  	\calR [\omega]   \Omega_F    \>_{\F_F}
=  \omega (x,y)   \  .
\end{equation}
In  words, the one-particle reduced density matrix of 	$\calR [\omega]  \Omega_F $ corresponds to

\item
$u^* = u^2 = u \ $
$ \ v^* = \overline v$.

\item
$u^* u + v^*v = \1  \ $
$ \ u \overline v = v u = 0 . $

\end{enumerate}

\end{lemma}

\begin{remark}
The unitary map
$\calR[\omega]$ is an example of the implementator of a Bogoliubov
transformation--that is, a map on Fock space that preserves the
Canonical Anticommutation Relations.
More precisely, consider the maps
\begin{equation}
\label{nu}
\nu
\begin{pmatrix}
u & \bar v \\
v & \bar u
\end{pmatrix} \
\qquad
\t{and}
\qquad
A(f,g)  \equiv  a(f) + a^* (\bar g)
\end{equation}
for all $f,g \in L^2(\R^d).$
In this context, $\nu$ is a \textit{Bogoliubov transformation}.
Namely,  it holds true that
for all $f,f_1 ,f_2 , g \in L^2(\R^d): $
\begin{align}
\{  A(\nu (f_1, g_1) , \nu (f_2, g_2) )   \} & = \{ A(f_1, g_1) ,    A(f_2 , g_2)   \}  \ , \\
A^*(\nu (f,g)    )  & =  A( \nu (\bar g , \bar f))   \, .
\end{align}
$\calR[\omega]$ implements $\nu$
on the Fock space $\F_F$,
in the sense that
for all $f,g\in L^2(\R^d)$
\begin{equation}
\label{conjugation 1}
\calR ^* [\omega ]A(f,g) 	\calR  [\omega]
A( \nu (f,g) )
\end{equation}
\end{remark}

Let us note that while the notion of Bogoliubov transformations is quite general,
in the physical situation at hand
it is sufficient to consider   particle-hole transformations,
which has an explicit representation.
This is because the initial state we consider  is a pure state
corresponding to a Slater determinant of $M$ particles.
its one-particle reduced density matrix is
a rank-$M$ orthogonal projection.
The situation is quite different in the positive temperature case,
when states are mixed, and no longer orthogonal projections.

%Let us present here a class of unitary transformations
%between Fock spaces that will be crucial in our analysis.
%Given any two  $ f , g \in L^2 (\R^d)$, we  introduce the following  two
%mutually adjoint operators in $\F_F  $
%	A(f,g)  \equiv  a(f) + a^* (\bar g)
%	\qquad
%	\t{and}
%	\qquad
%	A^*(f,g)  \equiv  a^*(f) + a (\bar g)   \ .
%We say that a  linear map
%$\nu :  L^2 (\R^d ) \oplus  L^2 (\R^d ) \rightarrow  L^2 (\R^d ) \oplus  L^2 (\R^d ) $ is   a \textit{Bogoliubov transformation} if
%the two following conditions are satisfied
%	\item $\{  A(\nu (f_1, g_1) , \nu (f_2, g_2) )   \} = \{ A(f_1, g_1) ,    A(f_2 , g_2)   \}$
%	\item $  A^*(\nu (f,g)    )  =  A( \nu (\bar g , \bar f)) $ .
%for all $f,f_i,g,g_i \in L^2(\R^d)$.
%Equivalently, one may verify that
% $\nu$ is a Bogoliubov transformation if there exist linear maps $u,v : L^2(\R^d) \rightarrow L^2 (\R^d)$ such that
%	\nu
%	=
%	\begin{pmatrix}
%		u & \bar v \\
%		v & \bar u
%	\end{pmatrix}
%and $u^* u + v^*v =1 $ and $u^* \bar v + v^* \bar u = 0 $.
%Here, $\bar u$ means  complex
%conjugation of the operator kernel.
%We say that $\nu$ is \textit{implementable} in $\F_F$ if there exists a unitary map $R  : \F_F \rightarrow \F_F$ such that
%	\label{bogoliubov property 1}
%	R^* A(f,g) R
%	=
%	A( \nu (f,g) ) \ , \forall  \qquad f,g \in L^2 (\R^d ) \ .
%It is known that this condition is equivalent to $ v $ being Hilbert-Schmidt.
%\textit{Pure states}.
%Let us consider now the special case in which one is given
%an   orthogonal projection
%$\omega: L^2(\R^d) \rightarrow  L^2( \R^d)$
%with $\t{Tr} \omega=M$.
%Then, there exists an orthonormal set $ \{  f_i\}_{i=1}^M \subset L^2 (\R^d)$
%such that
%	\omega =
%	 \sum_{ i = 1 }^M  \ket{f_i} \bra{f_i}
%	  \ .
%We now define the operators $u_\omega, v_\omega : L^2(\R^d) \rightarrow L^2(\R^d)$ according to
% $	u_\omega  \equiv  1 - \omega $
% and
% $
%	v_\omega
%	\equiv
%	 \sum_{ i = 1 }^M  \ket{ \overline{ f_i}} \bra{f_i} .$
%One may easily check that
% $
%	u_\omega = u_\omega^* = u_\omega^2
%	$
%$	v_\omega^* v_\omega = \omega $ .
%Thus, it holds true that
%$u_\omega^* u_\omega + v_\omega^*v_\omega =1 $
% $u_\omega^* \bar v_\omega + v_\omega^* \bar u_\omega = 0 $
%and   the associated map $\nu_\omega$ defined in \eqref{nu}
% is a Bogoliubov transformation.
%We denote by
%	R[\omega ] : \F_F \rightarrow \F_F
% the associated implementator, and by $R^*[\omega]$  its adjoint map.
% In particular, the conjugation of the creation- and annihilation operators
% with $R[\omega] $ gives the relations
% \begin{align}
% 	R^*[\omega]  a^*_x R[\omega]
% 	& =
% 	a^* ( u_{\omega,x})
% 	+
% 	a( \overline{  v_{\omega,x} } )		\\
% 	R^*[\omega] a_x R[\omega]
% 	& =
% 	a ( u_{\omega,x})
% 	+
% 	a^*( \overline{  v_{\omega,x} } )
% \end{align}
% where we denote
% $u_{\omega, x } (x') = u_\omega(x' , x)$
% and
% $v_{\omega, x } (x') = v_\omega(x' , x)$.
% These relations combined with the CAR  imply  that
% \begin{equation}
% 	 \< R[\omega]  \Omega_F    , a_{x'}^* a_x  R[\omega]   \Omega_F    \>_{\F_F}
% 	  =
% 	    \<  \overline{ v_{\omega, x'}     }     ,  \overline{ v_{\omega, x }}   \>_{L^2 }
% 	      =  \omega (x,x')   \, , \qquad
% 	      \forall x, x' \in \R^d
% \end{equation}
%In other words, the one-particle reduced density matrix of
%  $R[\omega]  \Omega_F $ corresponds to $\omega. $

\subsection{Bosons}
\label{subsection bosons}
Similarly, throughout this subsection
we emply the following notation
to denote the bosonic Fock space $_B$
\begin{equation}
\mathscr{F}_B
\equiv
\bigoplus_{n=0}^\infty
\mathscr{F}_B^{(n)}
\quad
\t{where}
\quad
\F_B^{(0)} \equiv \C
\quad
\t{and}
\quad
\mathscr{F}_B^{(n)} \equiv 	 L^2_s( \R^{dn} ) , \ \forall n \geq 1
\end{equation}
where $	L^2_s( ^dn )  $
corresponds to the subspace of  symmetric functions.
$_B$   is a Hilbert space when endowed with the inner product
\begin{equation}
\< \Phi_1,  \Phi_2 \>_{\F_B}
\equiv
\sum_{n=0}^\infty
\<\Phi_1^{(n) } , \Phi_2^{(n)}    \>_{\F_B^{(n)}}  \  ,
\qquad \forall \Phi_1 ,\Phi_2 \in \F_B \ .
\end{equation}
On the bosonic Fock space $_B$ one introduces the smeared-out
creation- and annihilation operators
as follows.
Given $f ∈L^2(^d) $,
we let   $b^*(f)$ and $b(f)$
be defined   for  $Φ∈_F$ as
\begin{align}
\big(		 b^*(f)  \Phi \big)^{(n)}
(y_1, \ldots, y_n)
& \equiv
\frac{1}{\sqrt n }
\sum_{i=1}^n
f( y_i) \Phi^{(n-1)}
(y_1 , \ldots, y_{i-1} , y_{i+1} , \ldots, y_n )	\ , \\
\big(	 b (f)	 \Phi \big)^{ (n)}
(y_1, \ldots, y_n)
& \equiv
\sqrt{n+1 }
\int_{\R^d }  \overline{f(y)} \Phi^{(n+1)} (y , y_1 , \ldots, y_n) \d y  \ .
\end{align}
In contrast to the fermions, they satisfy the Canonical Commutation Relations (CCR)
\begin{equation}
[  b(f)  , b^*(g)  ]   = \< f,g\>_{L^2}
\quad
\t{and}
\quad
[   b^\hash (f) , b^\hash (g)	] = 0  \ .
\end{equation}
In particular, they are unbounded operators on $_B$, but relatively bounded with respect to the bosonic number operator--see Lemma \ref{lemma boson estimates}.
They are connected to the  operator-valued distributions $b_y^*$ and $b_y$  by means of the formulae
\begin{equation}
b^*(f) = \int_{\R^d }  f(y)  b_y^* \d y
\qquad
\t{and}
\qquad
b(f) = \int_{	\R^d	} \overline{f(y)}
b_y \d y
\  , \qquad f \in L^2(   \R^d )    \ .
\end{equation}

\subsubsection{Operator estimates}
We proceed analogously as we did for fermions.
Namely,  given an operator $Ø$ in the one-particle space $L^2(^d)$,
we  considers its second quantization  $Γ̣_B [Ø]$
acting on $_B$  defined as
\begin{equation}
( 	\d\Gamma_B[\O]  )^{ (n)}
\equiv
\sum_{ i = 1 }^n
\mathbb{1}^{ i -1 } \otimes \O \otimes \mathbb{1}^{ n - i }  ,
\ n \geq 1
\qquad
\t{and}
\qquad
( 	 	\d\Gamma[\O]  )^{ (0)} \equiv 0  \ .
\end{equation}
Similarly, if $Ø$ has an operator kernel
$Ø(y,y')$,    one may write in terms of creation- and annihilation operators
\begin{equation}
\d\Gamma_B [\O ] = \int_{	\R^{2d}	} \O (y,y')      b_{y'}^*b_y     \d y \d y' \ .
\end{equation}

Analogously as we did in the case of fermions,
we now   introduce
the corresponding   relations
for the  { bosonic number operator}
\begin{equation}
\label{boson identity N}
\N_B  =   \bigoplus_{n=0}^\infty  n
\d \Gamma_B [\1 ] = \int_{	\R^d	} b_y^*b_y \d y  \ .
\end{equation}

Let us now  collect some basic results concerning estimates
for the second quantization of operators in bosonic Fock space.
For reference, see
\cite[Lemma 3.1]{ChenLeeSchlein2011}.

\begin{lemma}[Estimates for bosonic operators]
\label{lemma boson estimates}
Let $\O : L^2(\R^d)\rightarrow L^2(\R^d)$ be a bounded operator,
and let $ f  \in L^2(\R^d)$.
Then, the following holds true

\vspace{1mm}

\noindent (1) For all $\Phi \in \calD (\N_B^{\frac{1}{2}})$
\begin{align}
\|	  b( f )  \Phi		\|_{\F_b}			&  \leq \|	  f 	\|_{L^2}		\|   \N_B^{\frac{1}{2}}  \Phi		\|_{	\F_b	}			\\
\|	  b^*( f )  \Phi		\|_{\F_b}			&  \leq \|	  f 	\|_{L^2}		\|    (  \N_B+1)^{\frac{1}{2}}  \Phi		\|_{	\F_b	}		 \ .
\end{align}

\noindent (2) For all $\Phi \in \calD(\N_B )$
\begin{align}
\|	 \d \Gamma_B (\O )	 \Phi	\|_{\F_b} \leq 	  \| \O		\|_{B(L^2)}		\| \N_B	 \Phi		\|_{	 \F_b	} \ .
\end{align}
\end{lemma}

\subsubsection{Coherent states}
\label{subsubsection CH}

Analogously as we did for fermions, we introduce a unitary transformation that we shall make use of in the rest of the article.
Namely, for $f ∈L^2(^d )$ we introduce the \textit{Weyl operator} as
\begin{equation}
\W[f] \equiv \exp \Big(     b (f) - b^*(f)    \Big) : \F_B \rightarrow \F_B \ .
\end{equation}
Note that since the argument in the exponential is anti self-adjoint, $[f]$ is automatically a unitary map.
Its action on the vacuum vector creates a \textit{coherent state}.
That is,    thanks to the Baker-Camper-Hausdorff formula, one has
on $_B $
\begin{equation}
( 	\W[f] \Omega_B)^{(n) }
e^{   - \frac{\|	f\|^2  }{2}     }
\frac{ f^{\otimes n }}{\sqrt{n!}}\  , \quad \forall n \geq  0 \ .
\end{equation}
In particular, the probability of finding the system with
$n$ particles is given by
$ e^ -  f  	 f^n / n!  $,
which follows a Poisson distribution with parameter $λ= 	f_L^2$.

\vspace{1mm}

We collect in the following lemma
some properties of  the  Weyl operator
and coherent states.
For more details, we refer the reader to
\cite[Lemma 2.2]{RodnianskiSchlein2009}

\begin{lemma}[Properties of $\W$]
\label{lemma properties bosons}
Let $f\in L^2(\R^d)$, and $\W[f]$ be as above.
Then, the following statements hold true
\begin{enumerate}
\item
$\W[f]$ is a unitary transformation on $\F_B$,
$\W^*[f] = \W[-f]$.

\item
For all $y \in \R^d$
\begin{align}
\label{conjugation f}
\W^*[f] b_y  \W [f]
& 	=
b_y +    f(y) 	\ , \\
\W^*[f] b^*_y  \W [f]
& 				=
b^*_y  +   \overline{ f(y) } \   .
\end{align}

\item
For all $y,y' \in \R^d$
\begin{equation}
\< \W[f] \Omega_B , b_{y'}^* b_y \W[f] \Omega_B \>_{\F_B}
\overline{ f(y')} f(y) \ .
\end{equation}
In   words, the one-particle reduced density matrix of $\W[f] \Omega_B $
corresponds to the projector $ \ket{f} \bra{f}.$

\end{enumerate}

\end{lemma}

%\section{From Hartree to Vlasov}
%In this section, we are interested in building a roadmap for taking the $ \hbar \downarrow  0 $ limit of the coupled Hartree equations, introduced in the previous section.
%We will choose the following parameters for the theory
% \hbar^d = 	N = M = 1 / \lambda
%	\qquad
%	\t{and}
%	\qquad
%	m_f = m_b = 1 \ .
%Namely, we consider (we make the $\hbar$ dependence explicit)
%set of equations given by
%	\label{coupled hartree}
%	\begin{cases}
%		& i \hbar \partial_t \omega^\hbar
%		=
%		[   - (\hbar^2 / 2 )\Delta
%		+
%		(v  *  \rho_b^\hbar )  ,
%		\omega^\hbar
%		]			\\
%		&   i \hbar \partial_t  \vp^\hbar
%		=
%		- (\hbar^2 / 2 )
%		\Delta  \vp^\hbar
%		+
%		(v  *  \rho_f^\hbar )
%		\vp^\hbar
%		  	\\
%		& (\omega^\hbar  ,\vp^\hbar )(0)  = (\omega_0^\hbar  , \vp_0^\hbar ) \in X
%	\end{cases}     \  .
%Here, we denote the (time-dependent) fermionic and bosonic densities by
%\rho_b^\hbar   (x)  =      |   \vp^\hbar (x)  |^2
%\rho_f^\hbar   (x)
%\omega^\hbar (x;x)     ;
%\qquad \t{ where } x \in \R^d \ .
%In particular, they are both normalized to unity:
%$ \int_\R^d \rho^\hbar_b (x) \d x = \int_\R^d \rho^\hbar_f (x) \d x =1  $ .
%We shall make additional assumptions on the initial data $(\omega_0^\hbar , \vp_0^\hbar )$ regarding their
%\textit{semiclassical structure}, and \textit{regularity}.
%\subsection{What happens to fermions?}
%In order to take the $ \hbar \downarrow  0 $ limit, we have to study the semiclassical limit of the solutions of \eqref{coupled hartree}--this is one way to do it.
%Let us recall that position space is $\R^d = \R / 2 \pi \Z$, and we denote its dual lattice by
%$\R^d = \Z^d $.
%We introduce the space for momentum as the $\hbar$-dependent lattice
%	\R^d_\hbar^* : = ( \hbar \Z)^d  \equiv  \hbar \R^d
%	\quad
%	\t{ together  with }
%	\quad
%	\int_{\R^d_\hbar }
%	F(p)
%	 \d p
%	: =
%	\hbar^3 \sum_{ k \in \Z^d } F( \hbar k)
%for any summable function $F : \R^d_\hbar^* \rightarrow \C $.
%	Let $\hbar>0$.
%	\item 	We introduce the \textbf{Wigner transform} as the   map
%	\begin{equation}
%		W^\hbar :
%		L^2_{x,y} ( \R^d \times \R^d )
%		\rightarrow
%		L^2_{x,p} (  \R^d \times \R^d_\hbar  )
%	\end{equation}
%	initially defined on smooth functions by
%	\begin{equation}
%		W^\hbar [ \gamma ] (x,p)
%		=
%		\frac{1}{ (2 \pi)^{ \frac{d}{2} } }
%		\int_{\R^d}
%		\gamma
%		\Big(		 x + \frac{y}{2 } ,   x - \frac{y}{2}			\Big)
%		e^{	 - i \frac{ p \cdot y  }{ \hbar } 	}
%		\d y  \ , \qquad (x,p ) \in \R^d \times \R^d_\hbar
%	\end{equation}
%	and then closed.
% 	\item  We introduce the \textbf{Weyl quantization}
% 	 as the   map
% 	 	\begin{equation}
% 	  \mathrm{Op}^\hbar  :
% 	 	L^2_{x,p} (  \R^d \times \R^d_\hbar  )
% 	 	L^2_{x,y} ( \R^d \times \R^d )
% 	 \end{equation}
% initially defiend on smooth functions by
%	  	  \mathrm{Op}^\hbar [a] (x,y)
%	  	   =
%	  	   \frac{  \hbar^{-d }  }{ (2 \pi )^{d/2}}
%	  	   \int_{\R^d_\hbar }
%	  	    a \Big(	 \frac{x  +  y}{2}  ,  p 		\Big)
%	  	    	e^{	  + i \frac{ p \cdot (x-y)  }{ \hbar } 	}
%	  	    	\d   p
%	  	    	\ , \qquad (x,y ) \in \R^d \times \R^d
%and then closed.
%	As is well-known,  the  Wigner transform and the Weyl quantization are inverses of one another:
%	W^\hbar \circ  	  \mathrm{Op}^\hbar =  \1_{ 	 L^2_{x,p}	}
%	\qquad
%	\t{and}
%	\qquad
%		\circ
%		W^\hbar
%		 =  \1_{ 	 L^2_{x,y}	} 	 ;
%a fact that we will heavily exploit.
%Consider  now the fermionic one-particle density
%$\omega^\hbar $
%solving the coupled Hartree equation \eqref{coupled hartree}.
%Abusing notation and denoting its kernel with the same symbol $\omega^\hbar$,
%we introduce its (normalized) Wigner transform as
%	f^\hbar (t) : =
%	\frac{1}{N }
%	W^\hbar [  \omega^\hbar (t) ]  \in L^2_{x,p} (\R^d \times \R^d_\hbar^* )  \  , \qquad \forall t \in \R \ .
%Since $\omega^\hbar$ is trace-class and $\mathrm{Tr} \omega^\hbar = N$, it readily follows that
%	\int_{\R^d\times \R^d_\hbar}
%	f^\hbar (t , x , p ) \d x \d p
%	= 1  \ , \qquad \forall t \in \R    .
%We remark that
%$ f^\hbar $
% may take negative values and, therefore, we cannot   interpret it as a probability distribution.
% Noting that formally $\R^d \times \R^d_\hbar^* \rightarrow  \mathbb T ^d  \times \R^d$ as $\hbar \downarrow 0 $,
% it still makes sense to ask that $\omega^\hbar_0$
% is a sequence of initial data for which there exists $f_0 \in L^1_{x,p} \big(  \mathbb T ^d  \times \R^d \big)  $
% such that
% \begin{equation}
% 	f_0^\hbar \equiv f^\hbar (0) \longrightarrow  f_0  \  ,
% 	\qquad
% 	\t{ as }
% 	\hbar \downarrow  0  \  .
% \end{equation}
%We then ask if we can \textit{propagate} for $t> 0 $,  meaning that
%	f^\hbar (t)  \longrightarrow f (t)   \qquad \t { as }
%	\hbar \downarrow 0
%for some $f(t) \in L^1_{x,p} ( \T^d \times \R^d )$ solving an appropiate equation.
%\subsection{What happens to bosons?}
%Let us describe the heuristic behing our  study of the semiclassical limit $ \hbar \downarrow 0 $    of the boson field.
%First, recall that $\vp^\hbar$
%represents the one-particle wave function of a Bose-Einstein condensate that is experiencing full condensation:
%the $N$ bosonic particles of a gas  are occupying this state.
%Experimentally, such states are constructed by introducing a  strong enough external potential, i.e. a trapping force.
%Assuming this potential is  \textit{harmonic}: $V_{ \t{ext} } (x) = \kappa x^2 $ for some $ \kappa >0 $,
%one expects the boson field will tend to minimize the the energy functional
%	 \< \vp  ,   [  ( -\hbar^2 /2  )  \Delta     + V_{\t{ext}}     ] \vp  		 \>_{L^2 } ;
%the harmonic trap  $V_{\t{ext}}$ would dominate every other interaction of the system.
%It turns out that    the minimizer of the above functional corresponds to the ground state
% of the harmonic oscillator
% \begin{equation}
%  H  = - \frac{\hbar^2 }{2 } \Delta_x  +  \kappa x^2     \t{  on  } L^2_x ( \R^d )
% \end{equation}
%which  is  well-known to be a  \textit{coherent  state}, i.e.
% a minimizer of the position and momentum uncertainties:
% \begin{equation}
%\Delta_\vp (x) : = \< \vp , (x -   \< x\>_\vp    )^2 \vp  \>^{\frac{1}{2}  }
%\t{ and }
%\Delta_\vp (p ) : = \< \vp , (  p  -   \<  p \>_\vp    )^2 \vp  \>^{\frac{1}{2}  }
% \end{equation}
%where $ p = - i \hbar \nabla_x$ is the standard momentum operator in $L^2 $.
%In other words, coherent states saturate Heisenberg's uncertainty principle $\Delta_\vp (x) \Delta_\vp (p) \geq \hbar $
%by balancing out their error in $\hbar$
%	\Delta_\vp  (x)  = \Delta_\vp (p ) = \hbar^{1/2} \ .
%We assume that our initial data  $\vp_0^\hbar$ is   produced by  means of an external, harmonic trap and, therefore,
%approximately minimizes Heisenberg's uncertainty principle.
%Mathematically, we encode this property in the  following assumption.
%\begin{condition}[Boson initial data]
%	The sequence $\vp_0^\hbar \in L^2_x ( \R^d)$
%	of normalized
%	initial data for the boson field satisfies the following two estimates.
%There exists $(x_0 , p_0) \in \T^d \times \R^d $	 and $c>0$
%such that  for all $ \hbar > 0 $
% there holds
%	\begin{align}
%\label{variance position}
% &
% 	 \<  \vp , | x -x_0|^2 \vp    \>_{L^2_x }^{1/2}
% 		 \leq c \hbar^{1/2}   \  , \\
%\label{variance momentum}
% &
% 	 \<  \vp , | p  -p _0|^2 \vp    \>_{L^2_x }^{1/2}
%     \leq c \hbar^{1/2} \
%	\end{align}
%where we recall $ p = - i \hbar \nabla_x$
%is the standard momentum observable in $L^2_x (\R^d)$.
% \begin{remark}
% 	The triangule inequality easily implies that
% 	\begin{equation}
% |  		 \< \vp_0^\hbar  , \, x  \, \vp_0^\hbar  \>_{L^2 } -  x_0 | \leq c \hbar^{\frac{1}{2}}
% \qquad
% \t{and}
% \qquad
%   \Delta_{\vp_0^\hbar} (x) \leq  c \hbar^{\frac{1}{2}}
% 	\end{equation}
%where $\Delta_\vp(x)$ is the uncertainty (standard deviation) of the observable $x$ at state $\vp$.
%Analogous estimates hold for momentum $p$.
% \end{remark}
%	A simple example that satisfies the above two estimates is given by re-scaling
%	an arbitrary smooth profile $ g \in C^\infty (\T^d )$ according to
%	\begin{equation}
%		\vp_0^\hbar  ( x)
%		=
%		\hbar^{-d/4}
%		e^{ - i x\cdot  p_0 / \hbar}
%		g (\hbar^{-1/2}  (x - x_0)) \qquad \forall x \in \T^d \ .
%	\end{equation}
%The next natural question is how to study convergence of $\vp^\hbar$ as $\downarrow 0 $.
%This is one way to do it.
%Notice that as $\hbar \downarrow 0 $,  the variance of $\vp_0^\hbar$ goes to zero and,
%is approximated by a point in phase-space $(x_0 , p_0) \in \T^d \times \R^d$.
%We would then like to study the trajectory of this point accross phase space, when $t>0$.
%This can be accomplished as follows.
%For $t \in \R$,
%we introduce the
%position and momentum measures
% of the boson field $\vp^\hbar (t)$ as
%  \begin{align}
%  	\mu^\hbar_x (t , \d x )  & : = |  \vp^\hbar (t,x)  |^2 \d x    \   ,  \\
%  	  	\mu^\hbar_p (t , \d p  )  & : = |  \widehat \vp^\hbar (t,p )  |^2 \d p
%  \end{align}
%where we have introduced the $\hbar$-dependent Fourier transform
%	\widehat{\vp} (p)
%	:  =
%	\frac{1}{ (2 \pi \hbar )^{\frac{d}{2}} }
%	\int_{\R^d}
%	\vp(x) e^{ - i  x\cdot p/ \hbar 	}
%	\d  x   \  ,
%	\qquad
%	\forall p \in \R^d_\hbar^* \ .
%It turns out that   estimates of the form
% \eqref{variance position} and \eqref{variance momentum}
% --that is, variance estimates--
%can be understood as convergence
%of the measures $\mu^\hbar_x$ and $\mu^\hbar_p $ to Dirac measures centered in $(x_0 , p_0 )$,
%in a certain metric which we now explain.
%\begin{definition}[Wasserstein distance]
%Let  $ (E,d) $ be a metric space and $ m \in [1, \infty )$.
%	On the space  $\mathcal{P}_m(E)$
%	of   probability measures  over $E$ with $m$ finite moments,
%	we define the $m$-\textit{Wasserstein distance}
%	\begin{align}
%		\d_m  [\mu , \nu ]
%		: =
%		\inf_{	 \gamma \in \calC (\mu , \nu )	} \int_E  d (x,y)^m
%		\d \gamma(x,y)
%	\end{align}
%	where $\calC (\mu , \nu)$
%	is the space of all couplings of $\mu$ and $\nu$, i.e.
%	a probability measure on  $E \times E$
%	with first marginal  $\pi_1 \gamma = \mu $ and  second marginal $\pi_2 \gamma = \nu $.
%	In our case, we will consider $E = \T^d$ for position, and $E = \R^d$ for momentum; we do not make explicit the dependence of $\d_m$ on the underlying metric space $E$.
%	Consider $E  = \R$.
%	In general, if $ \nu = \delta_{x_0}$ is a Dirac measure,
%	one may show that there is only one coupling
% between $\mu$ and $\nu$.
% Consequently,
% their $m$-th Wasserstein distance takes the form
% \begin{equation}
% 	\d_m (\mu , \nu) = \int_\R | x - x_0|^m \d \mu (x) \ .
% \end{equation}
%The above formula clearly holds for general spaces.
%This then implies that for the position/momentum measures of the boson field, there holds at $t = 0 $
%	 \d_2
%	    [    \mu_x^\hbar  (0) , \delta_{x_0}   ]
%	  & =
%	  \int_{\T^d}
%	  |x - x_0|^2
%	  | \vp_0^\hbar (x) |^2 \d x
%  \leq c \hbar  \\
%	 \d_2  [    \mu_p^\hbar  (0) , \delta_{p_0}   ]
%	 & =
%	 \int_{\R^d}
%	 | p  - p_0|^2
%	 | \widehat \vp_0^\hbar ( p ) |^2 \d  p
%	   \leq c \hbar
%for some updated constant $c>0 \ . $
%One then hopes that for $ t >0$
% 	 \mu_x^\hbar (t)  (\d x )
%  &	 \longrightarrow \delta  ( x - x(t) )    \d x
%	 \qquad
%	 \t{as}
%	 \qquad
%	 \hbar \downarrow 0  \\
%	  	 \mu_p^\hbar (t)  (\d  p )
%	 &	 \longrightarrow \delta  ( p - \dot x(t) )    \d p
%	 \qquad
%	 \t{as}
%	 \qquad
%	 \hbar \downarrow 0
%for some trajectory $x(t) \in \T^d $, with convergence being measured in 2-Wasserstein distance.
%\subsection{Vlasov-Newton equation}
%Essentially, we have motivated the following question.
%Can we prove convergence from the one-particle reduced densities
%$(\omega^\hbar  , \vp^\hbar)$
%to some semi-classical object $ (f , x ) \in  L^\infty_t  \big( L^1_{x,p}  \times \T^d\big) $, as $\hbar \downarrow 0$?
%The most natural guess is yes, provided this object satisfies the following coupled \textit{Vlasov-Newton} equation
%	 & 	\big( 	 \partial_t + p \cdot \nabla_x + F_b ( t, x) \cdot \nabla_p	 \big)  f    =  0  \\
% & 	\ddot x (t)    =  F_f \big( t, x  (t)  \big)	\\
% & 	\big(	 f(0) , x (0) , \dot x(0)	\big)    =   \big(  f_0 , x_0 , p_0  \big)  \in L_{x,p}^1 \times \T^d \times \R^d
%where $F_f$ and $F_b$
%correspond to the forces arising from fermions, and bosons, respectively; they are given by
%  =
%-  \nabla v \big(x -  x(t) \big)
%\t{ and }
%F_f (t,x)	   =  - \int_{  \T^d \times \R^d  }   \nabla v  (x - y ) f ( t ; y,p ) \d  y \d p  \ .
%Note that these two forces
%are of the form $F (t)  = - \nabla (v * \rho (t) )$
%where $\rho$
%is either the fermionic, or the bosonic \textit{position} distributions:
%	\rho_b (t, x) = \delta \big( x- x(t) \big)
%	\qquad
%	\t{or}
%	\qquad
%	\rho_f(t , x) = \int_{\R^d  }
%	f(t,x ,p ) \ d p  \ .
%We hope we can prove the following result -- it is nothing but the analogous version of Theorem 2.5 contained in (BPSS2016).
%(I am more or less convinced we should change $ \R^d  = \T^d$ for $\R^d$; the following result would contain that modification)
%\begin{theorem}[From Hartree to Vlasov]
%	\label{theorem convergence 2}
%	Assume the interaction potential  $v \in L^1 (\R^d , \R   ) $  is  such that
%	\begin{equation}
%		 \|	 v	\|_{  W^{2,1} } : =  \int_{\R^d }  | \hat v (p) | \big(  1 + | p |^2 \big) \d p < \infty \ .
%	\end{equation}
%Let $ ( \omega_0^\hbar , \vp_0^\hbar )_{\hbar  \in (0,1)} $
%be a sequence of inital data for the coupled Hartree equation
%such that  the following   conditions are satisfied.
%	\item (Fermions)
%For all $\hbar \in (0,1)$,  	$\omega_0^\hbar$ is a trace-class operator on $L^2(\R^3 )$
%	such that $ 0 \leq  \omega_0^\hbar \leq  1$ and
%	$\mathrm{Tr} \omega_0^\hbar = N.$
%	Additionally, there exists $C>0$
%such that for all $ \hbar \in (0,1)$
%the following semi-classical structure is satisfied
%	\begin{equation}
%		\mathrm{Tr}
%		 \
%		\big|
%		\,
%		 [x , \omega_0^\hbar ]
%		 \,
%		\big|
%		\leq  C N \hbar
%		\quad
%		\t{ and }
%		\quad
%		\mathrm{Tr} \
%		\big| \,    [ p  , \omega_0^\hbar ] \,
%		 \big|
%		\leq  C N \hbar
%	\end{equation}
%where $ p = - i  \hbar \nabla_x $.
%Furthermore, denoting by  $f_0^\hbar  = W^\hbar [ \omega_0^\hbar  ]$
%the Wigner transform of $\omega_0^\hbar$,
%we assume that
%	\item $ \sup_{\hbar \in  (0,1)} \|	 f_0^\hbar	\|_{W^{1,1}} < \infty$.
%	\item There exists a probability density $f_0 \in W^{1,1}_{x,p}(\R^d \times \R^d  )$
%	and a sequence $( \delta_{\hbar })_{ \hbar \in (0,1)}$
%	such that
%	\begin{eqnarray*}
%		\|	 f_0^\hbar - f_0	\|_{L^1  } \leq \delta_{\hbar } \  , \qquad \forall \hbar \in  ( 0,1 ) \ .
%	\end{eqnarray*}
%\item (Bosons) For all $ \hbar \in (0,1)$, $\vp_0^\hbar \in L^2(\R^d )$ is such that $\|  \vp_0^\hbar 	\|_{L^2} =1 $.
%Additionally, there exists $(x_0  ,  p_0) \in \R^d \times \R^d $  and
%$C>0$ such that for all $ \hbar> 0 $
%the following semi-classical structure is satisfied
%	\<  \vp , | x -x_0|^2 \vp    \>_{L^2_x }^{  \frac{1}{2 } }
%	\leq C \hbar^{  \frac{1}{2}}
%	\quad
%	 \t{ and }
%	 \quad
%	\<  \vp , | p  -p _0|^2 \vp    \>_{L^2_x }^{ \frac{1}{2}}
%	\leq C \hbar^{ \frac{1}{2}}
%where $ p = - i  \hbar \nabla_x $.
%Let $(\omega^\hbar , \vp^\hbar)$
% be the solution of the coupled Hartree equation
%	\begin{cases}
%		& i \hbar \partial_t \omega^\hbar
%		=
%		[   - (\hbar^2 / 2 )\Delta
%		+
%		(v  *  \rho_b^\hbar )  ,
%		\omega^\hbar
%		]			\\
%		&   i \hbar \partial_t  \vp^\hbar
%		=
%		- (\hbar^2 / 2 )
%		\Delta  \vp^\hbar
%		+
%		(v  *  \rho_f^\hbar )
%		\vp^\hbar
%		\\
%		& (\omega^\hbar  ,\vp^\hbar )(0)  = (\omega_0^\hbar  , \vp_0^\hbar )
%	\end{cases}
%$ \rho_b^\hbar   (x)  =      |   \vp^\hbar (x)  |^2 $
%\rho_f^\hbar   (x)
%N^{-1 }
%\omega^\hbar (x;x)    ;
%and denote by
%	f^\hbar : = \frac{1}{N } W^\hbar [\omega^\hbar ]  \ ,
%			\mu^\hbar_x (\d x )  : = |  \vp^\hbar (x)  |^2 \d x
%			 \
%			\t{  and  }
%			 \
%				\mu^\hbar_p ( \d p  )    : = |  \widehat \vp^\hbar ( p )  |^2 \d p
%the Wigner transform of the fermionic one-particle density,
%and the  bosonic position and momentum distributions, respectively.
%Let  $ (f,x)  \in
% L^\infty_t L^1_{x,p}(\R_+ \times \R^d \times \R^d  )
% \times
% C^2(\R^d )
%be the solution of the coupled Vlasov-Newton equation
%	\begin{cases}
%		& 	\big( 	 \partial_t + p \cdot \nabla_x + F_b ( t, x) \cdot \nabla_p	 \big)  f  (t,x,p)   =  0  \\
%		& 	\ddot x (t)    =  F_f \big( t, x  (t)  \big)	\\
%		& 	\big(	 f(0) , x (0) , \dot x(0)	\big)    =   \big(  f_0 , x_0 , p_0  \big)  \in L_{x,p}^1 \times \R^d \times \R^d
%	\end{cases}
%$F_b  (t,x )   =
%	-  \nabla v \big(x -  x(t) \big)
% $
% and
% $
% 	F_f (t,x)	   =  - \int_{  \R^d \times \R^d  }   \nabla v  (x - y ) f ( t , y,p ) \d  y \d p  \ .
% $
% Then, there exists $C>0$
% such that for all $\hbar>0$  and $t \in \R $
% the following inequalities hold true
% \begin{align}
% 	 \|    f (t)  		  -  f^\hbar (t)	\|_{ W^{  - 2 , \infty} }  & 		 \leq  C  (  \hbar^\frac{1}{2}  + \delta_{\hbar }) \exp (C |t|)	 \   , \\
% 	 \d_1  \Big(  \mu_x^\hbar (t) ,  \delta_{x (t) }		\Big)   & 		 \leq  C  (  \hbar^\frac{1}{2}  + \delta_{\hbar }) \exp (C |t|) \  , 	\\
% 	 \d_1  \Big(  \mu_p^\hbar (t) ,  \delta_{ \dot x (t) }		\Big)   & 		 \leq  C  (  \hbar^\frac{1}{2}  + \delta_{\hbar }) \exp (C |t|)	  \   ,
% \end{align}
%where $\d_1$ is the 1-Wasserstein metric, and $\|    u  \|_{W^{ - k ,\infty}} : = \sup_{ ( \xi  , \eta ) \in \R^{2d}}    (1 + |\xi| + |\eta|)^{-2 }    | \hat   u (\xi, \eta )|$ .

\section{Second Quantization II: The Fluctuation Dynamics}
\label{section fluctuation}
The main goal of this section is to set up the proof of Theorem \ref{thm1}.
Namely, we will introduce and study the fluctuation dynamics around a state consisting of a   degenerate Fermi gas, and a Bose-Einstein condensate, evolving in time according to the mean-field equations \eqref{hartree}.
We prove that the number of fluctuations around this state is small, relative to the numbers $N$ and $M$.
In the next section, we show that these estimates imply Theorem \ref{thm1}.

\vspace{2mm}

This approach is nowadays considered standard,
and has been successfully employed in the derivation of several mean-field equations from many-particle systems.
The first work to use these techniques in the derivation of the Hartree equation
for bosons is [55],
while for fermions is [9].
Our proofs are heavily inspired by their ideas,
and actually borrows a few estimates.
Finally, let us also   refer the reader to the book [11]
for a cohesive treatment of the subject.

\vspace{2mm}

\textit{A note on domains}.
In what follows, we will be extensively
manipulating the number operators $_F$ and $_B$.
These are positive, unbounded self-adjoint
operators   on the Hilbert space $ℱ$,
with domains
$(_F) = { Ψ∈: ∑_n,m n^2  Ψ^(n,m) ^2_ _F^(n) ⊗_B^(m)      < ∞}$
$(_B) = { Ψ∈: ∑_n,m m^2  Ψ^(n,m) ^2_ _F^(n) ⊗_B^(m)      < ∞}$.
In order to simplify the exposition,
we shall
avoid making reference to the unbounded nature of these operators.
Let us note that there is no risk in doing so.
Indeed, in applications
all the states $Ψ∈ℱ$
that we manipulate
belong to the intersection $∩_k=1 ^∞(_F^k) ∩(_B^k)$,
and  all the dynamics   we consider
leave the above space invariant.

\subsection{Number Estimates}
Throughout this section, we denote by   $(ω, )$   the  pair of variables that solves
the Hartree-Hartree equations
\eqref{hartree}.

\vspace{2mm}

Let us   now introduce a fundamental  family of  unitary transformations in our analysis.
Namely, using the notation from Section \ref{section second quantization},
we define the following time-dependent
particle-hole transformation, and Weyl operator
\begin{align}
\label{definition W and R}
\mathcal{R}_t
\equiv
\mathcal{R} [\omega (t) ] : \F_F \rightarrow \F_F
\quad
\t{and}
\quad
\mathcal{W}_t \equiv
\mathcal{W} [\sqrt N \vp(t) ] : \F_B \rightarrow \F_B \ .
\end{align}
Note that both  of these transformations will   map
the respective vaccuum vectors
into states in Fock space whose one-particle reduced density matrices
correspond to  $ω(t)$ and $  N |(t)⟩ ⟨(t)|$, respectively.
In other words, for all $t ∈$ and $(x,y)∈^2d$
\begin{align}
\<
\mathcal{R}_t  \Omega_F , a_y^* a_x
\mathcal{R}_t \Omega_F
\>_{\F_F}
\omega(t ; x,y)
\quad
\t{and}
\quad
\<
\mathcal{W}_t  \Omega_b , b_y^* b_x
\mathcal{W}_t \Omega_b
\>_{\F_B}  =
\overline{  \vp(t, y ) }  \vp(t,x) .
\end{align}

\vspace{2mm}

We proceed to define the two-parameter
family of unitary transformations
$(t,s) : →$, which we refer to as  the \textit{fluctuation dynamics},
as follows
\begin{equation}
\label{unitary U}
\U(t,s)
\equiv
\big(
\mathcal{R}_t \otimes \mathcal{W}_t
\big)^*
\exp \bigg[ - i \frac{(t-s)}{\hbar} \H   \bigg]
\big(
\mathcal{R}_s \otimes \mathcal{W}_s
\big)  \ , \qquad t,s \in\R \ .
\end{equation}
The fluctuation dynamics measures,
in a quantitative way,
how far the many-body Schr\"odinger dynamics
is from the mean-field variables $(ω(t),(t))$.
In order to make this statement precise,
we 	recall that we have defined on $_F$ and $_B$ the
fermionic and bosonic number operators, respectively,
\begin{align}
\mathcal{N}_F    =
\int_{\R^d } a_x^* a_x  \d x
\qquad
\t{and}
\qquad
\mathcal{N}_b
\int_{\R^d }
b_y^* b_y \d y \  .
\end{align}
Unless confusion arises, we denote with the same symbols their natural extension to
$$. Finally, we introduce   on $$  the \textit{total number operator}
\begin{equation}
\N \equiv \N_F   +  \N_B \ .
\end{equation}
%An additional important observable that we will use in our analysis is the
% \textit{charge operator}, defined as
%	\label{charge}
%	\calQ (t)
%	 \equiv
%	  \d \Gamma
%	\big[
%	\1 - \omega(t)
%	\big] \otimes \1
%	-
%	\d \Gamma \big[
%	\omega(t)
%	\big]	  	       \otimes \1
%	\ , \qquad \forall t \in \R \ .
%As we shall see, in our present situation the states that we look at have zero initial charge.
%That is,
%we   consider states $\Psi \in \F $ such that  $\Q(0) \Psi = 0$.
%Furthermore, the charge operator is a conserved quantity, meaning that $\Q(t) \U(t,s) = \U(t,s) \Q(0)$--we prove this in Lemma \ref{lemma charge}.
%Hence, $\Q(t) \U(t,0)\Psi = 0 $ for all such states.

\vspace{2mm}

The main result of this section is the following theorem.
It contains estimates for the growth-in-time
for the expectations of $_F$ and $_B$
with respect to the fluctuation dynamics $(t,s)$.

\begin{theorem}[Number estimates]
\label{prop number estimates}
Let $(\omega,\vp)$ satisfy \eqref{hartree}, and assume that the assumptions in Theorem \ref{thm1} hold true.
Let $\U(t,s)$ be the fluctuation dynamics.
Then, the following statements hold true
\vspace{1mm}

\begin{enumerate}[leftmargin=.7cm]
\item
For all $\ell,  k \in \mathbb N$
there is a constant $C>0$ such that
for all $\Psi \in \F $ and $t,s\in \R $
\begin{align}
\langle      \Psi,  &  \U^*(t,s)    \,   \calN_F^\ell   	 \, \U(t,s)  \Psi	\rangle_\F  \\
\nonumber
\leq
K(t-s )
\Big[
\|	 \Psi	 \|_\F
\| (\N^\ell + \1 )	 \Psi	\|_\F
\Theta_{k,\ell}
\|  (\1 + M^{-1} \N_F)^\ell    \Psi		\|_\F
\|	 	  ( \N +\1)^{k + 3/2 	}    \Psi 		\|_\F
\Big]   \  .
\end{align}

\vspace{1mm}

\item
For all $\ell,  k \in \mathbb N$
there is a constant $C>0$ such that
for all $\Psi \in \F $ and $t,s\in \R $
\begin{align}
\langle   \Psi,  & \U^*(t,s)     \,   \N_B^\ell   \, 	\U(t,s)  \Psi	\rangle_\F 	\\
\nonumber
\leq
K(t-s )
\Big[
\|	 \Psi	 \|_\F
\| (\N^\ell+ \1 )	 \Psi	\|_\F
\Theta_{k,\ell}
\|  (\1 + N^{-1} \N_B)^\ell   \Psi		\|_\F 	\|	 	  ( \N+\1)^{k + 3/2  }    \Psi 		\|_\F
\Big]  \ .
\end{align}

\end{enumerate}
Here, we denote
$K(t)  \equiv  C 	{  \exp     [  C   \lambda  \sqrt{NM/\hbar}\,
( 1 + \sqrt{  \hbar M / N}    )	\exp |t|       ]		} $
together with
\begin{equation}
\Theta_{k,\ell}
\equiv
\frac{\lambda \sqrt N }{\hbar}
\frac{M^\ell}{(\hbar M)^k} \geq 0 \ .
\end{equation}
\end{theorem}

\begin{remark}[Evolution of  the vacuum vector]
\label{remark control omega}
In applications,
we will consider $\Psi = \Omega$.
In particular,  there holds  $\| ( \N_F + \1)^\ell\Omega	\|_\F =1 $
for all $\ell \in \mathbb N$.
Thus, given $\ell \in \mathbb N$, we may take $ k = k_\ell$  as in the statement of Theorem \ref{thm1}.
Hence,  there holds
$\Theta_{k,\ell }\leq 1$
uniformly in the physical parameters $\lambda$, $\hbar$, $N$ and $M$.
In this situation one obtains the following estimate,
provided one re-updates the constant $C> 0 $
\begin{equation}
\|	 (  \N+\1)^\ell  \U(t,s)  \Omega  	\|_{\F}
\leq
K(t - s)  \ , \qquad \forall t,s\in \R \ .
\end{equation}
\end{remark}

\begin{remark}[Boundedness of operators]
\label{remark boundedness of N}
An important consequence of      Theorem \ref{prop number estimates}
is that for all $\ell$,
there exists $K_\ell \geq 1$ such that
for all $\Psi \in \F$
\begin{equation}
\|	 (  \N+1)^\ell  \U(t,s)	 \Psi	\|_\F
\leq
K(t-s )
\|   (\N + \1 )^{K_\ell }	 \Psi	\|_{\F} \ , \qquad \forall t,s\in \R  \ .
\end{equation}
(\N+\1)^\ell \U(t,s) (\N+\1)^{-K_\ell}
is a bounded linear operator in $\F$,
and the same holds for its adjoint.
We record this in the following statement
\begin{equation}
\|		(\N+\1)^\ell \U(t,s) (\N+\1)^{-K_\ell}	\|_{B( \F)}
\|		(\N+\1)^{-K_\ell } \U(t,s) (\N+\1)^{\ell }	\|_{B( \F)}
\leq K(t-s ) \  .
\end{equation}
Here of course, we have used that
$\U^*(t,s) = \U(s,t)$
and the symmetry
$K(t-s) = K(s-t). $

\end{remark}

We dedicate the rest of this section to the proof of the above Theorem.

\subsection{The infinitesimal generator}
In order to establish the number estimates contained in Theorem \ref{prop number estimates}
we need to study the time evolution of the fluctuation dynamics.
To this end, we introduce
infinitesimal generator $ℒ(t) $ of  $(t,s) $
as  the time-dependent,
self-adjoint operator on $$
determined by the equation
\begin{equation} \label{L operator}
\begin{cases}
& i \hbar \partial_t \U (t,s )
= \mathcal L (t)
\U (t,s ) \ , 	\\
& \U (t, t ) = \1 \ .
\end{cases}
\end{equation}
The computation of $Ł(t)$ is
tedious, but can be carried out explicitly.
Let us record the result of the  calculation in Lemma \ref{lemma infinitesimal generator} below, and postpone the proof to an Appendix.

\vspace{2mm}

\textit{Notation.}
Before we state the explicit form of the infinitesimal generator, we   introduce
useful notations.
Recall that $(ω(t), (t))$
is a solution of the Hartree-Hartree equation
\eqref{hartree}.
\begin{enumerate}[leftmargin=0.7cm, label=$\square$]
\item
We denote by
and $h_B(t)$
the following  time-dependent
one-particle  Hamiltonians on $H^2(\R^d)$
\begin{align}
\label{hamiltonians}
& 	\equiv
- \frac{\hbar^2 }{2m_F}
\Delta
\lambda  N
\big(V * \rho_B (t)  \big)    \ ,   \\
& 	\equiv
- \frac{\hbar^2 }{2m_B}
\Delta
\lambda  M   \big(V * \rho_F (t)	\big)    \ .
\end{align}
$\rho_F (t,x) = M^{-1 } \omega(t ; x,x)$
$\rho_B(t,x) = |\vp(t,x)|^2$
are the corresponding fermionic and bosonic position densities.

\item
Upon decomposing $\omega(t) = \sum_{ i = 1 }^M \ket{\phi_i (t)} \bra{\phi_i(t)}$ we denote
with the same notation as in Section \ref{section second quantization},
for all $ t \in \R$
\begin{equation}
\label{u and v}
\textstyle
u (t ) \equiv  \1 - \omega(t)
\qquad
\t{ and }
\qquad
v(t) \equiv \sum_{ i = 1 }^M
\ket{  \overline{ \phi_i (t) } } \bra{ \phi_i(t)  }
\end{equation}

\item
For fixed
$x \in \R^d $ and $ t \in \R$,
we denote by $u_{x,t}$ and $v_{x,t}$
the distributions
given by
\begin{equation}
u_{t,x} ( x ' ) = u( t;  x, x')
\qquad
\t{and}
\qquad
v_{t,x} (x ') = v(t ; x,x')
\end{equation}
for all $x' \in \R^d . $
\end{enumerate}

\begin{lemma}
\label{lemma infinitesimal generator}
Let $\U(t,s)$ be the unitary transformation defined in \eqref{unitary U}, and let $\L(t)$ be its inifinitesimal generator.
Then, $\L(t)$ admits the following representation (modulo scalars)
\begin{equation}
\L(t)
\d \Gamma_F [h_F(t) ] \otimes \1
\1 \otimes \d \Gamma_B [h_B(t)]
+ \lambda \sqrt N \L_{2,1}(t) + \lambda \L_{2,2}(t)
\qquad
\forall t \in \R \ .
\end{equation}
Here, $h_F(t)$ and $h_B(t)$ are the one-particle Hamiltonians defined in Eq. \eqref{hamiltonians}.
The  time-dependent
operators  $\L_{2,1} (t) $ and $\L_{2,2} (t) $
are self-adjoint operators on $\F$,
are given by the expressions (here, we supress time labels for convenience)
\begin{align}
\nonumber
\L_{2,1}
& =
\int_{ \R^d \times \R^d  } v(x - y )
a^*(u_x) a(u_x) \otimes \big(
\vp (y) b_y^* +  \overline \vp (y) b_y
\big)
\ \d x \d y 	\\
\nonumber
& -
\int_{ \R^d \times \R^d  } v(x - y ) a^*( \overline{v _x}) a(\overline{v_x})
\otimes
\big(
\vp (y) b_y^* +  \overline \vp (y) b_y
\big) \ \d x \d y 	\\
\nonumber
& +
\int_{ \R^d \times \R^d  } v(x - y ) a^*(  {u _x}) a^*(\overline{v_x}) 	 \otimes
\big(
\vp (y) b_y^* +  \overline \vp (y) b_y
\big) \ \d x \d y 	\\
& +
\int_{ \R^d \times \R^d  } v(x - y ) a(  \overline{v _x}    ) a( u_x  ) 	 \otimes
\big(
\vp (y) b_y^* +  \overline \vp (y) b_y
\big) \ \d x \d y  \ ,
\end{align}
\begin{align}
\nonumber
\L_{2,2}
& =
\int_{ \R^d \times \R^d  } v(x - y ) a^*(u_x) a(u_x)
\otimes b_y^*b_y
\ \d x \d y 	\\
\nonumber
& -
\int_{ \R^d \times \R^d  } v(x - y ) a^*( \overline{v _x}) a(\overline{v_x})
\otimes b_y^*b_y
\ \d x \d y 	\\
\nonumber
& +
\int_{ \R^d \times \R^d  } v(x - y ) a^*(  {u _x}) a^*(\overline{v_x})  \otimes b_y^*b_y
\ \d x \d y 	\\
& +
\int_{ \R^d \times \R^d  } v(x - y ) a(  \overline{v _x}    ) a( u_x  ) \otimes b_y^*b_y
\ \d x \d y \ .
\end{align}
\end{lemma}

Unfortunately--as it happens often with   similar mean-field theories--the above generator
is not exactly to prove estimates as the ones contained in Proposition \ref{prop number estimates}.
We shall introduce a modified generator,
together with a new \textit{truncated} fluctuation dynamics,
and then prove that these two are close together.

\subsection{Truncated dynamics}
We introduce $(t,s)$
as the strongly continuous, two-parameter family of  unitary operators that solves
\begin{equation}
\begin{cases}
& i \hbar \partial_t \widetilde \U (t,s )
=   \widetilde \L (t)
\widetilde 	\U (t,s ) \ , 	\\
& \widetilde \U (s, s  ) = \1 \  ,
\end{cases}
\end{equation}
where the \textit{truncated infinitesimal generator}  $Ł(t)$
is defined as  follows
\begin{equation}
\label{truncated generator}
\widetilde 	\L(t)
\d \Gamma_F [h_F(t) ] \otimes \1
\1 \otimes \d \Gamma_B [h_B(t)]
+ \lambda \sqrt N \widetilde \L_{2,1}(t) + \lambda \L_{2,2}(t)
\qquad
\forall t \in \R \    .
\end{equation}
The first, second, and fourth terms in \eqref{truncated generator}
are identical to those found in Lemma \ref{lemma infinitesimal generator}
for $Ł(t)$.
However,  for the third term   we have
introduced a cut-off in one of the off-diagonal terms originally found in $Ł_2,1(t)$:
\begin{align}
\nonumber
\widetilde	\L_{2,1}
& =
\1 \otimes \chi(\N_B \leq \hbar M ) 		\int_{ \R^d \times \R^d  } v(x - y )
a^*(u_x) a(u_x) \otimes \big(
\vp (y) b_y^* +  \overline \vp (y) b_y
\big)
\ \d x \d y 	\\
\nonumber
& -
\1 \otimes  		\chi(\N_B \leq \hbar M )
\int_{ \R^d \times \R^d  } v(x - y ) a^*( \overline{v _x}) a(\overline{v_x})
\otimes
\big(
\vp (y) b_y^* +  \overline \vp (y) b_y
\big) \ \d x \d y 	\\
\nonumber
& +
\int_{ \R^d \times \R^d  } v(x - y ) a^*(  {u _x}) a^*(\overline{v_x}) 	 \otimes
\big(
\vp (y) b_y^* +  \overline \vp (y) b_y
\big) \ \d x \d y 	\\
& +
\int_{ \R^d \times \R^d  } v(x - y ) a(  \overline{v _x}    ) a( u_x  ) 	 \otimes
\big(
\vp (y) b_y^* +  \overline \vp (y) b_y
\big) \ \d x \d y  \ .
\end{align}
The cut-off $_B ≤ħM $  is specifically tailored for the problem at hand,
and is introduced so that one can ``close the inequalities''
when   running  a Gr\"onwall-type argument.

\vspace{1mm}

Let us recall that we have introduced the number operators
$_F$ and $_B$.
In this subsection,
we shall give estimates for the growth-in-time
of expectations of these observables with
respect to the truncated fluctuation dynamics $(t,s).$
The main result of this subsection is the following proposition, containing
relevant number estimates.
\begin{proposition}
[Number estimates for the truncated dynamics]
\label{prop number estimates truncated}
Assume that the solution of the mean-field equations $(\omega,\vp)$
satisfies the following apriori bound
\begin{eqnarray}
\|       [e^{i \xi \cdot x}   ,\omega(t)  ]       \|_{\tr}
\leq C \exp(C t)  M \hbar \<\xi \> \  \ , \qquad \forall \xi \in \R^d  \ .
\end{eqnarray}
Then, for every $k \in \mathbb N $  there exists a constant
$C>0$ such that for all   $t,s \in \R$
$\Psi \in\F$
there holds
\begin{align}
\langle      \Psi,  \widetilde \U^*(t,s)  \N^k  	  \widetilde 	\U(t,s)  \Psi	\rangle_\F
\nonumber
& 		\leq
{  \exp     \bigg[  C   \lambda  \sqrt \frac{NM}{\hbar }
\bigg( 1 + \sqrt{\frac{\hbar M}{  N}}    \bigg)	\exp |t-s|      \bigg]		}
\<  \Psi ,  (\calN^k  + \1 ) \Psi	\>_\F \  .
\end{align}
\end{proposition}

The proof of the above estimates is based on the Gr\"onwall inequality, and the commutator estimates contained in the following lemma.

\begin{lemma}[Commutator estimates]
\label{lemma commutator estimates}
Under the same assumptions of Proposition \ref{prop number estimates truncated},
for all $\theta \in \R $ there is $C>0$ such that for all $t,s\in\R$   for all $\Psi , \Phi \in \F$
there holds
\begin{align}
|	 \langle	\Psi , [  \widetilde \L (t) , \calN_F	] \Phi  	\rangle_\F		|
& \leq
C e^{Ct }
\lambda
\Big(  \sqrt{NM \hbar }   + M \hbar \Big)
\bigg(
\|	 (\N + \1 )^{	 \frac{1 + \theta}{2}	}	 \Psi	\|_{\F}^2
%\< \Psi ,  (\N  +1  )^{1+\theta } \Psi \>_\F
\|	 (\N + \1 )^{	 \frac{1 + \theta}{2}	}	 \Phi	\|_{\F}^2
\bigg)  \ ,
\label{prop estimate fermion}
\\
|	 \langle	\Psi , [  \widetilde \L (t) , \calN_B	] \Psi 	\rangle_\F		|
& \leq
C e^{Ct }
\lambda
\Big(  \sqrt{NM \hbar }   + M \hbar \Big)
\bigg(
\|	 (\N + \1 )^{	 \frac{1 + \theta}{2}	}	 \Psi	\|_{\F}^2
%\< \Psi ,  (\N  +1  )^{1+\theta } \Psi \>_\F
\|	 (\N + \1 )^{	 \frac{1 + \theta}{2}	}	 \Phi	\|_{\F}^2
\bigg)  \ .
\label{prop estimate boson}
\end{align}
\end{lemma}

%In the proof, we shall be using the following well-known estimates for operators in second quantization.
%\begin{lemma}[Fermion Estimates]
%	\label{lemma fermion estimates}
%	Let $\O: L^2(\R^d) \rightarrow L^2(\R^d)$ 	 be a bounded operator.
%	Then, the following holds true.
%	\vspace{1mm}
%	\noindent (1) For all $\Psi , \Phi \in \calD(\N_F)$
%	\begin{align}
%		\|	 \d \Gamma (    \O	  	) \Psi	 \|_{\F_f}			&
%			 \leq	\|	   \O  	\|_{B(  L^2 )}
%		\|	  \N_F  \Psi 	\|_{\F_f}		 \ , \\
%| 		\< \Psi , \d  \Gamma  (	 \O 	)\Phi \>_{\F_f}	|
%	 \leq	\|	   \O  	\|_{B(  L^2 )}
% \| 			\N_F^{\frac{1}{2}}	 		\Psi 	\|_{\F_f}
%  \| 				\N_F^{\frac{1}{2}}				 \Phi 	\|_{\F_f} \ .
%	\end{align}
%%	Additionally, if $\mathbb P$ is an orthogonal projection on $L^2(\R^d)$, for all $\Psi_1, \Psi_2 \in\F_f$		it holds true
%%	\begin{equation}
%%		| 	\< 	  \Psi_1 , \d \Gamma(	 \mathbb P \O\mathbb P 		)	 \Psi_2 		\>_{\F_f} |
%%		\leq
%%		\|	 \O	\|_{B(L^2)	}
%%		\| 		\d \Gamma (\mathbb  P )^{\frac{1}{2}}		 \Psi_1 	\|_{	\F_f	}
%%		\| 		\d \Gamma (\mathbb  P )^{\frac{1}{2}}		 \Psi_2	\|_{	\F_f	} \ .
%%	\end{equation}
%	\noindent (2) If $\O$ is Hilbert-Schmidt   with kernel $\O(x,y)$, then for all $\Psi \in \calD (\N_F^{\frac{1}{2}})$
%	\begin{align}
%		\|		 \d \Gamma (\O)	 \Psi 	\|_{\F_f}		& \leq 		\|	 \O	\|_{HS}		\|		 \N_F^{\frac{1}{2}}	\Psi 	\|_{\F_f}		\\
%		\|		\int_{	 \R^d \times \R^d 	}	 \O(x_1 , x_2 )  a_{x_1} a_{x_2} \d x_1 \d x_2  		\|_{\F_f}				 & \leq 		 \|	 \O	\|_{HS}		\|		 \N_F^{\frac{1}{2}}	\Psi 	\|_{\F_f}		\\
%		\|		\int_{	 \R^d \times \R^d 	}	 \O(x_1 , x_2 )  a^*_{x_1} a^*_{x_2} \d x_1 \d x_2  		\|_{\F_f}				 & \leq 		 \|	 \O	\|_{HS}		\|	( 	 \N_F+\1)  ^{\frac{1}{2}}	\Psi 	\|_{\F_f}	 \ .
%	\end{align}
%	\vspace{1mm}
%	\noindent (3) If $\O$ is trace-class with kernel $\O(x,y)$, then for all $\Psi \in \F_f$
%	\begin{align}
%		\|		 \d \Gamma (\O)	 \Psi 	\|_{\F_f}		& \leq 		\|	 \O	\|_{Tr}		\|		 	\Psi 	\|_{\F_f}			\\
%		\|		\int_{	 \R^d \times \R^d 	}	 \O(x_1 , x_2 )  a_{x_1} a_{x_2} \d x_1 \d x_2  		\|_{\F_f}			 & \leq 			\|	 \O	\|_{Tr}		\|		 	\Psi 	\|_{\F_f}			\\
%		\|		\int_{	 \R^d \times \R^d 	}	 \O(x_1 , x_2 )  a^*_{x_1} a^*_{x_2} \d x_1 \d x_2  		\|_{\F_f}			 & \leq 			\|	 \O	\|_{Tr}		\|		 	\Psi 	\|_{\F_f}		 \ .
%	\end{align}

%	All the estimates are proved in BPS2018, with the exception of the one containing $\mathbb P$.
%	The proof goes as follows.
%	For fixed $ n \in \mathbb N$, and  $ i \leq n  $, we	let $\O^{ (i)} \equiv \1^{i -1 } \otimes \O \otimes \1^{n-i}$
%	act on $\bigwedge_{i=1}^n L^2(\R^d)$.
%	We find
%	\begin{align}
%		\nonumber
%		| 	\< 	  \Psi_1 , \d \Gamma(	 \mathbb P \O\mathbb P 		)	 \Psi_2 		\>_{\F_f} |
%		& 	\leq
%		\sum_{n\geq1}
%		\sum_{i=1}^n
%		| 	 \< 	  \Psi_1,   \mathbb P ^{(i)}  \O^{(i)}   \mathbb P ^{(i)}  \Psi_2     	 \>_{\F_f}		|		\\
%		\nonumber
%		&  \leq
%		\|	 \O	\|_{B(L^2)}
%		\sum_{n\geq1}
%		\sum_{i=1}^n
%		\|	  \mathbb P ^{(i)} \Psi_1	\|_{\F_f}
%		\|	  \mathbb P ^{(i)} \Psi_2	\|_{\F_f}		\\
%		\nonumber
%		&	\leq
%		\|	 \O	\|_{B(L^2)}
%		\sum_{n\geq1}
%		\Big(
%		\sum_{i=1}^n
%		\|	  \mathbb P ^{(i)} \Psi_1	\|_{\F_f}^2
%		\Big)^{\frac{1}{2}}
%		\Big(
%		\sum_{i=1}^n
%		\|	  \mathbb P ^{(i)} \Psi_2	\|_{\F_f}^2
%		\Big)^{\frac{1}{2}}		\\
%		\nonumber
%		%&  =
%		%\|	 \O	\|_{B(L^2)}
%		% 	\sum_{n\geq1}
%		%\Big(
%		%\sum_{i=1}^n
%		%\|	  \mathbb P ^{(i)} \Psi_1	\|_{\F_f}^2
%		%\Big)^{\frac{1}{2}}
%		%\Big(
%		%\sum_{i=1}^n
%		%\|	  \mathbb P ^{(i)} \Psi_2	\|_{\F_f}^2
%		%\Big)^{\frac{1}{2}}		\\
%		\nonumber
%		&
%		\leq
%		\|	 \O	\|_{B(L^2)}
%		\Big(
%		\sum_{ n \geq 1 }
%		\sum_{i=1}^n
%		\|	  \mathbb P ^{(i)} \Psi_1	\|_{\F_f}^2
%		\Big)^{\frac{1}{2}}
%		\Big(
%		\sum_{ n \geq 1 }
%		\sum_{i=1}^n
%		\|	  \mathbb P ^{(i)} \Psi_2	\|_{\F_f}^2
%		\Big)^{\frac{1}{2}}		\\
%		&  =
%		\|	 \O	\|_{B(L^2)}
%		\<	   \Psi_1, \d \Gamma (\mathbb P)	\Psi_1				\>_{\F_f}^{		\frac{1}{2}	}
%		\<	   \Psi_2, \d \Gamma (\mathbb P)	\Psi_2				\>_{\F_f}^{		\frac{1}{2}	}\  .
%	\end{align}
%	In view of $\|	  \d \Gamma (\mathbb P )^{\frac{1}{2}} \Psi 	 \|_{\F}^2  =  \< \Psi, \d \Gamma (\mathbb P) \Psi 		 \>_\F$, the proof is then finished.

%\begin{lemma}[Boson estimates]
%	\label{lemma boson estimates}
%	Let $\O : L^2(\R^d)\rightarrow L^2(\R^d)$ be a bounded operator,
%	and let $\phi \in L^2(\R^d)$.
%	Then, the following holds true
%	\vspace{1mm}
%	\noindent (1) For all $\Phi \in \calD (\N_B^{\frac{1}{2}})$
%	\begin{align}
%		\|	  b(\phi)  \Phi		\|_{\F_b}			&  \leq \|	 \phi	\|_{L^2}		\|   \N_B^{\frac{1}{2}}  \Phi		\|_{	\F_b	}			\\
%		\|	  b^*(\phi)  \Phi		\|_{\F_b}			&  \leq \|	 \phi	\|_{L^2}		\|    (  \N_B+1)^{\frac{1}{2}}  \Phi		\|_{	\F_b	}		 \ .
%	\end{align}
%	\noindent (2) For all $\Phi \in \calD(\N_B )$
%	\begin{align}
%		\|	 \d \Gamma (\O )	 \Phi	\|_{\F_b} \leq 	  \| \O		\|_{B(L^2)}		\| \N_B	 \Phi		\|_{	 \F_b	} \ .
%	\end{align}

First, we turn to the proof of the above commutator estimates.
Subsequently, we show how they imply Proposition \ref{prop number estimates truncated}.

\begin{proof}[Proof of Lemma \ref{lemma commutator estimates}]
First, we prove the estimates for $\N_F$.
Secondly, we prove the estimates for $\N_B$.
Here, we only prove the case for $\theta = 0 $.
For general $\theta $ it suffices to insert an indentity
$\1 = (\N + 4  )^\theta (\N + 4 )^{1-\theta}$
and use the pull through formulas
$  \N a_x = a_x (\N + \1 )$,
and $\N b_y = b_y (\N+ \1 )$
on each term of $\widetilde \L (t). $

\vspace{1mm}

\noindent 	\textit{Proof of Eq. \eqref{prop estimate fermion}}.
Using the relations  $[\N_F, a^*(g)] = + a^*(g)$ and $[\N_F , a(g)] = - a(g)$
one is able to calculate the commutator
\begin{align}
\nonumber
[  \widetilde \L(t) ,\N_F]
&  +
\lambda \sqrt N
[  \L_{2,1}(t), \N_F]
\lambda
[\L_{2,2} (t), \N_F ]
\end{align}
\begin{align}
\nonumber
[  \L_{2,1}(t), \N_F]
&  +
\int_{	 \R^d \times \R^d 	} V(x-y )
a^*(u_x) a^*(	 \overline v _x) \otimes \big(	  \vp(y) b_y^* + \overline{\vp(y)}  b_y					\big)	\d x \d y 		\\
\label{prop eq 1}
&  -  	 2
\int_{	 \R^d \times \R^d 	} V(x-y )
a(	 \overline v _x)		a (u_x)  \otimes \big(	  \vp(y) b_y^* + \overline{\vp(y)}  b_y					\big)	\d x \d y 		\  , 	\\
\nonumber
[\L_{2,2 }(t) , \N_F ]
= &   +
\int_{	 \R^d \times \R^d 	} V(x-y )
a^*(u_x) a^*(	 \overline v _x) \otimes  b_y^* b_y 	\d x \d y 		\\
\label{prop eq 2}
&  -  	 2
\int_{	 \R^d \times \R^d 	} V(x-y )
a(	 \overline v _x)		a (u_x)  \otimes	 b_y^* b_y 	\d x \d y 	 \  .
\end{align}
Let us first estimate the terms in \eqref{prop eq 1}.
To this end, we use a Fourier decomposition
$V(x)  =   (2\pi)^{-d/2}	\int_{\R^d }	 \hat V (\xi)	 e_\xi(x) 	 \d \xi 	$,
where $e_{\xi} (x) \equiv \exp[ i x \cdot \xi ]$ ,
to find that
\begin{align}
[\L_{2,1}(t) & ,   \N_F] \\
\nonumber
& 	=
\frac{2}{(2\pi)^{\frac{d}{2}}}
\int_{	\R^d	}
\hat V ( 		\xi)
\Big(    \int_{	 \R^d \times \R^d 	}   [   u e_\xi v^*   ] (x_1 , x_2)	 a_{x_1}^* a_{x_2}^* \d x_1 \d x_2 	 \Big)
\otimes
\big(
b^*[ e_{-\xi} \vp 		]
b[ e_{ \xi} \vp 		]
\big) \d \xi 		\\
\nonumber
&  -
\frac{2}{(2\pi)^{\frac{d}{2}}}
\int_{	\R^d	}
\hat V ( 		\xi)
\Big(    \int_{	 \R^d \times \R^d 	} [   v e_\xi u   ] (x_1 , x_2)	 a_{x_1} a_{x_2} \d x_1 \d x_2 	 \Big)
\otimes
\big(
b^*[ e_{-\xi} \vp 		]
b[ e_{ \xi} \vp 		]
\big) \d \xi 	  \  .
\end{align}
Thus, we use the estimates contained in Lemma \ref{lemma fermion estimates} and \ref{lemma boson estimates}
and the Cauchy-Schwarz inequality for $\<  \cdot , \cdot \>_\F $
to find that
there exists a constant $C>0$
such that
for all $\Psi ,\Phi \in \F$
there holds
\begin{align}
\nonumber
|			 \langle	 \Psi , [  \L_{2,1} (t),  & \N_F]  \Phi  	\rangle_\F		|    \\
\nonumber
&  \leq
\|	 V	\|
\sup_{\xi \in\R^d}   \<\xi\>^{-1  }
\|	 u e_\xi v^*	\|_{HS}
\|   (\N_F+\1 )^{\frac{1}{2}} \otimes \1 	 \Psi 	\|_\F
\|	 \vp 	\|_{L^2} \| \1 \otimes (\N_B+\1)^{\frac{1}{2}}	 \Phi 	\|_\F  \  \\
& \leq
\|	 V	\|
\sqrt{M\hbar} \exp(Ct)
\|     (\N_F+\1 )^{\frac{1}{2}} \otimes \1 	 \Psi 	\|_\F   \| \1 \otimes   (\N_B+\1 )^{\frac{1}{2}}	 \Phi 	\|_\F  \  	 \ ,
\label{L 21 NF}
\end{align}
where in the last line we have used
$	 \|	 u e_\xi v^*	\|_{HS}
\leq    C \sqrt{M\hbar}   \< \xi \>  e^{Ct } $
and $\|	  \vp	\|_{L^2} =1 $.
Similarly, the second term in \eqref{prop eq 2}
may be expanded in its Fourier  coefficients--we find that
\begin{align}
[\L_{2,2}(t)  ,   \N_F]
\nonumber
& 	=
\frac{2}{(2\pi)^{\frac{d}{2}}}
\int_{	\R^d	}
\hat V ( 		\xi)
\Big(    \int_{	 \R^d \times \R^d 	}   [   u e_\xi v^*   ] (x_1 , x_2)	 a_{x_1}^* a_{x_2}^* \d x_1 \d x_2 	 \Big)
\otimes
\d\Gamma[e_{-	 \xi 	}]
\d \xi 		\\
&  -
\frac{2}{(2\pi)^{\frac{d}{2}}}
\int_{	\R^d	}
\hat V ( 		\xi)
\Big(    \int_{	 \R^d \times \R^d 	} [   v e_\xi u   ] (x_1 , x_2)	 a_{x_1} a_{x_2} \d x_1 \d x_2 	 \Big)
\otimes
\d\Gamma[e_{-	 \xi 	}]
\d \xi 	  \  .
\end{align}
Thus, we use the estimates contained in Lemma \ref{lemma fermion estimates} and \ref{lemma boson estimates}
and the Cauchy-Schwarz inequality for $\<  \cdot , \cdot \>_\F $
to find that
there exists a constant $C>0$
such that
for all $\Psi \in \calD (\N_F) \cap \calD(\N_B)$
there holds
\begin{align}
\nonumber
|			 \langle	 \Psi , [  \L_{2,2} (t),  \N_F]  \Phi  	\rangle_\F		|
& =
|			 \langle      (\N_B + \1 )^{\frac{1}{2}}  	  \Psi , [  \L_{2,2} (t),  \N_F]  	  (\N_B + \1 )^{ - \frac{1}{2}}  	\Phi	\rangle_\F		|     \\
\nonumber
&  \leq
\|	 V	\|
\sup_{\xi \in\R^d}   \<\xi\>^{-1  }
\|	 u e_\xi v^*	\|_{Tr}
\|   	  (\N_B + \1 )^{\frac{1}{2}}  \Psi 	\|_\F
\|	     \vp 	\|_{L^2}
\| \1 \otimes \N_B    (\N_B + \1 )^{- \frac{1}{2}}   \Phi 	\|_\F  \  \\
& \leq
\|	 V	\|
M  \hbar  \exp(Ct )
\|   	(\N_B + \1 )^{	\frac{1}{2}	} \Psi 	\|_\F
\|  \1 \otimes ( \N_B  +\1)^{	 \frac{1}{2}	} \Phi 	\|_\F   \  ,
\label{L 22 NF}
\end{align}
where in the last line we have used  $	 \|	 u e_\xi v^*	\|_{Tr}   \leq C M  \hbar \<\xi \>e^{Ct } $
and $\|	  \vp	\|_{L^2} =1 $.

\vspace{1mm}

We gather the estimates contained in Eqs. \eqref{L 21 NF} and \eqref{L 22 NF},
and use     the basic estimates $\N_F \otimes \1 \leq \N $,
$ \1 \otimes \N_B \leq \N $
together with Young's inequality $ab \leq a^2/2 + b^2/2$ for $a,b\geq0. $
This finishes the proof.

\vspace{1mm}

\noindent 	\textit{Proof of Eq. \eqref{prop estimate boson}}.
Using the relations  $[\N_B, b^*(g)] = + b^*(g)$ and $[\N_B , b(g)] = +  b(g)$
one is able to calculate the commutator
\begin{align}
\nonumber
[  \widetilde\L(t) ,\N_B]
&  +
\lambda \sqrt N
[\widetilde \L_{2,1}(t), \N_B]
\end{align}
\begin{align}
\nonumber
[ \widetilde \L_{2,1}(t), \N_B]
&	=
\chi(\N_B \leq \hbar M )
\int_{	 \R^d \times \R^d 	} V(x-y )
a^*(u_x)  a(u_x )\otimes \big(	  \vp(y) b_y^* + \overline{\vp(y)}  b_y					\big)	\d x \d y 		\\
\nonumber
&   +
\chi(\N_B \leq \hbar M )
\int_{	 \R^d \times \R^d 	} V(x-y )
a^*(	 \overline v _x)		a(	 \overline v _x)	 \otimes \big(	  \vp(y) b_y^* + \overline{\vp(y)}  b_y					\big)	\d x \d y 	 	\\
\nonumber
\int_{	 \R^d \times \R^d 	} V(x-y )
a^*(u_x) a^*(	 \overline v _x) \otimes \big(	  \vp(y) b_y^* + \overline{\vp(y)}  b_y					\big)	\d x \d y 		\\
&   +  	 2
\int_{	 \R^d \times \R^d 	} V(x-y )
a(	 \overline v _x)		a (u_x)  \otimes \big(	  \vp(y) b_y^* + \overline{\vp(y)}  b_y					\big)	\d x \d y 		\ .
\end{align}
Similarly as before, a Fourier decomposition
for the interaction potential yields
\begin{align}
[\widetilde \L_{2,1}(t), \N_B]
&	=
\label{prop eq 3}
\frac{2}{(2\pi)^{\frac{d}{2}}}
\int_{	\R^d	}
\hat V ( 		\xi)
\d \Gamma [u e_\xi u ]
\otimes		\chi(\N_B \leq \hbar M )
\big(
b^*[ e_{-\xi} \vp 		]
b[ e_{ \xi} \vp 		]
\big)
\d \xi 	\\
\label{prop eq 4}
&   +  	 	\frac{2}{(2\pi)^{\frac{d}{2}}}
\int_{	\R^d	}
\hat V ( 		\xi)
\d \Gamma [ v^* e_\xi v ]
\otimes
\chi(\N_B \leq \hbar M )	\big(
b^*[ e_{-\xi} \vp 		]
b[ e_{ \xi} \vp 		]
\big)
\d \xi 	\\
\label{prop eq 5}
\frac{2}{(2\pi)^{\frac{d}{2}}}
\int_{	\R^d	}
\hat V ( 		\xi)
\Big(    \int_{	 \R^d \times \R^d 	}
[   u e_\xi v^*   ] (x_1 , x_2)	 a_{x_1}^* a_{x_2}^* \d x_1 \d x_2 	 \Big)
\otimes
\big(
b^*[ e_{-\xi} \vp 		]
b[ e_{ \xi} \vp 		]
\big) 	\d \xi 	\\
\label{prop eq 6}
&   +  	 	\frac{2}{(2\pi)^{\frac{d}{2}}}
\int_{	\R^d	}
\hat V ( 		\xi)
\Big(    \int_{	 \R^d \times \R^d 	} [   v e_\xi u   ] (x_1 , x_2)	 a_{x_1} a_{x_2} \d x_1 \d x_2 	 \Big)     \otimes  	\big(
b^*[ e_{-\xi} \vp 		]
b[ e_{ \xi} \vp 		]     		 \big) 	\d \xi 	\ .
\end{align}
Up to an overall minus sign, the terms contained in
Eqs. \eqref{prop eq 5} and \eqref{prop eq 6}
have already been estimated above.
Thus, it suffices to estimate the two terms contained in
Eqs. \eqref{prop eq 3} and \eqref{prop eq 4}.
Let us start from the one in Eq. \eqref{prop eq 3} and, for simplicity, let us only show how to bound the
contribution arising from $b[e_\xi \vp]$--the other one is analogous.
Let us fix $\xi \in \R^d$, and let  $\Psi , \Phi \in \F $.
Denoting $\chi_M \equiv 		\chi(\N_B \leq \hbar M )$,
we find using
and   \textit{pull-through formula}
$\N_B b_x = b_x (\N_B + \1)$
\begin{align}
|	 \<   \Psi  ,
\d \Gamma[   u e_\xi u  ]  \otimes \chi_M
b[e_\xi \vp ]
\Phi  	 \>_\F 		|
& 		=
\nonumber
|	 \langle  \chi_M   \N_B^{\frac{1}{4}	}   \Psi  ,   \d \Gamma[   u e_\xi u  ]  \otimes \chi_M
b[e_\xi \vp ]   (\N_B+\1)^{-\frac{1}{4}	} 	 \Phi  	 \rangle_\F 		| \\
\nonumber
\leq
\|		\N_F^{\frac{1}{2}} 		 \chi_M    \N_B^{\frac{1}{4}	}   \Psi  	 	\|_\F
\|			\N_F^{\frac{1}{2}} 		  \chi_M   b[e_\xi \vp ]   (\N_B+\1)^{-\frac{1}{4}	} 	 \Phi  	\|_\F    \\
& \leq
\sqrt{ M \hbar }
\|	 \N_F^{\frac{1}{2}}	 \otimes \1 		\Psi	\|_{\F}
\|	 \N_F^{\frac{1}{2}}		 \otimes \1 	\Phi	\|_{\F}
%		& = \nonumber
%				\|		 \chi_M    \N_B^{\frac{1}{4}	}      	\d \Gamma[u]^{\frac{1}{2}}     \Psi  	 	\|_\F
%		\|		 \chi_M   b[e_\xi \vp ]   (\N_B+\1)^{-\frac{1}{4}	}  	\d \Gamma[u]^{\frac{1}{2}}   	 \Phi  	\|_\F    \\
%		& \leq
%	C
%	\sqrt{M \hbar }
%	 \< 	 \Psi , \d \Gamma[u]\Psi	\>_\F \  \\
%	 & \leq
%	 	C
%	 \sqrt{M \hbar }
%\|  \N_F^{\frac{1}{2}}\otimes\1 	 \Psi	\|_\F^2 \ ,
\end{align}
where we used Lemma \ref{lemma fermion estimates} and \ref{lemma boson estimates}
together with
$		\|  u e_\xi u 	 \|_{B(L^2)} \|  \leq 1 $.
Similarly, we find
\begin{equation}
|	 \<   \Psi  ,   \d \Gamma[   v^* e_\xi v	  ]  \otimes \chi_M  b[e_\xi \vp ]	 \Psi 	 \>_\F 		|
\leq
\sqrt{ M \hbar }
\|	 \N_F^{\frac{1}{2}}	 \otimes \1 		\Psi	\|_{\F}
\|	 \N_F^{\frac{1}{2}}		 \otimes \1 	\Phi	\|_{\F}
\end{equation}
The proof if finished, once we gather our estimates and use
the elementary  inequalities
$\N_F \otimes \1 \leq \N $,
$ \1 \otimes \N_B \leq \N $
together with Young's inequality $ab \leq a^2/2 + b^2/2$ for $a,b\geq0. $
\end{proof}

%Before we turn to the proof of Proposition \ref{prop number estimates}, we need to establish
%the following result that shows that the charge operator $\Q(t)$
%is a conserved quantity.
%\begin{lemma}[Conservation of charge]
%	\label{lemma charge}
%	Let $\U(t,s)$ be the fluctuation dynamics, and let $\calQ (t)$
%	be the charge operator defined in Eq. \eqref{charge}.
%	Then, it holds that for all $t,s\in\R$ that
%	\begin{equation}
%		\Q(t)   \U(t,s) = \U(t,s) \Q (s) \ .
%	\end{equation}
%	Let us fix $t,s\in\R$. First, a straightforward computation shows that
%	\begin{equation}
%		\label{conjugation of N}
%		\calR_t^*  \N_F  \calR_t = \Q(t) +  M \ .
%	\end{equation}
%	Next, we recall that the original many-body dynamics
%	preserves the number of fermions in the system.
%	In other words, there holds
%	\begin{equation}
%		e^{	   \frac{i (t-s)}{\hbar }	\H	}
%		( \N_F \otimes \1)
%		e^{	   \frac{ - i (t-s)}{\hbar }	\H	}	 =
%		( \N_F \otimes \1)  \ .
%	\end{equation}
%	We apply $\calR_s^* \W_s^*$ to the left and $\calR_s \W_s $ to the right to find that
%	\begin{eqnarray}
%		\U^*(t,s)
%		\calR_t^*	( \N_F \otimes \1) \calR_t
%		\U(t,s )
%		=
%		\calR_s^*	( \N_F \otimes \1) \calR_s \ .
%	\end{eqnarray}
%	It suffices now to use the relation \eqref{conjugation of N}, cancel the constant terms, and apply $\U(t,s)$ on both sides of the resulting equation.

\begin{proof}[Proof of Proposition \ref{prop number estimates truncated}]
Let us fix $ k \in \mathbb N  $ and   recall that
$\i \hbar \partial_t \widetilde \U(t,s ) = \widetilde\L(t)\widetilde \U(t,s)$.
Hence, we may compute the time derivative for $t,s \in\R $
\begin{align}
\nonumber
\frac{d}{dt}
\widetilde \U^* (t,s) \N^k 	\U(t,s)
&  =
\frac{1}{ \i \hbar }
\widetilde \U^* (t,s) [  \widetilde \L (t ) , \N^k   ] 		\widetilde \U (t,s) 	 \\
&   =
\frac{1}{ \i \hbar }
\sum_{i =1 }^k
\widetilde \U^* (t,s)  \N^{ i -1 }  [  \widetilde \L (t ) , \N   ]  \N^{ k - i }		\widetilde \U (t,s)  \ .
\end{align}
Taking $\Psi \in \F$, we can then estimate that
\begin{align}
\frac{d}{dt}	 \<	  \Psi , 	\widetilde \U^* (t,s)  \N^k \widetilde \U (t,s) \Psi 	 \>_\F
\leq
\frac{1}{\hbar }
\sum_{i =1 }^k
| \langle
\N^{ i -1 }   \widetilde \U (t,s) \Psi,     [  \widetilde \L (t ) , \N   ]  \N^{ k - i }		\widetilde \U (t,s)  \Psi \rangle_\F  | \ .
\end{align}
Fix $ i  =1, \ldots , k $ and let $\theta = 1 + k -2 i \in \R $.
In view of Lemma
\ref{lemma commutator estimates} there exists $C > 0 $ such that
\begin{align}
\nonumber
| \langle
& 	\N^{ i -1 }     \widetilde \U (t,s) \Psi	 		,     [  \widetilde \L (t ) , \N   ]
\N^{ k - i }		\widetilde \U (t,s)  \Psi
\rangle_\F  |  	\\
\nonumber
& 	\leq
C e^{Ct }
\lambda
(  \sqrt{NM \hbar }   + M \hbar )
\Big(
\|	 (\N + \1 )^{	 \frac{1 + \theta}{2}	} 	\N^{ i -1 }   \widetilde \U (t,s) \Psi 	\|_{\F}^2
%\< \Psi ,  (\N  +1  )^{1+\theta } \Psi \>_\F
\|	 (\N + \1 )^{	 \frac{1 - \theta}{2}	}		 \N^{ k - i }		\widetilde \U (t,s)     \Psi	\|_{\F}^2
\Big)   , \\
\nonumber
& \leq
C e^{Ct }
\lambda
(  \sqrt{NM \hbar }   + M \hbar )
\|	 (\N + \1 )^{\frac{k}{2}	} 	  \widetilde \U (t,s) \Psi 	\|_{\F}^2 	  \\
&   =
C e^{Ct }
\lambda
(  \sqrt{NM \hbar }   + M \hbar )
\langle 		\Psi , \widetilde \U^*(t,s)  ( \N^k+ \1)  \widetilde \U(t,s) \Psi 	 \rangle_\F 	  \ .
\end{align}
The proof is then finished after we apply the Gr\"onwall inequality and
use the initial condition $\widetilde\U(s,s)=\1 . $
\end{proof}

\subsection{Proof of Theorem \ref{prop number estimates}}
Let us now give a proof of Theorem \ref{prop number estimates}.
Essentially, we shall compare the number estimates generated by   $ (t,s)$
and $(t,s)$.
The latter have already been established in Proposition \ref{prop number estimates}.
Additionally, we shall need the following \textit{weak bounds}
on the growth of particle number with respect to the original fluctuation dynamics.

\begin{lemma}[Weak number estimates]
\label{lemma weak estimates}
Let $\U(t,s)$ be the fluctuation dynamics, defined in  \eqref{unitary U}.
Then, the following statements hold true.
\begin{enumerate}[leftmargin=*]
\item
For all $\ell \in \mathbb N$ there is a constant $C>0$ such that for all $t,s\in \R$ and $\Psi \in \F$
there holds
\begin{equation}
\|    ( \N_F+M)^\ell \U(t,s) \Psi \|_\F
\leq
\|    ( \N_F+M)^\ell   \Psi \|_\F
\end{equation}

\item
For all $\ell \in \mathbb N$ there is a constant $C>0$ such that for all $t,s\in \R$ and $\Psi \in \F$
there holds
\begin{equation}
\|    ( \N_B+N )^\ell \U(t,s) \Psi \|_\F
\leq
\|    ( \N_B+N)^\ell   \Psi \|_\F
\end{equation}
\end{enumerate}
\end{lemma}

\begin{proof}
We only give a proof for the fermion number operator since the proof for bosons is similar;
we refer the reader to  \cite[Lemma 3.6]{RodnianskiSchlein2009}
for the situation in which only bosons are present (interactions do not change the proof).

\vspace{1mm}
Let us then  consider the particle-hole transformation $\calR_t$ in the definition of $\U(t,s)$ and notice that
\begin{equation}
\calR_t (\N_F+M)^\ell \calR_t^* = [ \calR_t (\N_F+M) \calR_t^* ]^\ell  =
\d \Gamma[v(t)]
\d\Gamma[u (t)]
+ 2 M
)^\ell  \ .
\end{equation}
The terms
$\d\Gamma[u (t)] $ and $ \d \Gamma[v(t)] $
can be estimated using Lemma \ref{lemma fermion estimates}.
Namely,  it follows that there exists $c_0$ such that
for all $\Phi\in \F$
\begin{equation}
\|
\d \Gamma[v(t)]
\d\Gamma[u (t)]
+ 2 M
) \Phi
\|_\F
\leq
c_0 \|
(\N_F + M )
\Phi\|_\F \ .
\end{equation}
Consequently, since $ [ \d\Gamma_F[J], \N_F]=0$
for $J=u(t)$ and $J=v(t)$,
an $\ell$-fold application of the previous estimate yields
\begin{align}
\label{l fold application}
\|
( \N_F     +M)^\ell
\calR_t^*
\Phi
\|_\F
& \leq
\|       (\N_F + M )^\ell \Phi     \|_\F \  ,
\qquad \forall t \in \R\  .
\end{align}
Here  we have used the fact that $\calR_t$ is a unitary transformation on Fock space.
We conclude using the fact that $(\N_F + M)$
commutes with the Schr\"odinger dynamics $\exp(- i(t-s) \H / \hbar)$,
and the bosonic  Weyl operator $\W_t$.
That is
\begin{align}
\nonumber
\|       (\N_F + M )^\ell   \U (t,s) \Psi        \|_\F
&  =
\|             (\N_F + M )^\ell    \calR_t^*        \W_t^*e^{-i(t-s)\H /\hbar} \calR_s \W_s \Psi            \|_\F        \\
\nonumber
& \leq
\|        (\N_F + M)^\ell  \calR_s \Psi    \|_\F \\
& \leq
c_0^{2\ell }
\|        (\N_F + M)^\ell   \Psi    \|_\F \  .
\end{align}
where in the last line we have used that $\calR_s = \calR_{-s}^*$
together with the    estimate \eqref{l fold application}.
This finishes the proof.
\end{proof}

\begin{proof}[Proof of Theorem \ref{prop number estimates}]
Since the proof is essentially the  same one for fermions and bosons,
we   only present the proof of the former in full details,  and point out the differences with respect to the latter.
To this end,  	let us now fix  $\Psi \in \F$   together with $t,s \in \R $
and $k,\ell \in \mathbb N. $
We start by computing that
\begin{align}
\nonumber
\| \N_F^{\frac{\ell}{2}	} \U(t,s) \Psi		\|_\F^2
& 	 =
\langle
\Psi , \U (t,s )
\N_F^\ell
\widetilde \U(t,s)
\Psi 		\rangle_\F
\langle 	 \Psi , \U (t,s )
\N_F^\ell
\Big( 				\U(t,s) 	- 			\widetilde \U(t,s) 				\Big)
\Psi 		\rangle_\F 	\\
& \leq
\|	 \Psi 	\|_\F
\|	 	 \N_F^\ell   \widetilde \U(t,s) \Psi 		\|_\F
\|	 \N_F^\ell 		\U (t,s) 		\Psi 	\|_\F
\|	 	 (	  \widetilde \U(t,s)			-  \U(t,s))		 \Psi 		\|_\F
\label{eq 1} \ .
\end{align}
We now estimate  the two terms in \eqref{eq 1}.

\vspace{1mm}