For the first one, we use the fact that the truncated dynamics satisfies  the estimate
contained in Proposition \ref{prop number estimates truncated}.
Namely,  for a constant $C>0 $
\begin{eqnarray}
\|	 	 \N_F^\ell   \widetilde \U(t,s) \Psi 		\|_\F
\leq
%			{  \exp     \bigg[  C   \lambda  \sqrt \frac{NM}{\hbar }
%				\bigg( 1 + \sqrt{\frac{\hbar M}{  N}}    \bigg)	\exp |t-s|      \bigg]		}
\|	 	  ( \N+\1)^\ell     \Psi 		\|_\F  \ .
\end{eqnarray}

\vspace{1mm}

For the second term  in \eqref{eq 1}
we use the weak number estimates from Lemma \ref{lemma weak estimates}.
Namely, for a constant $C = C(\ell)$
\begin{equation}
\|	\N_F^\ell \U(t,s)	 \Psi \|_\F
\leq
\|	    (\1 + \N_F / M)^\ell \Psi	\|  \  , \qquad \forall t, s \in \R \   .
\end{equation}
Next,  we study the difference between
the original and the truncated dynamics.
That is,
we use Duhamel's formula
to find that
\begin{align}
\widetilde \U(t,s)  -  	\U(t,s )
-i  \delta
\int_s^t
\U(t, r  )
\Big(
\widetilde \L_{2,1} (r )  - \L_{2,1}(r )
\Big)
\widetilde \U(r ,  s ) \d r   ;
\end{align}
where we have used the fact that
$ \widetilde \L (t)   - \L(t) =  \lambda \sqrt N  ( \widetilde \L_{2,1} (t)  -  \L_{2,1} (t))$,
and we have collected $\delta = \lambda \sqrt N  / \hbar $ as a pre factor.
Next, for $\Phi = \U(t,s)\Psi  \in \F$ one may estimate using Lemma  \ref{lemma fermion estimates}
and \ref{lemma boson estimates} that
for all  $ k \geq 0 $
(here, we omit time labels for notational convenience)
\begin{align}
\nonumber
\|	  (	 \widetilde \L_{2,1}	  	-   \L_{2,1}  	)	 \Phi	\|_{\F}
& \leq
\int_{\R^d}
|	 \hat V(\xi)	|
\Big\|
\d \Gamma [u e_\xi u ]	 \otimes  \chi( \N_B \geq \hbar M 	)\Big(	 b[e_{\xi} \vp]  + b^*[ e_{-\xi} \vp ]		\Big)	  \Phi
\Big\|_\F
\d \xi
\\
\nonumber
& +
\int_{\R^d}
|	 \hat V(\xi)	|
\Big\|
\d \Gamma [ v^* e_\xi v ]	 \otimes \chi( \N_B \geq \hbar M 	)\Big(	 b[e_{\xi} \vp]  + b^*[ e_{-\xi} \vp ]		\Big) \Phi 		\Big\|_\F  \d \xi			\\
\nonumber
& \leq 2
\|	 \hat V	\|_{L^1 }
\|	   \N_F \otimes 	 \chi( \N_B \geq \hbar M 	)(\N_B +\1)^{\frac{1}{2}} \Phi	\|_\F 		\\
\nonumber
& \leq
\frac { 2
\|	 \hat V	\|_{L^1 }}{		(\hbar M)^k 		}
\|	   \N_F \otimes 	 \chi( \N_B \geq \hbar M 	)(\N_B +\1)^{k+ \frac{1}{2}} \Phi	\|_\F 	  \\
& \leq
\frac { C
\|	 \hat V	\|_{L^1 }}{		(\hbar M)^k 		}
\|	  (\N + \1 )^{k + \frac{3}{2}} \Phi	\|_\F 	 \   .
\end{align}
We put our last three estimates together to find that
thanks to Proposition \ref{prop number estimates truncated}
\begin{align}
\|	 	\Big(  \widetilde \U (t,s)      - \U(t,s) \Big)  		\Psi 		\|_\F
\leq
\frac{\lambda \sqrt N }{\hbar}
\frac{1}{( \hbar M)^k}
\|	 (1 + \N)^{k+3/2} \Psi	\|_{\F}  \ .
\end{align}
Putting our estimates together, we find that
there exists $C>0 $ such that
\begin{align}
\nonumber
\| \N_F^{\frac{\ell}{2}	} \U(t,s) \Psi		\|_\F^2
& 	  \leq
K(t-s )
\Big[
\|	 \Psi	 \|_\F  \| (\N + \1 )	 \Psi	\|_\F
\Theta_{k,\ell}
\|  \big(\1 + \frac{\N_F}{ M } \big)^\ell   \Psi		\|_\F
\|	 	  ( \N+\1)^{k + \frac{3}{2}	}    \Psi 		\|_\F
\Big]
\end{align}
where $\Theta_{k,\ell}$ is as in the statement of the theorem.

\vspace{1mm}

As for the bosons, the only modification comes from the weak number estimates
obtained from Lemma \ref{lemma weak estimates},
In the bosonic case,
one finds that for a constant $C = C(\ell)$
\begin{equation}
\|	\N_B^\ell \U(t,s)	 \Psi \|_\F
\leq
\|	    (\1 + \N_B /  N )^\ell \Psi	\|  \  , \qquad \forall t, s \in \R \  ,
\end{equation}
This finishes the proof. 	\end{proof}

\section{Second Quantization III: Proof of Theorem \ref{thm1}}
\label{section proof thm 1}

Let us now  turn to the proof our first main result.
We shall reduce the problem of  proving the estimates contained in Theorem \ref{thm1}, to that of proving number estimates for the fluctuation dynamics, as defined in Section \ref{section fluctuation}.
Then, we apply Theorem \ref{prop number estimates}.
The proofs we present are heavily inspired by the works
[55, 9]
and also
[15, 16],
and are adapted to the case at hand.
Let us remark here
that our proofs are shorter, because we do not aim at obtaining the optimal convergence rates
$1 / M $ for fermions, and $1/N$ for bosons.
Here instead, we content ourselves with $1/√(M)$ and $1/√(N)$, respectively.

\vspace{2mm}

Let us first introduce some notation.
Letting $(t,s)$ be the fluctuation dynamics,
we consider  the following
\textit{fluctuation vectors}
in Fock space $$
\begin{align}
\Omega_1 (t)
\equiv
\U(t,0) \Omega
\quad
\t{and}
\quad
\Omega_2(t)
\equiv
\U(t,0)
\Omega_F \otimes \W^*[\sqrt N \vp_0] \frac{b^*(\vp_0)}{\sqrt  N } \Omega_B \ .
\end{align}
They will be extremely useful in the proof of Theorem \ref{thm1}.
determines the fluctuations
of the system when the initial data
has a bosonic  coherent state.
On the other hand, $Ω_2(t)$
describes the fluctuations
when the initial data has a bosonic factorized state.
Let us collect  in the following lemma
some results concerning these vectors;
for reference, see [16]

\begin{lemma}[Properties of $\Omega_1$ and $\Omega_2$]
Let us denote from now on
\begin{equation}
d_N \equiv \frac{\sqrt{N!}}{N^{N/2} e^{-N/2}}
\end{equation}
Then, the following statementes hold true
\begin{enumerate}
\item  $P_N \W[\sqrt N \vp_0]\Omega_B
\frac{1}{d_N}
\frac{b^*(\vp_0)}{\sqrt{N!}}\Omega_B$
$P_N \equiv \1 (\N_B =N)$.

\item
\<   \Omega_2(t) , \Omega_1(t) \>_\F
\<   \Omega_2(0 ) , \Omega_1(0) \>_\F
\frac{1}{d_N} \ .

\item
There exists a constant $C>0$ such that
\|	 (\N_F+\1)^{-1/2}	 \Omega_2(0)	\|_\F
\leq
\frac{C}{d_N}  \ . $

\end{enumerate}
\end{lemma}

\begin{proof}[Proof of Theorem \ref{thm1}]
We recall  recall here
that one-particle reduced densities
have been re-cast in creation-
and annihilation operators in
\eqref{definition kernels}.

\vspace{1mm}
\noindent
\textit{Proof of \eqref{thm1 fermion estimate}}
First, we establish the result for fermions.
If we denote the time evolution of operators by
$a_x (t) = e^{it/\hbar \H} a_x e^{-it/\hbar \H}$
then we obtain
\begin{align}
\nonumber
\gamma_F(t;x_1 ,x_2)
&  =
\< \Psi(t) , a_{x_2}^*a_{x_1}\Psi (t) \>_\F 		\\
\nonumber
& =
\<
\calR_0 \otimes \frac{b^*(\vp_0)^N}{\sqrt{N!}} \Omega
a_{x_2}^* (t)  a_{x_1}  (t)
\calR_0   \otimes \frac{b^*(\vp_0)^N}{\sqrt{N!}} \Omega
\>_\F  \\
\nonumber
& =
\<
\calR_0   \otimes \frac{b^*(\vp_0)^N}{\sqrt{N!}} \Omega
a_{x_2}^* (t)  a_{x_1}  (t)
\calR_0   \otimes
\W[\sqrt N \vp_0 ]
\Omega
\>_\F  \\
& =
\<  \Omega_2(t)  ,
\calR_t^*
a_{x_2}^*  a_{x_1}
\calR_t
\Omega_1 (t)
\>_\F
\label{thm 1 eq 1 fermion}
\end{align}
%	As a reminder, we will use the conjugation relations
%	\begin{align}
%\label{conjugation fermions}
%		R_t^* a^*_x R_t
%		 =
%		a^* ( u_{t,x})
%		+
%		a( \overline{  v_{t,x} } )
%		\quad \t{and}\quad
%		R_t^* a_x R_t
%		 =
%		a ( u_{t,x})
%		+
%		a^*( \overline{  v_{t,x} } )
%	\end{align}
%for fermions, and for bosons.
%\label{conjugation bosons}
%	\W_t^* b_x  \W_t
%	=
%	b_x + \sqrt N \vp(t,x)
%	\qquad
%	\t{and}
%	\qquad
%	\W_t^* b^*_x  \W_t
%	=
%	b^*_x +  \overline{ \vp(t,x) } \ .
Next, we look at the conjugation relations
\begin{align}
\label{conjugation fermions}
\calR_t^*
\calR_t
a^* ( u_{t,x})
a( \overline{  v_{t,x} } )
\quad \t{and}\quad
\calR_t^* a_x \calR_t
a ( u_{t,x})
a^*( \overline{  v_{t,x} } )
\end{align}
and obtain, using $\<\bar v_x , \bar v_y \>= \omega(x,y)$
$\<\Omega_1 , \Omega_2\> = d_N$
\begin{align}
\nonumber
\gamma_F( t; x_1 , x_2)
-  \omega ( t;  x_1, x_2)
& 			=
d_N  \<   \Omega_2 (t ) ,  a^*(u_{x_2,t})  a(u_{x_1,t})  \Omega_1 (t)			\>_\F  \\
\nonumber
& 	-
\<   \Omega_2 (t ) ,  a^*(\overline{v}_{x_1,t})  a(\overline{v}_{x_2,t})  \Omega_1 (t)			\>_\F 	\\
\nonumber
& 	 			 +  d_N
\<   \Omega_2 (t ) ,  a^*(u_{x_2,t})  a^*(v_{x_1,t})  \Omega_1 (t)			\>_\F  \\
&  +
\<   \Omega_2 (t ) ,  a(\overline{v}_{x_2,t})a( u_{x_1,t}) \Omega_1 (t)\>_\F    \    .
\end{align}
Next, we let $\O $ be a compact operator in $L^2(\R^d)$
with kernel $\O(x_1 ,x_2) $.  The above identity now implies that
\begin{align}
\nonumber
\mathrm{Tr}  \O \big( \gamma_F (t) - \omega (t)	\big)
&   =
\<	 \Omega_2 (t) ,
\,  \d \Gamma_F [u_t \O u_t ]
\Omega_1 (t) 	 \>_\F  \\
\nonumber
& 			-  d_N
\<	 \Omega_2 (t) ,
\,  \d \Gamma_F [v_t^* \overline{\O^*} v_t  ]
\Omega_1 (t) 	 \>_\F 	\\
\nonumber
\<	 \Omega_2 (t) ,
\int_{ \R^d \times \R^d }
[v_t \O   u_t] (x_1, x_2)
a_{x_1}    a_{x_2}   \,
\Omega_1 (t)	\>_\F  \\
& +
\<	 \Omega_2  (t) ,
\int_{ \R^d \times \R^d }
[v_t \O   u_t]^* (x_1, x_2)
a_{x_1}^*    a_{x_2}^*
\,
\Omega_1 (t)	\>_\F
\end{align}
Next, we make use of Remark \ref{remark boundedness of N}.
In particular,
we consider $K_{1/2} \geq 1 $
to be the integer that controls of $\N^{1/2}$.
The Cauchy-Schwarz inequality
Lemma \ref{lemma fermion estimates}
then gives
\begin{align}
\nonumber
\big|
\mathrm{Tr}  \O \big( \gamma_F (t) - \omega (t)	\big)
\big|
& 	\leq
\|		 u \O  u 	\|_{B(\F)}
\| (\N + \1 )^{-K_{1/2}}	\Omega_2(t)	\|_{\F}
\|	  (\N +\1 )^{K_{1/2} + 1  }	  \Omega_1 (t)	\|_{\F } \\
\nonumber
& 	+
\|		 v \O  v 	\|_{B(\F)}
\| (\N + \1 )^{-K_{1/2}}	\Omega_2(t)	\|_{\F}
\|	  (\N +\1 )^{K_{1/2} + 1  }	  \Omega_1 (t)	\|_{\F }	\\
\nonumber
& 	+
\|		 v \O  u 	\|_{HS}
\| (\N + \1 )^{-K_{1/2}}	\Omega_2(t)	\|_{\F}
\|	  (\N +\1 )^{K_{1/2} + 1/2  }	  \Omega_1 (t)	\|_{\F }	\\
& 	+
\|		(vOu)^*	\|_{HS}
\| (\N + \1 )^{-K_{1/2}}	\Omega_2(t)	\|_{\F}
\|	  (\N +\1 )^{K_{1/2} + 1/2  }	  \Omega_1 (t)	\|_{\F } \ .
\end{align}
First, we control $\Omega_2(t)$ and then, we control $\Omega_1(t)$.
Indeed, Remark \ref{remark boundedness of N} implies that
\begin{align}
\nonumber
\| (\N + \1 )^{-K_{1/2}}	\Omega_2(t)	\|_{\F}
& 		  \leq
\| (\N + \1 )^{-K_{1/2}}   \U(t,0) (\N + \1 )^{1/2}	\|_{B(\F)}
\| (\N+\1)^{-1/2} 	\Omega_2 (0)	\|_{\F} \\
& 	  	 \leq
K(t) / d_N
\label{thm1 eq 2}
\end{align}
where we have used
\|	 (\N_F+1)^{-1/2}	 \Omega_2(0)	\|_\F
\leq
c/d_N   \ . $
On the other hand, Remark \ref{remark control omega} immediately
implies that
\begin{equation}
\label{thm1 eq 3}
\|	  (\N +\1 )^{K_{1/2} + 1  }	  \Omega_1 (t)	\|_{\F }
\leq K (t)
\end{equation}
provided one re-updates $K(t)$.
Finally, we use
$\|	 u \O u 	\|_{B(\F)} \leq \|	 \O	\|_{B(\F)}$,
$\|	 v \O v 	\|_{B(\F)} \leq \|	 \O	\|_{B(\F)}$,
$\|	 v \O u 	\|_{HS)} \leq    \sqrt M   \|	 \O	\|_{B(\F)}$,
$\|	   ( v \O u)^*	\|_{HS} \leq    \sqrt M   \|	 \O	\|_{B(\F)}$;
and the fact that the space of compact operators is dual to the space of trace-class operators to conclcude that
\begin{equation}
\|	 \gamma_F (t)  - \omega(t)		\|_{\tr}
\leq
K(t) \sqrt M \ .
\end{equation}
This finishes the proof of the estimate for fermions.

\vspace{1mm}
\noindent
\textit{Proof of \eqref{thm 1 boson estimate}.}
As for the bosons,  we start by looking at the identity
that is the analogous of
\eqref{thm 1 eq 1 fermion}.
Namely,  a similar argument shows that
\begin{align}
\gamma_B(t; y_1, y_2)
\< \Omega_1(t) ,   \W_t^* b_{y_2}^* b_{y_1}^* \W_t  \Omega_2(t)\>_\F \  .
\end{align}
Next, the conjugation relations
\begin{align}
\label{conjugation bosons}
\W_t^* b_x  \W_t
b_x + \sqrt N \vp_t(x)
\qquad
\t{and}
\qquad
\W_t^* b^*_x  \W_t
b^*_x +  \overline{ \vp_t(x) }
\end{align}
combined with the identity
$\<\Omega_2, \Omega_1\>= 1/ d_N $ now imply that
\begin{align}
\nonumber
\gamma_B(t; y_1, y_2)
\vp_t(y_1 )
\overline{\vp_t(y_2)}
&  =
\< 	\Omega_2(t) ,				b_{y_2}^*  b_{y_1}					\Omega_1(t) 		\> \\
\nonumber
&  +
\sqrt N
\< 	\Omega_2(t) ,			b^*_{y_2}			\vp_t(y_1) 			\Omega_1(t) 		\> \\
&  +
\sqrt N
\< 	\Omega_2(t) ,			b_{y_1}  \overline{ \vp_t(y_2)	} 					\Omega_1(t) 		\>  \ .
\end{align}
Next, we consider a trace-class operator
$\O$ with kernel $O(x,y)$
and integrate over $y_1, \ y_2$ to find that
\begin{align}
\nonumber
\Big|
\tr \Big(
\O (	\gamma_B(t) - N \ket{\vp_t} \bra{\vp_t}	)
\Big)
\Big|
& 	\leq
\<	 \Omega_2(t) , \d \Gamma_B [\O ] \Omega_1 (t)		 \>  \\
& + d_N
\sqrt N
\<	 \Omega_2(t) ,
\big(    b^*[\O \vp_t] + b[\O \vp_t ]    \big)
\Omega_1 (t)		 \>   \ .
\end{align}
Similarly as we did for fermions,
we consider $K_{1/2}$
to be the integer from Remark \ref{remark control omega}.
Thus,  the Cauchy-Schwarz inequality,
the pull-through formula
$\N b_y = b_y ( \N +1)$
Lemma \ref{lemma boson estimates} now imply
\begin{align}
\nonumber
\Big|
\tr \Big(
\O (	\gamma_B(t) - N  & \ket{\vp_t} \bra{\vp_t}	)
\Big)
\Big| 			\\
\nonumber
& 	\leq
\|	(\N+\1)^{-K_{\frac{1}{2}}}  \Omega_1(t)	\|_\F
\|	(\N+\1)^{K_{	\frac{1}{2}		}}    \d \Gamma_B[\O]      \Omega_1(t)	\|_\F \\
\nonumber
& +
d_N \sqrt N
\|	(\N+\1)^{-K_{\frac{1}{2}}}  \Omega_1(t)	\|_\F
\|	(\N+\1)^{K_{	\frac{1}{2}		}}      b^*[\O \vp_t]       \Omega_1(t)	\|_\F \\
\nonumber
& +
d_N \sqrt N
\|	(\N+\1)^{-K_{\frac{1}{2}}}  \Omega_1(t)	\|_\F
\|	(\N+\1)^{K_{	 \frac{1}{2}		}}    b[\O \vp_t ]       \Omega_1(t)	\|_\F  \\
\nonumber
& \leq
\|	\O 	\|_{B(L^2)}
\|	(\N+\1)^{-K_{\frac{1}{2}}}  \Omega_1(t)	\|_\F
\|	(\N+\1)^{K_{	\frac{1}{2}	} +1 }         \Omega_1(t)	\|_\F \\
& +
2 d_N
\sqrt N \|	 \O \vp_t	\|_{L^2}
\|	(\N+\1)^{-K_{\frac{1}{2}}}  \Omega_1(t)	\|_\F
\|	(\N+\1)^{K_{	\frac{1}{2}	}+ \frac{1}{2}   }         \Omega_1(t)	\|_\F  \ .
\end{align}
In order to conclude,
we control the right hand side of the above equation
with the estimates \eqref{thm1 eq 2} and \eqref{thm1 eq 3}.
Using that $\|	 \O \vp	\|_{L^2} \leq \|	 \O	 \|_{B(L^2)}	$
recalling the space of compact operators is dual to the space
of trace-class operators, we find
\begin{equation}
\|    \gamma_B(t)   -  N    \ket{\vp_t}  \bra{\vp_t}        		\|_{\tr}
\leq
K(t) \sqrt N \ .
\end{equation}
This finishes the proof of the Theorem.
\end{proof}

\section{Quantum Optimal Transportation: Proof of Theorem \ref{thm2}}
\label{section QOT}

In this section, we address the semi-classical limit of the coupled
Hartree-Hartree system first introduced in \eqref{hartree hartree}.
Here, we adopt the scaling regime presented in \eqref{scaling regime}.
Thus, we study the
of the equation
\begin{align}
\label{hartree hartree 2}
\begin{cases}
& i\hbar \partial_t \omega^\hbar
= \Big[ -  \frac{1}{2}\hbar^2 \Delta + V*\rho_B , \omega^\hbar    \Big]  \\
& i \partial_t \vp^\hbar   =  -  \frac{1}{2}\Delta \vp^\hbar  + (V*\rho_F) \vp^\hbar  \ .
\end{cases}
\end{align}
where $ρ_B(t,x) =  |^ħ(t,x)|^2$
$ρ_F (t,x) = M^-1 ω^ħ(t;x,x)$
are the bosonic and fermionic position densities,
with some prescribed initial data
$(ω_0^ħ, _0^ħ) ∈ℒ ^1(L^2(^d))×L^2(^d)$.
We always assume that $ω_0^ħ= M = ħ^-d$
and $ (ω_0^ħ)^* = ω_0^ħ≥0$;
these properties are of course preserved by the flow generated
in \eqref{hartree hartree 2}.

\vspace{2mm}

In order to analyze the macroscopic limit,
given  $ħ>0$ we denote the   \textit{Wigner transform} of a
reduced density matrix $ω∈ℒ^1(L^2(^d))$ by
\begin{equation}
W^\hbar [\omega] (x,p)
\equiv
\frac{ \hbar^d }{(2 \pi)^{ d}}
\int_{\R^d }
\omega
\Big(
x +    \frac{\hbar }{2 } \xi  , x -   \frac{\hbar }{2 } \xi
\Big)
e^{		 - i  \xi \cdot p	}		 \d \xi    \ ,
\qquad    ( x,p )\in \R^{2d} \ .
\end{equation}
We understand  the above map as a unitary transformation
$W^ħ: L^2(^d_x ×^d_x') →L^2 ( ^d_x ×^d_p )$.
Additionally, we record here the associated inverse map.
Namely, for a regular enough phase-space $f = f(x,p)$ distribution,
we consider the operator on $L^2(^d)$
whose kernel is defined as
\begin{equation}
\t{Op}_w^\hbar[f] (x,x' )
\equiv
\hbar^{-d}
\int_{\R^d}
f\Big(     \frac{x+x' }{2} , p         \Big)
e^{ i p \cdot (x-x')/ \hbar}  \d p  \ ,
\qquad    ( x,y)\in \R^{2d}
\end{equation}
to be the \textit{Weyl quantization} map.
%In what follows, we let $(\omega^\hbar , \vp^\hbar)$
%be the solution of the coupled Hartree equation  \eqref{hartree}.
%We denote by $f^\hbar  (t) = W^\hbar [\omega^\hbar (t)]$
%the Wigner transform of the fermion component of such system, at time $t \in \R$.
%On the other hand, we denote by $ ( f(t ) , \vp(t)) $
%the solution of the Vlasov-Hartee system, \eqref{vlasov hartree}.
%We wish to compare $f^\hbar (t)$ and $f(t)$ as $\hbar \downarrow 0 $.

\vspace{2mm}
Intuitively, two phenomena happen when the $ħ↓0 $ limit is taken.
First, the \textit{dynamics} (i.e. the time evolution of the system)      is changing.
Secondly, the \textit{initial data} of the systems is converging from the quantum-mechanical, to the classical regime.
We shall study each process separately.
%In order to do this, we introduce an intermediate dynamics:
%it corresponds to the Hartree dynamics, but with the quantization of the classical initial data.
%Namely, we consider here  the solution to the equation
%	\label{hartree intermediate}
%	\begin{cases}
%		& i \hbar \partial_t \widetilde \omega^\hbar
%		=
%		[   - (\hbar^2 / 2  )
%		\Delta
%		+   \,
%		(V  * \widetilde \rho_B )  ,
%		\omega
%		]			\\
%		&   i   \partial_t   \widetilde \vp
%		=
%		-  (1/2)\Delta  \widetilde \vp
%		+     \,
%		(V  *   \widetilde \rho_F)
%	\widetilde 	\vp  	\\
%		& ( \widetilde \omega , \widetilde \vp)(0)  = ( \widetilde \omega_0^\hbar , \vp_0)
%	\end{cases}
%where $\widetilde \omega_0  \equiv \t{Op}_w^\hbar [f_0]$ is the Weyl quantization    of the limit $f_0 \in L^1 $, introduced by Assumption \ref{assumption 2}.
%Note as well that the initial data for the boson field does not depend on $\hbar$ anymore.

\subsection{Stability of Hartree-Hartree}
Here and in the rest of the article, we shall be using the following notation
for a special class of semi-classical operators on $L^2(^d)$.
\begin{equation}
\label{semiclassical operator}
\O_{\xi , \eta} \equiv \exp
\big(
i   \xi \cdot x
i 	  \eta \cdot  p
\big)
,  \qquad \xi , \eta \in \R^d
\end{equation}
where $p = - i ħ∇_x$.
If $ω∈ℒ^1(  L^1(^d)  )$ is a trace-class operator,
then it is possible to define the additional (weaker) norm, which we shall use to measure distances
\begin{equation}
||| \omega  |||_s
\equiv
\sup_{\xi, \eta \in \R^d }
(1 + |\xi |  + |\eta|)^{-s }
| 	\t{Tr} (  \O_{\xi , \eta} \omega )|
\ , \qquad s  \geq 0 \ .
\end{equation}
The motivation for the introduction  of such norm comes from the following observation.
If we denote by $ f = W^ħ[ω]$ its Wigner transform,  then it is well-known that
\begin{equation}
\hat{f}(\xi , \eta ) =  \frac{1}{ M }	\t{Tr} (  \O_{\xi , \eta} \omega ) \ , \qquad \forall \xi, \eta \in \R^d
\end{equation}
and, consequently,
\begin{equation}
\label{formula isomorphism}
|f|_s  =\frac{1}{M}  ||| \omega |||_s \ , \qquad \forall s \geq  0 \  ;
\end{equation}
we remind the reader that the norm $ | ·|_s $ has been defined in   \eqref{fourier based norm definition}.
It is important to note that  if $ |ω| ≤N$, then  $|f|_s ≤1 $.

\vspace{3mm}

Our first result towards the prove of Theorem \ref{thm2}
is the stability of the Hartree dynamics
with respect to norm  $  ( |||  ·|||_1 , 	 ·_L^2)$.
Let us note here that this result uses the ideas
presented in [8],
although the result itself is nor stated or proved.

\begin{proposition}[Stablity of Hartree-Hartree]
\label{proposition stability hartree}
Let $(\omega_1^\hbar, \vp_1^\hbar)$ and
$(\omega_2^\hbar , \vp_2^\hbar )$
be two solutions of the coupled Hartree system
\eqref{hartree hartree}
with intial data
$ (\omega_{1,0}^\hbar , \vp_{1,0}^\hbar)$
$ (\omega_{2,0}^\hbar , \vp_{2,0}^\hbar) $,
$\tr \omega_{i,0}^\hbar  =M$
$( \omega_{i,0}^\hbar)^*
\omega_{i,0}^\hbar
\geq0$
for $i=1,2$.
Then, there exists $C>$ such that
\begin{align}
\frac{1}{M}	|||   \omega_1^\hbar  (t) -  \omega_2^\hbar (t) |||_1
\| \vp_1^\hbar (t)  & - \vp_2^\hbar (t)		\|_{L^2}  \\
\nonumber
& 	  \leq
\exp ( \exp C|t|)
\Big(
\frac{1}{M}	|||   \omega_{1,0}^\hbar  -  \omega_{2,0}^\hbar  |||_1
\| \vp_{1 , 0 }^\hbar  - \vp_{2,0}^\hbar 	\|_{L^2}
\Big) \
\end{align}
for all $  t \in \R$
\end{proposition}

\begin{proof}
In what follows, we shall drop the $\hbar$ superscript in order to ease the notation.
Let us  then
consider the unitary evolution groups associated to the solutions
$\{ (\omega_i , \vp_i)\}_{i=1,2}$.
That is, for $i=1,2$,
in the notation of Appendix \ref{section analysis PDEs},
we consider
\begin{equation}
\omega_i(t)
U_{F,i}^* (t) \omega_{i.0}
U_{F,i} (t)
\qquad
\t{and}
\qquad
\vp_i(t) = U_{B,i}(t)
\vp_{i,0}
\end{equation}
where $t\in\R$.
Then, a straightforward computation
using the generators of $U_{F,1}$
$U_{F,2}$ yields
\begin{align}
i \hbar \frac{d}{dt}
\Big(
\omega_{1} (t) -  \omega_{2} (t)
\Big)
U_{F,1}^*  (t)
\Big[
V* \big(
\rho_{B,1} (t)  -  \rho_{B,2} (t)
\big)
, \omega_2(t)
\Big]
U_{F,1}(t)  \ .
\end{align}
Here, $\rho_{B,i} (t,x) = | \vp_i(t,x)|^2$.
Thus, take a semi-classical observable $\O_{\xi , \eta}$ and compute the trace
\begin{align}
\nonumber
\tr \O_{\xi , \eta}
\Big(
\omega_1 (t) -  & \omega_2 (t)
\Big) \\
& =
\tr \O_{\xi , \eta}
\big(
\omega_{1,0} -  \omega_{2,0}
\big) \\
\nonumber
&  \qquad
\frac{ - i }{\hbar}
\int_0^t
\tr
\Big(
\O_{\xi , \eta}
\big[
V* \big(
\rho_{B,1} (s)  -  \rho_{B,2} (s)
\big)
, \omega_2(s)
\big]
\Big)
\d s  \ .
\end{align}
| \tr \O_{\xi , \eta}
\omega_{1,0}-  \omega_{2,0}
) |
\leq
(1 + | \xi|   +  | \eta|  )  \   |||    \omega_1 (0) - \omega_2 (0) |||_1
For the second term,
we use ciclicity of the trace
$\tr AB = \tr BA$,
combined with
the following  general
Fourier series expansion
\begin{equation}
[ V * \rho , \omega ]
\int_{\R^d}
\hat V (k)
\hat\rho(k)
[e_k , \omega] \d k
\qquad
\t{where}
\qquad
\equiv (2\pi)^{-d/2}
e^{i k \cdot x} \ .
\end{equation}
We obtain
\begin{align}
\nonumber
| 	\tr
\Big(
& 	\O_{\xi , \eta}
\big[
V* \big(
\rho_{B,1} (s)  -  \rho_{B,2} (s)
\big)
, \omega_2(s)
\big]
\Big)
|  \\
\nonumber
& \leq
\int_{	\R^d	}  | \hat V (k)| \
|  \hat  \rho_{B,1} (s,k) - \hat \rho_{B,2} (s,k)    |
\  |	 \tr  \Big(
\big[   	e_k ,
\O_{\xi , \eta}
U_{F,1} (t-s)  \big]
\Big)	 \omega_2 (s) 			| \  d k  \\
\nonumber
& \leq
\hbar M  (1 + |\xi| + |\eta| )
\exp (C |t-s |)
\int_{	\R^d	} |k | |\hat V (k) | \
|   \hat \rho_{B,1 }(s,k) -  \hat \rho_{B,2} (s,k)    | \ d k  \ .  \\
& \leq \hbar M  (1 + |\xi| + |\eta| )
\exp (C |t-s |)
\|	 |k| \hat V	\|_{L^1}
\|	\hat \rho_{B,1} (s) - \hat\rho_{B,2 }(s)	\|_{L^\infty  } \ .
\end{align}
In the third  line, we have used the commutator estimate  from Lemma \ref{lemma commutator 1},
$U_\rho  =  U_{F,1}$.
Next, we use
\|	\hat \rho_{B,1}  - \hat\rho_{B,2 }	\|_{L^\infty  }
\leq
\|	  \rho_{B,1}  - \rho_{B,2 }	\|_{L^1 }  \leq  2 \|	 \vp_1 - \vp_2\|_{L^2} $ and analyze the boson fields.
A similar argument shows that
\begin{align}
\nonumber
\|	 \vp_1 (s) - \vp_2(s)	\|_{L^2}
& 	\leq
\|	 \vp_1 ( 0 ) - \vp_2( 0 )	\|_{L^2}
\int_0^s
\|    V*(  \rho_{F,1} (r) -  \rho_{F,2} (r)  )    \vp_2 (r)		\|_{L^2} d r   \\
\nonumber
& \leq
\|	 \vp_1 ( 0 ) - \vp_2( 0 )	\|_{L^2}
\|    |k| \hat V  \|_{L^1}
\int_0^s \d r   \sup_{ k \in \R^d }
\big|
\hat  \rho_{F,1} (r,k ) -  \hat \rho_{F,2} (r,k)
\big| 	   \\
& \leq
\|	 \vp_1 ( 0 ) - \vp_2( 0 )	\|_{L^2}
\frac{ 	\|    |k| \hat V  \|_{L^1}	}{M}
\int_0^s
|||  \omega_{1} (r)   - \omega_2 (r)    |||_{1 } \d r
\end{align}
where we used $\hat \rho_F (k) = M^{-1 } \tr \O_{k,0} \omega. $
We can now close the inequalities and apply the Gronwall inequality.
\end{proof}

\subsubsection{Borrowed Commutator Estimates}
In this subsubsection, we state a result that we used in the proof of the above stability estimate.
They concern   the propagation-in-time
of a certain class of commutator estimates.
These results were originally proved in
\cite[Lemma 4.2]{BenedikterPortaSaffirioSchlein2016}
an interacting system of fermions in the  combined semi-classical
and mean-field regime.
Since the proofs can be easily adapted to the present case, we shall omit them.

\vspace{2mm}

In order to state    it, we introduce the following unitary dynamics.
given a time-dependent position density
$ρ: →L^1(^d) $ satisfying
$ρ(t,x )≥0$
and $∫ρ(t,x) dx =1 $
we consider
\begin{equation}
\label{fermion flow}
\begin{cases}
&  i \hbar \partial_t
(t,s )
\Big(
- \frac{\hbar^2}{2} \Delta
\Big)
U_\rho  (t,s) \\
U_\rho   (t,t) = \1  \ .
\end{cases}
\end{equation}
In the interacting fermion system, one chooses
$ρ(t,x) = N^-1  ω(t ; x,x)$
whereas in our case,
the density corresponds to that of bosons
$ρ(t,x) = | (t,x)|^2$.
Since the estimates are quite robust with respect to the choice of density,
we state here the result for arbitrary $ρ$.  One only needs
$sup_t  |ρ̂(t )|_L^∞< ∞$.

\begin{lemma}
\label{lemma commutator 1}
Assume that
$\int_{\R^d}  \< \xi\>^2 |\hat V (\xi) |\d \xi  < \infty$
and let $ U_\rho (t,s) $
be the unitary flow defined in \eqref{fermion flow}.
Then, there exists $C>0$ such that
for all $(\xi , \eta ) \in \R^d \times \R^d $
and $t,s \in \R $
there holds
\begin{equation}
\sup\Biggl\{
\frac{1}{ | k | }
\Big|
\mathrm{Tr}
\big[
e^{ i k \cdot x },
U_\rho^* (t,s)
\O_{\xi , \eta}
U_\rho  (t,s)
\big]
\omega
\Big|
\ : \
k \in \R^d \ ,
\mathrm{Tr}|\omega| \leq 1
\Biggr\}
\leq
\hbar ( | \xi| + |\eta |)
e^{ C  |t -s |   }
\end{equation}
where $ \O_{\xi , \eta}$ is the semi-classical observable   defined in \eqref{semiclassical operator}.
\end{lemma}

%We will use the above stability result to regularize the Hartree-Hartree dynamics.
%Namely, we introduce
%	\label{regularization}
%	f_{0 }^{(\hbar )} \equiv f_0 * \mathscr{G}_\hbar
%where for $(x,p)\in\R^{2d}$ we let
%	\mathscr{G}_\hbar (x,p) = \hbar^{-2d} \mathscr{G}_1 ( \hbar (x,p)  )
%	\quad
%	\t{and}
%	\quad
%	\mathscr{G}_1  (x,p)  = \pi^{-d/4} \exp(  -x^2 - p^2  )
%as the approximation of the identity,
%with Gaussian profile.
%In particular, under Assumption \ref{assumption 2},
%$f_{0}^{ (\hbar ) }$
%is an   of $f_0$, meaning that
%	 | 	 f_0 -   f_{0}^{(\hbar)}  |_1   \leq C \hbar
%for some $C>0$.
%The fact that we have $L^1$ convergence will be crucial.
%Indeed, in this case, the Vlasov-Hartee equation can be regularized $ |\cdot|_s$ for $s \geq 1. $
%\vspace{1mm }
%In the upcoming proof, we take
%$\omega_{1,0}  =  Op_w^\hbar[f_0^\hbar]$
%$\omega_{2,0} =  Op_w^\hbar[ f_0*\mathscr{G}_\hbar  ] $.
%We then show that under Assumption 2, $f_0^\hbar$ and $f_0*\mathscr{G}_\hbar $
%are close with respect to the $ | \cdot  |_1$ norm.

% \subsection{Stability of Vlasov-Hartee}
% The next step towards the proof of Theorem \ref{thm2},
% is to study the stability of the Vlasov-Hartee system with respect to an appropiate metric.
% We will prove stability with respect to the first Wasserstein distance, defined in our context as follows
% \begin{equation}
% 	d_1 (f_1, f_2)
% 	\equiv
% 	\sup_{ h : |h|_{Lip} \leq 1 }
% 	\int_{	\R^{2d}	} h(x,p)\Big(
% 	f_1 (x,p) - f_2(x,p)
% 	\Big)  \ d x d p
% \end{equation}
% where $f_1 , f_2  \in L_+^ 1 (  \R^{2d}_{x,p })$  have finite first moments.
% The   argument we present is standard, but we include it in order to incorporate the boson field.
%\begin{lemma}[Stability of Vlasov-Hartee]
%	\label{lemma regularization Hartree}
%	Let
%	$ (f_1 , \vp_1 )$
%$ (f_2, \vp_2)$
%be two solutions of the Vlasov-Hartee equation
%with initial data
%$(f_{1,0} , \vp_{1,0})$
%$(f_{2,0} , \vp_{2,0})$.
%Then,  there is $C>0$ such that
% 	\begin{align}
%d_1 \Big(	 f_1 (t), f_2(t) 		\Big)
% 	+
% 	\|      \vp_1 (t)  & - \vp_2(t)		\|_{L^2}  \\
% 	\nonumber
% 	& 	  \leq
% 	C
% 	\exp (   C|t|)
% 	\Big(
%d_1 \Big(	 f_{1,0}  ,  f_{2,0 }	\Big)
% 	+
% 	\| \vp_{1 0 } - \vp_{2,0}	\|_{L^2}
% 	\Big)
% \end{align}
%for all $t \in \R$.
% We denote by $\Phi_{	 \vp } (t,x,p)$ the solution map of the ODE associated
% to the time dependent force field $F_\vp (t) \equiv - \nabla V * |\vp(t)|^2$.
% Then, we may write for the difference of the fermion densities
% \begin{align}
% 	f_1 (t,x,p) - f_2 (t,x,p ) =
% 	f_{1,0} \Big(	 \Phi_{\vp_1(t)}^{ -1 } (t,x,p ) 	\Big)
% 	-
% 	f_{2,0} \Big(	 \Phi_{\vp_2(t)}^{-1 } (t,x,p	) \Big) \ .
% \end{align}
%Thus, for a test function $h : \R^{2d} \rightarrow \C$ we find
%	\<  h ,   f_1  (t) - f_2 (t)\>
%	=
%	\< h \circ \Phi_{ \vp_1 (t)} , f_{1,0}    -  f_{2,0} \>
%	+
%		\< h \circ \Phi_{ \vp_1 (t)}      - h \circ \Phi_{ \vp_2 (t)}    ,    f_{2,0} \>  \ .
%For the first term, we use the fact that
%$h \circ \Phi_{ \vp_1 (t)} $ is a Lipschitz function with
%$ |  h \circ \Phi_{ \vp_1 (t)}   |_{Lip}
% \leq
%  | h |_{Lip} | \nabla \Phi_{\vp_1 (t)}|_{L^\infty }
%  \leq
%  C |t| \,  | h |_{Lip} $.
%On the other hand, for the second term, we use the fact that $\Phi$ is stable with respect to perturbations of the boson field. We find for these terms
%	|	\< h \circ \Phi_{ \vp_1 (t)} , f_{1,0}    -  f_{2,0} \>|  & \leq C | h|_{Lip}  |t | d_1 (f_{1,0} , f_{2,0}) \\
%	|		\< h \circ \Phi_{ \vp_1 (t)}      - h \circ \Phi_{ \vp_2 (t)}    ,    f_{2,0} \>    |
%	& \leq |h|_{Lip} \| \vp_1 (t) - \vp_2(t) \|_{L^2}  \|	 f_{2,0}  \|_{L^1} \ .
%Thus, we now examine the boson field. Duhamel's formula gives
%	\|	\vp_1 (t) - \vp_2(t)	\|_{L^2}
%	\leq
%		\|	\vp_1 ( 0 ) - \vp_2( 0 )	\|_{L^2}
%		+
%	 \int_0^{|t|} d s
%	 \|	 V* (   \rho_{1} (s) - \rho_2(s)   )	\|_{L^\infty} \|	 \vp_{2} (s)	\|_{L^2 } \ .
%In the last line, we can use the fact that
%the map $h_x (y,p ) = V(x-y)$ is Lipshchitz continuous with $|h_x|_{Lip} \leq |V|_{Lip}$.
%Hence, we can now close the inequalities and use the Gronwall inequality.
%%In our proof of Theorem \ref{thm2}, we shall use
%%$f_{1,0}  = f_0$
%%$f_{2,0}  = f_0 * \mathscr{G}_\hbar$.
%%These are close together with respect to the Wasserstein metric, in view of Assumption \ref{assumption 2}.
%%Note that the above is \textit{not} a stability estimate when compared to the one in the last section.
%%Here instead we lose regularity. However, it is enough for us to regularize the Vlasov-Hartee equation
%%in the $|\cdot|_1$ norm.
%%On the other hand, the Hartree equation is easier to regularize,
%%in view of Lemma \ref{lemma stability hartree},
%%  the formula   \eqref{formula isomorphism} and $  |f|_1 \leq \|	f\|_{L^1}$.
%%  We state it here as a Corollary.
%%  \begin{corollary}[Regularization of Hartree-Hartree]
%%  	Let $(f_0 , \vp_0)$ be as in Assumption \ref{assumption 2},
%%  	and let $f_0^{(\hbar)}$ be its regularization, introduced in
%%  	\eqref{regularization}.
%%  	Let  $  ( \widetilde \omega(t), \widetilde \vp (t)) $ and
%%  	$ ( \widetilde \omega^{ (\hbar) } (t)  , \widetilde \vp^{(\hbar)} (t)   )  $ solve the Hartree dynamics, with initial data
%%  	$ ( \t{Op}_w^\hbar  (f_0) , \vp_0) $
%%  	and
%%  	$ ( \t{Op}_w^\hbar  (f_0^{(\hbar)}) , \vp_0) $, respectively.
%%  	Then
%%  	\begin{equation}
%%\frac{1}{N}  			|||   \widetilde \omega (t) -   \widetilde \omega^{ (\hbar) } (t) |||_1
%%  		+
%%  		\|  \widetilde \vp (t)  -  \widetilde \vp ^{ (\hbar)  } (t)		\|_{L^2}
%%  		\leq C \hbar \exp ( C \exp C t )
%%  	\end{equation}
%%for all $ t \in \R \ . $
%%  \end{corollary}
%\subsection{Comparison of dynamics}
%Finally, the heart of the matter.
%Now, we compare the Hartree-Hartree and Vlasov-Hartee dynamics,
%when they are both evolving  regular initial data.
%The last two stability results will justify such regularization process.
%This has the advantage of having extra regularity at our disposal, and allows to show that the following bound holds true
%	\t{Tr}   |   \t{Op}_w (   f^{(\hbar) }(t)  ) |
%		\leq C N
%		\sum_{k=0}^4
%		 \hbar^k
%		 \|	   f^{  (\hbar)  } (t) 	\|_{H_4^k}
%where $ \| f  \|_{H_a^s}^2 \equiv \sum_{|\beta| \leq s }  \int  \< x,p\>^{2a} |\nabla^\beta f(x,p)|^2 dx dp .$
%It is \textbf{possible} to propagate Sobolev regularity along the solutions of the Vlasov-Hartee equation
%and conclude that
%	\t{Tr}   |   \t{Op}_w (   f^{(\hbar) }(t)  ) |
%	\leq C
%	\exp(C_0 t ) N
%	\sum_{k=0}^4
%	\hbar^k
%	\|	   f^{  (\hbar)  }_0 	\|_{H_4^k}
%Here, we assume extra decay of $\hat V$  so that the constant $C_0$
%in the exponent \textbf{does not} depend on any Sobolev norm of $f_0$--this is contrary to BPSS, where they forgot about this extra dependence.
%Thus,  using   $f_0^{(\hbar)  } = f_0 * \mathscr{G}_\hbar $
%we expect that
%$ \|        f^{  (\hbar)  }_0 	\|_{H_4^k} \leq   C \hbar^{-k}$  for $k=0,1,2,3,4$,
%and thus the above equation  becomes bounded above by    $  \lesssim \exp(C_0 t )  N $.
% In particular, under Assumption \ref{assumption 2}
%   the following chain of inequalities hold true   (say, $k=1$)
% \begin{align}
% 	\|	 f_0^{(\hbar)}		\|_{H_m^1}
%%  & =
%%  \int_{	\R^{2d}	}
%%  |   \nabla  (f_0  * \mathscr{G}_\hbar)     |^2 (z)
%%    \<  z \>^{2m}   \d z        \\
%%    & =
%%      \int_{	\R^{2d}	}
%%    |     (f_0  *   \nabla   \mathscr{G}_\hbar     |^2 (z)
%%    \<  z \>^{2m}   \d z        \\
%    & =    \|	 f_0 *    \nabla   \mathscr{G}_\hbar  	\|_{L^2 ( \< z\>^m d z )} \\
%    & \leq
%    \|	 f_0	\|_{L^2 ( \< z\>^m d z )}
%    \|	  \nabla   \mathscr{G}_\hbar  	\|_{L^1 ( \< z\>^m d z )} \\
%    & =
% \hbar         \|	 f_0	\|_{L^2 ( \< z\>^m d z )}
% \|	    \nabla \mathscr{G}_1	\|_{L^1 ( \<  \hbar z\>^m d z )} \\
% & \leq
%C   \hbar         \|	 f_0	\|_{L^2 ( \< z\>^m d z )}  \
% \end{align}
%and very similarly for higher-order derivatives.
%\begin{lemma}[Comparison of the regularized dynamics]
%	Let $\omega^{(\hbar)}(t) $
%	and $f^{(\hbar)}(t) $ solve the Hartree   and Vlasov-Hartee dynamics, respectively,
%	with initial data
%	$\omega^{(\hbar)  }(0 ) = \t{Op}_w (f_0^{(\hbar)})$
%	and
%	$f_0^{(\hbar)}$.
%		Then,  there is $C>0$ such that
%		  	\begin{equation}
%		\frac{1}{N}
%			|||     \omega^{(\hbar)} (t) -
%			 \t{Op}_w (f^{(\hbar)}(t))   |||_2
%			+
%			\|    \vp_1 (t)  -    \vp_2   (t)		\|_{L^2}
%			\leq C \hbar \exp ( C \exp C t )
%		\end{equation}
%	Same as BPSS. We may need to add a number of derivatives to $V$ to have everything corrected.
%\subsection{Proof of the main result}
%	[Proof of Theorem \ref{thm2}]
%	Put the three above lemmas together, with the right choice of initial data.
%\subsection{Borrowed Commutator Estimates}
%In this subsection, we state some results that we used in the proof of Theorem \ref{thm2}.
%They concern   the propagation-in-time
%of a certain class of commutator estimates.
%These results were originally proved in BPSS for
%an interaction system of fermions in the  combined semi-classical
%and mean-field regime.
%Since the proofs can be easily adapted to the present case, we shall omit them.
%In order to state these results, we introduce the following unitary dynamics.
%given a time-dependent position density
%$\rho : \R \rightarrow  L^1(\R^d) $ satisfying
% $\rho\geq0$
%and $\int\rho =1 $
%we consider
% \begin{equation}
%	\label{fermion flow}
%	\begin{cases}
%		&  i \hbar \partial_t \U_\rho  (t,s )
%		=
%		\Big(
%		- \frac{\hbar^2}{2} \Delta
%		+
%		V*\rho(t)
%		\Big)
%		\U_\rho  (t,s) \\
%		& \U_\rho   (t,t) = \1  \ .
%	\end{cases}
%In the interacting fermion system, one chooses
%$\rho(t,x) = N^{-1 } \omega (t ; x,x)$
%whereas in our case,
%the density corresponds to that of bosons
%$\rho(t,x) = | \vp(t,x)|^2$.
%Since the estimates are quite robust with respect to the choice of density,
%we state here the result for arbitrary $\rho$.  One only needs
%$\sup_{t}  |\hat \rho (t )|_{L^\infty}< \infty$.
%\begin{lemma}[BPSS2016 Lemma 4.2]
%	\label{lemma commutator 1}
%	Assume that
%	$\int_{\R^d}  \< \xi\>^2 |\hat V (\xi) |\d \xi  < \infty$
%	and let $ \U_\rho (t,s) $
%	be the unitary flow defined in \eqref{fermion flow}.
%	Then, there exists $C>0$ such that
%	for all $(\xi , \eta ) \in \R^3 \times \R^3 $
%	and $t,s \in \R $
%	there holds
%	\begin{equation}
%		\sup\Biggl\{
%		\frac{1}{ | k | }
%		\Big|
%		\mathrm{Tr}
%		\big[
%		e^{ i k \cdot x }, \U_\rho^* (t,s)
%		\O_{\xi , \eta}
%		\U_\rho  (t,s)
%		\big]
%		\omega
%		\Big|
%		\ : \
%		k \in \R^d \ ,
%		\mathrm{Tr}|\omega| \leq 1
%		\Biggr\}
%		\leq
%		\hbar ( | \xi| + |\eta |)
%		e^{ C  |t -s |   }
%	\end{equation}
%	where $ \O_{\xi , \eta}$ was defined in \eqref{semiclassical operator}.
%\begin{lemma}[BPSS2016 Lemma 5.1]
%	Assume that
%$\int_{\R^d}  \< \xi\>^3 |\hat V (\xi) |\d \xi  < \infty$
%and let $ \U_\rho (t,s) $
%be the   unitary flow defined in \eqref{fermion flow}.
%	Then,
%	there exists $C>0$ such that
%	for all $(\xi , \eta ) \in \R^3 \times \R^3 $
%	and $t,s \in \R $
%	there holds
%	\begin{align}
%		\sup\Biggl\{
%		\Big|    \,
%		\mathrm{Tr}
%		\big[ x_i ,  \big[ x_ j
%		, \U_\hbar^* (t,s)
%				\O_{\xi , \eta}
%		\U_\rho  (t,s)
%		\big] \big ]
%		\omega
%		\, 	\Big|
%		\ : \
%		i,j = 1 , \ldots, d , \
%		& 	\mathrm{Tr}|\omega| \leq 1
%		\Biggr\} \\
%		& \leq C \hbar^2
%		( | \xi| + |\eta |)^2
%		e^{ C  |t -s |   }			\nonumber
%	\end{align}
%	where $ \O_{\xi , \eta}$ was defined in \eqref{semiclassical operator}.

\subsection{Quantum Optimal Transportation}
A lot of technology has been developed in the last decade in the context of Quantum Optimal Transporation,
by means of the introduction of relevant pseudo-metrics.
In this paper, we apply and adapt some of these results to the study of the problem at hand.

\vspace{2mm}
In this context, one is given a probability measure $f ∈𝒫(^2d)$
describing a classical state, and wishes to compare it to a quantum state $ω$,
belonging to the space
\begin{equation}
\label{P h}
\mathcal P (L^2(\R^{d}))
\equiv
\Bigl\{
\omega \in \mathscr L ^1  (  L^2(\R^d) )
\, : \,
\omega = \omega^* \geq 0
,  \
\tr \omega = M = \hbar^{-d}
\Bigr\} \ .
\end{equation}
Let us immediately note here that our normalization is different than most results on the literature,
where one considers trace-class operators $$
with  $= 1 $.
This is of course only a matter of scaling,
and the passage from one to the other is   given by  $= ħ^d ω$.
In particular, we choose the normalization in \eqref{P h}
because it is more natural for the problem at hand.

\vspace{2mm}
Let us now introduce the concept of a \textit{coupling}, between classical and quantum states
\begin{definition}[Coupling]
Given $f \in \mathscr P (\R^{2d})$
\omega \in
\mathcal P (L^2(\R^{d}))
we say that the operator-valued function
$Q : \R^{2d}_{x,p} \rightarrow \mathscr{B}( L^2(\R^{2d}) )$
is a \textbf{coupling} for $f$ and $\omega$, if the following is satisfied
\begin{enumerate}
\item  For almost every $(x,p) \in \R^{2d}$ there holds $Q(x,p) = Q(x,p)^* \geq 0 $.
\item $\int_{\R^{2d} } Q(x,p) \d x \d p = \omega .$
\item $\tr Q(x,p) = \hbar^{-d } \, f(x,p).$
\end{enumerate}
The set of all couplings between $f $ and $\omega$ is denoted by $\mathcal C (f,\omega). $
\end{definition}

Throughout this section,
we denote by $x̂ : 𝒟 (x̂) ⊂L^2(^d) →L^2(^d)$
the standard multiplication operator on $L^2$,
and similarly for $p̂ = - i ħ∇_x   $ on $H^1(^d).$
We introduce the following cost function,
taking values
in the space of unbounded self-adjoint
operators on $L^2(^d)$
\begin{equation}
c_\hbar (x,p) \equiv \frac{1}{2} (x - \hat x)^2 +  \frac{1}{2} ( p - \hat p)^2  \ , \qquad (x,p)\in \R^{2d}  ,
\end{equation}
initially defined on $𝒮 (^d)$ and then closed in $L^2$.
Let us recall  that  we denote by $𝒫_n(^2d)$
the space of measures with finite $n$ moments.
Similalry, we denote
𝒫 _n (L^2(^d))
the space of
quantum states $ω∈𝒫  (L^2(^d))$
such that
$ [ √(ω)( x̂^2 + p̂ ^2)^n   √(ω)]  < ∞ .  $

\begin{definition}[Quantum Wasserstein]
For all $f \in \mathscr P _2 (\R^d)$
$\omega \in 		\mathcal P _2 (L^2(\R^{d})) $
we define the \textbf{second quantum Wasserstein distance}
as the quantity
\begin{equation}
E_\hbar (f, \omega)
\equiv
\hbar^{d/2}
\inf_{Q \in \mathcal C (f,\omega)}
\Big(
\int_{	\R^{2d}	} \tr Q(x,p) c_\hbar (x,p)
\d x \d p
\Big)^{1/2}
\in [0,\infty] \ .
\end{equation}
\end{definition}

Up to scaling, 		the functional $E_ħ$ has been the object of several studies in recent years.
In particular, it has been proven that it is a natural object to study when comparing the dynamics
of the Hartree and Vlasov equation, for system of interacting fermions.
We adapt the proof of
\cite[Theorem 2.5]{GolsePaul2017}
for the present case of interest.

\begin{theorem}
\label{proposition wasserstein}
Assume $V\in C^{1,1} (\R^d,\R)$.
Let $(\omega_\hbar (t),\vp_\hbar (t))$ solve the Hartree-Hartree equation \eqref{hartree hartree 2} with initial data
$(\omega_0^\hbar  , \vp_0^\hbar)$.
Further, let
$(f (t) , \vp(t))$ solve the Vlasov-Hartee equation, with initial data
$(f_0 , \vp_0) $.
Then, there exists $C = C (V)>0$ such that for all $ t \in \R$ there holds
\begin{equation}
\label{theorem 3 rhs}
E_\hbar  \big(  f(t) , \omega_\hbar (t) \big)
\|	 \vp(t)   -  \vp_\hbar (t)	\|_{L^2}
\leq
C\exp(C t^2 )
\Big(
(f_0 , \omega_0^\hbar )
\|	 \vp_0   -  \vp_0^\hbar	\|_{L^2}
\Big) \ .
\end{equation}
\end{theorem}

Two questions remain.
The first one is: what norms are    bounded above by  $E_ħ(f,  ω_ħ) $?
The second one is: is the right hand side of \eqref{theorem 3 rhs} small when $ħ↓0 $?
Neither of these questions has a trivial answer.
Fortunately, they have already been answered quite recently in the literature.
In order to state it, we must introduce
the following Gaussian mollifer at scale $ħ>0$
\begin{equation}
\mathscr G_{\hbar} (z)
\equiv \hbar^{ -d/4}
\mathscr G_1( \hbar^{-1 /2} z)
\quad
\t{with}
\quad
\mathscr G_1 (z)
\equiv
\pi^{-d/4}
\exp(-z^2/2)
\end{equation}
for $z = (x,p) ∈^2d$.
In our notation and scaling, we record the relevant results in the following lemma.

\begin{lemma}
\label{lemma wasserstein}
The following statements are true.
\begin{enumerate}[leftmargin=*]
\item

\cite[Corollary 1.1]{Lafleche2023}
Let $d \geq 2 $,
$\omega \in \mathcal P_2 (L^2(\R^d))$
$f \in \mathscr P_2 (\R^{2d})$.
Assume that $  \omega \leq 1$.
Then, there exists $C = C(d)$ such that
\begin{equation}
\|	 f -  W_\hbar [\omega ]		\|_{ \dot H ^{-1}   }
\leq
E_\hbar [f , \omega]
C \sqrt \hbar \ .
\end{equation}

\item
\cite[Theorem 2.4]{GolsePaul2017}
For all $ f \in \mathscr{P}_2(\R^{2d})$
there holds
\begin{equation}
E_\hbar [f,  \t{Op}_w^\hbar (f * \mathscr{G}_\hbar )  ] \leq \sqrt{\hbar}.
\end{equation}

%	\item Let $n>2$ and let $f\in \mathscr P _n (\R^{2d})$ and $\omega \in \mathscr L_n^1( L^2(\R^d))$ uniformly in $\hbar $.
%	Then, there exists $C = C (n,f, \omega) $ independent of $\hbar$ such that
%	\begin{equation}
%		E_\hbar [f, \omega]
%		\leq C
%		\|f -  W_\hbar [\omega] * \mathscr{G}_\hbar	\|_{\dot W ^{-1,1}  (\R^{2d}) }^\theta
%		+ (\tr \Delta \omega )^{1/2} \hbar
%	\end{equation}
%where $\theta = \frac{1}{2}(1 - \frac{1}{n-1})$.
\end{enumerate}
\end{lemma}

\begin{remark}[Anti-Wick Quantization]
So far, we have introduced the Weyl quantization map, as
the inverse of the Wigner transform.
As is well-known, this quantization map does not preserve positivity.
On the other hand, the anti-Wick (or, Toeplitz) quantization does, and is defined as follows.
Given a phase-space distribution
$\mu \in \mathscr P(\R^{2d})$,
one defines
on $L^2(\R^{d})$ the bounded, self-adjoint operator
\begin{equation}
\t{Op}^{\hbar}_{\t{aw}}
\equiv
\t{Op}^\hbar_w
(\mu * \mathscr G _\hbar )
\hbar^{-d}
\int_{\R^{2d}}
\ket{      f_{x,p}^\hbar   }
\bra{f_{x,p}^\hbar }
\d \mu(x,p)
\end{equation}
where $f_{x,p}^\hbar(y) \equiv  \hbar^{-d/4} g (  \hbar^{-1/2} (x-y)  ) \exp(i p \cdot y / \hbar)$
is a coherent state at scale $\hbar>0$,
with an $L^2$-normalized Gaussian profile $g $.
Note that preservation of positivity follows immediately from the last formula,
since it is the convex combination of positive operators (in this case, orthogonal projections).
In particular, it follows from the previous Lemma
$E_\hbar \big[      f ,   \t{Op}^{\hbar}_{\t{aw}}
(f)   \big] \leq \hbar^{1/2}$.
\end{remark}

Let us now prove Theorem  \ref{proposition wasserstein}.
Since the   proof is a simple adaptation of that of
\cite[Theorem 2.5]{GolsePaul2017}
in which we replace the self-interacting term,
with the interaction with the boson field $(t)$,
we only provide the sketch of the proof and refer to the original reference for the details.

\begin{proof}[Proof of Theorem \ref{proposition wasserstein}]
First,  we   introduce some notation we shall make use of throughout the proof.
Namely, letting
$( \omega^\hbar (t) , \vp^\hbar (t))$
$(f(t) ,\vp(t))$
be as in the statement of the Theorem, we let
the bosonic densities be
\begin{equation}
\rho_f(t,x) \equiv
|\vp (t,x)|^2
\quad
\t{and}
\quad
\rho_\omega (t,x)
\equiv
|\vp^\hbar  (t,x)|^2
\ , \qquad \forall (t,x) \in \R\times \R^d  \ .
\end{equation}
The notation is made such that they enter the equations for $f$ and $\omega$, respectively.

\vspace{2mm}

\textit{Step 1. Propagation of moments}.
Let $f_0 \in \mathscr{P}_2(\R^{2d})$ be the initial datum of the Vlasov-Hartee system.
Then,  it holds true
that $f(t) \in \mathscr{P}_2 (\R^{2d})$.
Indeed, let $\Phi_\vp$ be the ODE characteristics map
from Theorem \ref{theorem wp 2}.
Then, the claim follows from the formula $f(t,z) = f_0 \circ \Phi_{\vp}^{-1}(t,z)$ and changing variables
$z = (x,p)\in \R^{2d}$ in
\begin{equation}
\int_{	\R^{2d}	} |z|^2  f(t,z) \d z
\int_{	\R^{2d}	}   |    \Phi_\vp (t,z) |^2  f_0(z) \d z
\leq
\kappa(t)^2
\int_{	\R^{2d}	}  |z|^2     f_0(z) \d z
< \infty \ .
\end{equation}
Here, we have used that the estimate
| \Phi_\vp (t,z)| \leq  \kappa(t)  |z|,
see Remark \ref{remark grow}.
\vspace{2mm}

\textit{Step 2. Choice of  Coupling}.
Given $Q_0^\hbar \in \mathcal C (f_0 , \omega_0)$
we consider $Q^\hbar(t)$ to be the time-dependent coupling solving the PDE
\begin{equation}
\label{coupling dynamics}
\begin{cases}
& 	\partial_t Q^\hbar
\{
\frac{1}{2}  p^2  + V*\rho_f (t,x)  ,
\}
\frac{i}{\hbar}
\big[
\frac{1}{2}\hat p^2
V* \rho_\omega (t,\hat x)  , Q_\hbar
\big] = 0 \   ,   \\
& Q^\hbar (0)  =  Q^\hbar_0 \ .
\end{cases}
\end{equation}
Here, we denote by $ \{ F,G\}  = \nabla_x F \nabla_p G - \nabla_p F \nabla_x G$ the Poisson bracket of two classical observables,
and we recall that $\hat p = - i \hbar \nabla_x$ and $\hat x$ are
the standard momentum and position observables in $L^2(\R^d)$.
One may define $Q^\hbar (t)$  by means of   conjugation of a unitary map, and composition of a vector field--see for details.
In particular, it follows from  such representation and
\cite[Lemma 4.2]{GolsePaul2017}
that this procedure actually defines a coupling between $f(t)$ and $\omega^\hbar (t)$.
Namely,  $Q_\hbar (t) \in \mathcal{C} (f(t) , \omega^\hbar  (t))$ for all $t\in \R $.

\vspace{2mm}

\textit{Step 3. Dynamics of the coupling.}
We now estimate the growth of the second order moments of the coupling $Q_\hbar (t)$.
Namely, we define
\begin{equation}
\mathcal E _\hbar (t)
\equiv
\hbar^d
\int_{	\R^{2d}	} \tr \big[
c_\hbar (x,p) Q_\hbar (t;x ,p)
\big] \d  x \d  p \  , \qquad \forall t \in \R \
\end{equation}
and compute its time derivative as follows.
In view of \eqref{coupling dynamics}
we find that
\begin{align}
\nonumber
\hbar^{-d}
\frac{d}{dt}
\mathcal E _\hbar (t)
& =
\int_{	\R^{2d}	}
\tr
\Big[
Q_\hbar (t;x,p)
\Big\{
\frac{1}{2} p^2  + V*\rho_f (t,x)  ,
c_\hbar (x,p)
\Big\}
\Big]
\d x \d p  \\
&   +
\frac{i}{\hbar}
\int_{	\R^{2d}	}
\tr
\Big[
Q_\hbar (t;x,p)
\big[
\frac{1}{2}\hat p^2
V* \rho_\omega (t,\hat x)  , c_\hbar (x,p)
\big]
\Big]
\d x \d p
\end{align}
where we have integrated by parts and put both  brackets on the cost function.
[37]
we calculate explicitly the above brackets and find that
\begin{align}
\big\{
\frac{1}{2}p^2 + V*\rho_f (t,x) , c_\hbar (x,p)
\big\}
& 	 =
p \cdot (x - \hat x)
- \nabla V * \rho_f (t,x) (p - \hat p )  \\
\frac{i}{\hbar}	\big\[
\frac{1}{2}\hat p ^2 + V*\rho_\omega (t,\hat x ) , c_\hbar (x,p)
\big\]
& 	 =
- \frac{1}{2} (x - \hat x)  \hat p  -  \frac{1}{2}  \hat p (x - \hat x )  \\
\nonumber
&  +   		  \frac{1}{2} \nabla V*\rho_\omega(t,\hat x) \cdot  ( p - \hat p )
+     \frac{1}{2}		( p - \hat p )	\cdot  \nabla V*\rho_\omega(t,\hat x)    \ .
\end{align}
Thus, a straightforward manipulation,
combined with the   identities
for self-adjoint operators
$AB + BA \leq A^2 + B^2  $
\begin{align}
\nonumber
\hbar^{-d}
\frac{d}{dt}
\mathcal E _\hbar (t)
& =
\int_{	\R^{2d}	}
\tr
\Big[
Q_\hbar (t;x,p)
\frac{1}{2}
\Big(
(x - \hat x) \cdot (p - \hat p)
(p - \hat p)\cdot  (x - \hat x)
\Big)
\Big]
\d x \d p
\\
\nonumber
&   +
\int_{	\R^{2d}	}
\tr
\Big[
Q_\hbar (t;x,p)
\frac{1}{2}
\Big(
\nabla V * \rho_\omega (t,\hat x)
\nabla V  * \rho_f (t,x)
\Big)
\cdot(p-\hat p )
\Big]
\d x \d p \\
&  +
\int_{	\R^{2d}	}
\tr
\Big[
Q_\hbar (t;x,p)
\frac{1}{2}
(p-\hat p )
\cdot
\Big(
\nabla V * \rho_\omega (t,\hat x)
\nabla V  * \rho_f (t,x)
\Big)
\Big]
\d x \d p \\
\nonumber
& \leq
\hbar^{-d}
\mathcal E _\hbar (t)
\frac{1}{2}
\int_{	\R^{2d}	}
\tr  \Big[
Q_\hbar (t;x,p)
|\nabla V * \rho_\omega  (t,\hat x)     -    \nabla V *\rho_f(t,x)  |^2
\Big] \d x \d p  \ .
\end{align}
The second	term in the last displayed  equation  can be estimated as follows.
We use the triangle  inequality and the fact that $Q_\hbar (t)$ is a coupling for $f (t)$
to find that
\begin{align}
\nonumber
\int_{	\R^{2d}	}
\tr  \Big[
&Q_\hbar (t;x,p)
|\nabla V * \rho_\omega  (t,\hat x)     -    \nabla V *\rho_f(t,x)  |^2
\Big] \d x \d p    \\
\nonumber
& 	\leq
\int_{	\R^{2d}	}
\tr  \Big[
Q_\hbar (t;x,p)
|\nabla V * \rho_\omega  (t,\hat x)     -    \nabla V *\rho_\omega (t, x)  |^2
\Big] \d x \d p   \\
\nonumber
& 	+
\int_{	\R^{2d}	}
\tr  \Big[
Q_\hbar (t;x,p)
\Big]
\,
| \nabla V * \rho_\omega  (t,  x)     -    \nabla V *\rho_f(t,x)  |^2
\d x \d p    \\
\nonumber
& \leq
| \nabla V   |_{Lip}^2   \|  \rho_\omega(t)	\|_{L^1}^2  \hbar^{-d } \mathcal E _\hbar (t)
\|	 \nabla V * \rho_\omega  (t)     -    \nabla V *\rho_f(t)  	\|_{L^\infty}^2
\|	f (t)	\|_{L^1 }   \\
& \leq
| \nabla V   |_{Lip}^2   \hbar^{-d }   \mathcal E _\hbar (t)
\|   \nabla V \|_{L^\infty }^2
\|	 \vp(t) -  \vp^\hbar (t)	\|_{L^2}^2
\end{align}
We can now collect the previous estimatesa and   integrate in time
to find that
\begin{equation}
\label{gronwall 1}
\hbar^{-d}
\mathcal{E}_\hbar (t)
\leq
\hbar^{-d}	\mathcal{E}_\hbar (0)
\int_0^t
\Big(
\hbar^{-d}
\mathcal{E}_\hbar (s )
\|	 \vp(s) -  \vp^\hbar (s)	\|_{L^2}^2
\Big) \d s \ .
\end{equation}
It suffices now to estimate the difference in the $L^2$ norm of the boson fields.

\vspace{2mm}

\textit{Step 4. Estimates on the boson densities}.
The equations for the boson fields can be written
in mild formulation as follows
\begin{align}
\vp (t)  & =
e^{  -  i t \Delta /2  }
\vp_0
- i
\int_0^t
e^{    - i (t-s) \Delta /2    }
V*  \widetilde \rho _f (s )   \vp(s) \d s  \\
\vp^\hbar (t)  & =
e^{  -  i t \Delta /2  }
\vp_0^\hbar
- i
\int_0^t
e^{    - i (t-s) \Delta /2    }
V*  \widetilde \rho _\omega (s)   \vp^\hbar (s) \d s
\end{align}
where we denote for the fermion densities
$\widetilde \rho_f (t,x) = \int_{\R^d} f(t;x,p) \ d p $
$\widetilde \rho _\omega (t,x) = (1/N) \omega(t;x,x)$.
We estimate the $L^2$ norms  as follows
\begin{align}
\nonumber
\|	 &  \vp(t)  -  \vp^\hbar (t)     \|_{L^2}
\\
\nonumber
& 	\leq
\|	 \vp_0 -  \vp^\hbar_0    \|_{L^2}
\int_0^t
\|		  (V* \widetilde \rho_f (s)   -   V* \widetilde \rho_\omega (s)      )   \vp(s) 	   \| \d s  	 +
\int_0^t
\|		    V* \widetilde \rho_\omega (s)     (   \vp(s) 	 - \vp^\hbar (s))    \| \d s  \  \\
\nonumber
& \leq
\|	 \vp_0 -  \vp^\hbar_0    \|_{L^2}
\int_0^t \|		V* \big(     \widetilde \rho_f (s) -  \widetilde \rho_\omega (s)  	 		\big)   	 \|_{L^\infty  }  \d s
\|	 V	\|_{L^\infty }
\int_0^t
\|	 \vp(s)  -  \vp^\hbar (s)     \|_{L^2} \d s \ .
\end{align}
Finally, we notice that because $Q_\hbar (t)$ is a coupling between $f(t)$ and $\omega(t)$,
we obtain that  for all $X \in \R^d$
\begin{align}
V* \widetilde \rho_f (s,X)
V* \widetilde \rho_\omega (s,X)
&   =     \int_{	\R^{2d}	}
\tr
\Big[
Q_\hbar (t;x,p)
\big(
\nabla V(X -x )
\nabla V(X - \hat x )
\big)
\Big] \d x \d p
\end{align}
and, consequently, using the Cauchy Schwarz inequality
\footnote{In the following form
\int_{	\R^{2d}	} \d z  \tr ( A(z)^* B(z)  ) \d z
\leq
(\int_{	\R^{2d}	}   \tr (A^*(z)A(z))  dz     )^{1/2}
(\int_{	\R^{2d}	}   \tr (B^*(z)B(z))  dz     )^{1/2}
we find
\begin{align}
\nonumber
| V* \widetilde \rho_f (s,X)
V* \widetilde \rho_\omega (s,X)  |^2
& \leq
\Big| 	   \int_{	\R^{2d}	}
\tr
\Big[
Q_\hbar (t;x,p)
\big|
\nabla V(X -x )
\nabla V(X - \hat x )
\big|^2
\Big] \d x \d p  \Big| 		 \\
\nonumber
& \leq
\|	 \nabla V	\|_{Lip}^2
\int_{	\R^{2d}	}
\tr
\Big[
Q_\hbar (t;x,p)  |x - \hat x |^2
\Big] \d x \d p \\
& 	   \leq  2
\|	 \nabla V	\|_{Lip}^2
\hbar^{-d }
\mathcal E_\hbar (t) \ .
\end{align}

We conclude that for some $C = C(V)$ there holds
\begin{equation}
\label{gronwall 2.0}
\|	   \vp(t)  -  \vp^\hbar (t)     \|_{L^2}
\leq
\|	 \vp_0 -  \vp^\hbar_0    \|_{L^2}
\int_0^t
\Big(
\hbar^{ - \frac{d}{2}}				\mathcal E _\hbar (t)^{\frac{1}{2}} +
\|	   \vp(s)  -  \vp^\hbar (s)     \|_{L^2}
\Big) \d s \  .
\end{equation}
In order to compare it with the inequality found in
\eqref{gronwall 1}, we take the square of both sides in \eqref{gronwall 2.0}
and use the Cauchy-Schwarz inequality to find that for some $C = C(V)>0$ there holds
\begin{equation}
\label{gronwall 2}
\|	   \vp(t)  -  \vp^\hbar (t)     \|_{L^2}^2
\leq
\Bigg(
\|	 \vp_0 -  \vp^\hbar_0    \|_{L^2}
\int_0^t
\Big(
\hbar^{-d }			\mathcal E _\hbar (t)
\|	   \vp(s)  -  \vp^\hbar (s)     \|_{L^2}^2
\Big) \d s
\Bigg)
\  .
\end{equation}

\textit{Step 5. Conclusion.}
In order to conclude our argument, we  put the inequalities \eqref{gronwall 1} and \eqref{gronwall 2} together
to find that
\begin{align}
\nonumber
\hbar^{-d }				\mathcal E _\hbar (t)
\|	   \vp(t)  -  \vp^\hbar (t)     \|_{L^2}^2
& 	\leq
\Big(
\hbar^{-d }				\mathcal E _\hbar ( 0 )
\|	   \vp_0  -  \vp^\hbar_0   \|_{L^2}^2
\Big)  \\
& 				+
C (1+t )
\int_0^t
\Big(
\hbar^{-d }				\mathcal E _\hbar (s)
\|	   \vp(s)  -  \vp^\hbar (s)     \|_{L^2}^2
\Big) \d s \  .
\end{align}
Thus, we may  apply the Gr\"onwall inequality to find that
there exists $C>0$ such that
\begin{equation}
\hbar^{-d }				\mathcal E _\hbar (t)
\|	   \vp(t)  -  \vp^\hbar (t)     \|_{L^2}^2
\leq
C \exp(Ct^2 )
\Big(
\hbar^{-d }				\mathcal E _\hbar ( 0 )
\|	   \vp_0  -  \vp^\hbar_0   \|_{L^2}^2
\Big)  \ .
\end{equation}
In order to conclude,
let us recall that
$   				E_\hbar (f(t)  ,  \omega (t))^2
\leq
\hbar^d 			\mathcal E _\hbar (t)  .
we recall  that $	 		\mathcal E_{\hbar }(0) $ is defined in terms of the initial datum
$Q^\hbar (0) = Q_0^\hbar \in \mathcal C (f_0 , \omega_0^\hbar)$, corresponding to an arbitrary coupling.
we  minimize the right hand side over all couplings,
and take the square root,  to finally find that
\begin{equation}
E_\hbar (f(t)  ,  \omega (t))
\|	   \vp(t)  -  \vp^\hbar (t)     \|_{L^2}
\leq
C \exp(C
\Big(
(f_0 , \omega_0^\hbar )
\|	 \vp_0 -  \vp^\hbar_0    \|_{L^2}
\Big)  \ .
\end{equation}
This finishes the proof.
\end{proof}

\subsection{Proof of Theorem \ref{thm2}}
In this subsection, we combine the results  previously established in
Proposition \ref{proposition stability hartree} and
Theorem \ref{proposition wasserstein}.

\begin{proof}[Proof of Theorem \ref{thm2}]
Let us consider $(\omega_0^\hbar, \vp_0^\hbar)$,
$f_0^\hbar \equiv W^\hbar [\omega_0^\hbar]$,
and $(f_0, \vp_0)$
as in the statement of Theorem \ref{thm2}.
The proof is divided into two steps.
In the first step, we consider the evolution given by the  Hartree-Hartree dynamics,
and use Proposition \ref{proposition stability hartree}
to change the initial data from $\omega_0^\hbar = \t{Op}_w^\hbar [ f_0^\hbar ]$
to the intermediate initial data given by
\begin{equation}
\label{intermediate}
\widetilde	\omega_0^\hbar
\equiv
\t{Op}_w^\hbar [ f_0 * \mathscr{G}_\hbar] \ .
\end{equation}
Here, $\mathscr{G}_\hbar (z )= \hbar^{-d/4}   \mathscr{G}_1 (\hbar^{-1} z ) $
is the Gaussian mollifier at scale $\hbar$.
In the second step, we use Proposition \ref{proposition wasserstein}
and Lemma  \ref{lemma wasserstein}
to go from the Hartree-Hartree dynamics, to the Vlasov-Hartee dynamics.
Let us note that these steps will involve different metrics when measuring the distance of the fermion density.
In order to conclude, we put these distances together
by restricting  our collection of test functions.

\vspace{2mm}
\textit{Step 1.}
Let $( \omega^\hbar (t) , \vp^\hbar(t))$ solve the Hartree-Hartree dynamics with initial data $(\omega_0^\hbar, \vp_0^\hbar)$,
and let $( \widetilde \omega ^\hbar  (t),  \widetilde \vp^\hbar (t) )$
solve the Hartree-Hartree dynamics with initial data
$ ( \widetilde \omega_0^\hbar  ,  \vp_0^\hbar)$, where
the fermion component has been defined in \eqref{intermediate}.
We consider its Wigner transform as
\begin{equation}
\label{tilde f}
\widetilde f^ \hbar (t) \equiv W^\hbar [\widetilde \omega ^\hbar (t)] \ , \qquad \forall t \in \R \ ,
\end{equation}
which we shall refer to as the \textit{intermediate dynamics}.
A direct application of  the stability estimate contained in
Proposition \ref{proposition stability hartree},
together with
the isometric property \eqref{formula isomorphism} yields
\begin{align}
\label{proof thm2 eq 1}
|	   f^\hbar (t)  -  \widetilde f ^\hbar (t) 				|_1
\| \vp^\hbar  (t)   - \widetilde \vp^\hbar (t)		\|_{L^2}
\leq  C \exp ( C \exp C |t|)
|	   f^\hbar_0  -   f_0* \mathscr G_\hbar 			|_1
\end{align}
for a constant $C>0$, and all $t\in \R$.
It suffices to estimate the rigth hand side of Eq. \eqref{proof thm2 eq 1}.
In particular, the triangle inequality gives
\begin{align}
\label{proof thm 2 eq 2}
|	   f^\hbar_0  -   f_0* \mathscr G_\hbar 			|_1
\leq
|	   f^\hbar_0  -   f_0  			|_1
|	   f_0  -   f_0* \mathscr G_\hbar 			|_1    \ .
\end{align}
The first term on the right hand side of Eq. \eqref{proof thm 2 eq 2} is already contained in the estimate of Theorem \ref{thm2},
so it suffices to estimate the second term.
Indeed, we find that
for $\zeta \in \R^{2d}$
\begin{align}
|	\hat f_0 (\zeta) -  \widehat{ f_0 * \mathscr G _\hbar }  (\zeta )		|
|	 1 -  \widehat{\mathscr G_1 } (\hbar \zeta )	|	\, |	 \hat f_0 (\zeta) 	|
\leq  	 \t{Lip}
(  \widehat{ \mathscr G _1}  )				\hbar |\zeta|  \ .
\end{align}
In the last line we have used the fact that
$  \|	 \hat f_0	\|_{L^\infty }    \leq    \|	f_0\|_{L^1} \leq 1$,
and $\widehat{\mathscr{G} _1 }(0)=1$.
Upon taking supremum over $\zeta \in \R^{2d}$
one finds that
$|	f_0 - f_0*\mathscr G _\hbar	|_1 \leq C \hbar$.
Putting everything together, we find
\begin{equation}
\label{proof thm 2 eq 3}
|	   f^\hbar (t)  -  \widetilde f ^\hbar (t) 				|_1
\| \vp^\hbar  (t)   - \widetilde \vp^\hbar (t)		\|_{L^2}
\leq
\exp ( C \exp (C |t|))
\Big(
|	   f^\hbar_0  -   f_0  			|_1
\hbar
\Big)
\end{equation}
for a constant $C>0$ and all $t \in \R$.

\vspace{2mm}
\textit{Step 2.}
Let $(f(t), \vp(t))$ be the solution of the Vlasov-Hartree system with initial data
$(f_0 , \vp_0)$.
Similarly,  let $(	\widetilde \omega^\hbar (t),   \widetilde \vp ^\hbar (t)		)$
be the solution of the Hartree-Hartree syste with
initial data $(\widetilde \omega_0^\hbar , \vp_0^\hbar)$, as defined in  Eq. \eqref{intermediate}.
Then, Proposition \ref{proposition wasserstein} immediately implies that
\begin{equation}
\label{proof thm2 eq 3}
E_\hbar \big(
f(t), \widetilde \omega _\hbar (t)
\big)
\|	  \vp(t) -   \widetilde \vp ^\hbar (t)	\|_{L^2}
\leq
C\exp(C t^2 )
\Big(
(f_0 ,  \widetilde \omega _0^\hbar )
\|	 \vp_0   -  \vp_0^\hbar	\|_{L^2}
\Big) \ .
\end{equation}
Consequently, letting $\widetilde f ^\hbar (t) = W^\hbar [\widetilde \omega^\hbar (t)]$
be the Wigner transform of the intermediate dynamics,
we combine Eq. \eqref{proof thm2 eq 3} with
the two estimates found in  Lemma \ref{lemma wasserstein}
to conclude that
\begin{equation}
\label{proof thm 2 eq 4}
\|	 f(t) - \widetilde f ^\hbar (t)	\|_{\dot H ^{-1}}
\|	  \vp(t) -   \widetilde \vp ^\hbar (t)	\|_{L^2}
\leq
C\exp(C t^2 )
\Big(
\hbar^{1/2}
\|	  \vp_0 -   \vp_0^\hbar 	\|_{L^2}
\Big) \ .
\end{equation}
\vspace{1mm}
\textit{Conclusion.}
First, we estimate the density of the fermions.
Namely, let $f^\hbar(t)$ and $f(t)$ be as in the statement of Theorem \ref{thm2},
and let $\widetilde f^\hbar (t)$
be the intermediate dynamics we have previously introduced.
Further, let us  consider a test function $h : \R^{2d} \rightarrow \C$
such that
$  \|	  \< \eta\> \hat h 	 \|_{L^1}$
$  \|	  \< \eta\> \hat h 	 \|_{L^2}$
are both finite.
Then, the triangle inequality  combined
with Eqs. \eqref{proof thm 2 eq 3} and \eqref{proof thm 2 eq 4}
imply that
\begin{align}
\nonumber
|	 \<   h , ( f(t)   -  f^\hbar (t))  \>		|
& \leq
|	 \langle    h , ( f^\hbar (t)   -   \widetilde f ^\hbar (t))  \rangle 		|
|	 \<   h , ( f(t)   -  \widetilde f ^\hbar (t))  \>		|	 	\\
& \leq
\|	  \<\zeta\>  \hat h  	\|_{L^1}
|	 f^\hbar (t)  - 	  \widetilde f ^\hbar (t)  		|_1
\|	   |\zeta|  \hat h  	\|_{L^2}
\|	 f(t) - \widetilde f ^\hbar (t)	\|_{\dot H ^{-1}  } \\
\nonumber
& \leq
\|	  \<\zeta\>  \hat h  	\|_{L^1}
\Big(
|	   f^\hbar_0  -   f_0  			|_1
\hbar
\Big)
\|	   |\zeta|  \hat h  	\|_{L^2}
\Big(
\hbar^{1/2}
\|	  \vp_0 -   \vp_0^\hbar 	\|_{L^2}
\Big)
\end{align}
$C_1(t) \equiv
C\exp(C t^2 )
$C_2(t) \equiv C \exp ( C \exp C |t|)$.
Similarly, for the boson fields we find
\begin{align}
\nonumber
\| \vp(t)     -  \vp^\hbar (t)		\|_{L^2 }
& 	\leq
\| \vp(t)     -  \widetilde \vp ^\hbar (t)		\|_{L^2 }
\|   \vp^\hbar (t)	 -   \widetilde \vp^\hbar (t)    	\|_{L^2 }\\
&  \leq
\Big(
|	   f^\hbar_0  -   f_0  			|_1
+ \hbar
\Big)
C_1 (t)
\Big(
\hbar^{1/2}
\|	  \vp_0 -   \vp_0^\hbar 	\|_{L^2}
\Big) \ ,
\end{align}
where $C_1(t)$ and $C_2(t)$ are as above. This finishes the proof of the theorem.
\end{proof}

%\noindent \textit{Evolution of the Wigner transform}. First, we compute the time evolution of $f^\hbar (t)$. A straightforward calculation shows that
%	\partial_t f^\hbar
% & 	 =
%	 W^\hbar [\partial_t   \omega^\hbar]   \\
%	 & =   \frac{1}{i\hbar}
%		 W^{\hbar }  \big[       [h_F , \omega^\hbar]      \big] \\
%		 		 	 & =
%\frac{  i  \hbar }{2}
%		 W^{\hbar }  \big[    [        \Delta , \omega^\hbar      ]\big]
%		 -
%		  \frac{i }{\hbar}
%		 W^\hbar \big[		  [    V* \rho_B , \omega^\hbar     ]  			\big] \ .
%		 \label{eq 2}
%It is a standard calculation to show that
%$	 i \hbar/2  W^{\hbar }  \big[    [        \Delta , \omega^\hbar      ]\big]  = -     p \nabla_x f^\hbar$.
%Let us momentarily denote
%$V_B \equiv V * \rho_B$.
%For the commutator
%$[V_B , \omega] (x,x' ) = ( V_B(x)  - V_B(x'))  \omega(x,x') $
%we use the following Taylor expansion
%V_B(x + \hbar \xi /2 )  - V_B( x - \hbar \xi /2 )
% =  \hbar \xi \nabla V_B(x)
%   +
%   \hbar^2
%  r_\hbar ( x , \xi )
%where $r_\hbar (x,\xi)$
%is a second-order reminder term that satisfies,
%for a constant $C>0$
%	|    r_\hbar ( x , \xi )   | \leq C |\xi|^2 \|	 \nabla^2 V_B	\|_{L^\infty}
%	 \qquad
%	 \forall (x,\xi)\in \R^{2d} , \ \forall \hbar >0 \ .
%We use the above decomposition to write for the second term in \eqref{eq 2}
%	\frac{i }{\hbar}
%	W^\hbar \big[		  [    V_B , \omega^\hbar     ]  			\big]
%	=
%	  \nabla V_B  \nabla_p f^\hbar
%	   -
%	  \hbar
%	  R_\hbar[f^\hbar ]
%where the second term in the above equation corresponds to a remainder term given by
%	R_\hbar [f ](x,p)
%	=
%	i \int_{	 \R^d \times \R^d 	}
%	  f(x , p')
%	  e^{ - i \xi (p - p') }
%	  r_\hbar (x,\xi ) 	\d p'   \d \xi
%where in the derivation of the above equation we have used the inverse relation
%$\omega (x,x' ) = N  \int_{	\R^d	} f \big(\frac{x+x'}{2} , p \big) e^{ i p(x-x')/\hbar } \d p $.
%We conclude that the dynamics  for $(f^\hbar , \vp^\hbar)$
%is given by the system
%	\label{vlasov hartree 3}
%	\begin{cases}
%		& 		(\partial_t  + p  \nabla_x + F_B (t,x)   \nabla_p   ) f^\hbar   = \hbar R_\hbar[f^\hbar ]	\\
%		&		  i \partial_t \vp^\hbar
%		=
%		- \frac{1 }{2} \Delta   \vp^\hbar  +   ( V*\rho_F   )   \vp^\hbar 	\\
%		&  (f^\hbar,\vp^\hbar) ( 0 )  = (f_ 0^\hbar , \vp_0^\hbar)   \
%	\end{cases}
%where $F_B  = \nabla V * |\vp^\hbar|^2 $
%$\rho_F (x)= \int_{	\R^d	} f^\hbar (x,p) \d p $
%depend on $\hbar$ but we do not make it explicit in the notation.
%\noindent \textit{Comparison with an intermediate dynamics}.
%Let us now compare $(f^\hbar, \vp^\hbar)$
%with the solution of an intermediate dynamics, in which we have removed the remainder term, but we do not change the initial data.
%Namely, let us consider now
%$(\widetilde f ^\hbar, \widetilde \vp ^\hbar )$
%as the solution of the system
%	\label{vlasov hartree 4}
%	\begin{cases}
%		& 		(\partial_t  + p  \nabla_x + \widetilde F_B (t,x)   \nabla_p   ) \widetilde f ^\hbar    = 0 	\\
%		&		  i \partial_t 	\widetilde \vp ^\hbar
%		=
%		- \frac{1 }{2} \Delta   \widetilde \vp ^\hbar  +   ( V*\widetilde \rho _F   )   \widetilde \vp ^\hbar 	\\
%		&  ( \widetilde f ^\hbar,  \widetilde \vp ^\hbar) ( 0 )  = (f_ 0^\hbar , \vp_0^\hbar)   \
%	\end{cases}
%where $\widetilde F_B  = \nabla V * |\widetilde \vp ^\hbar|^2 $
%$\widetilde \rho _F (x)= \int_{	\R^d	} \widetilde f ^\hbar (x,p) \d p $
%also depend on $\hbar$ but we do not make it explicit in the notation.
%In order to establish the comparison between $f^\hbar$ and $\widetilde f^\hbar$ we proceed as follows.
%First, we introduce an important semigroup.
%Namely, given an arbirary field boson field  $\vp : \R \rightarrow  L^2(\R^3 )$ we introduce the ODE problem
%	\begin{align}
%	\label{ODE 2}
%	\begin{cases}
%		& 			 d/dt  \, x (t) = v(t) 		\\
%		& d/dt \, v(t)	 = - \nabla(V * |\vp(t)|^2 )( x(t) )	\\
%		&  (x,v) (0) =  z_0 \equiv (x_0 , v_0 )\in \R^{2d}.
%	\end{cases}  \ .
%We denote by $\Phi_\vp (t,z) \equiv (x (t) , v(t))$ the associated solution of the ODE, at time $t \in \R$, with intial data $z$.
%One may verify that the map $z\mapsto \Phi_\vp(t,z)$ is a $C^1$ volume preserving diffeomorphism.
%In particular,  the solution map
%induces a norm-preserving map in $L^1 (\R^{2d})$ as follows
%	 S_\vp (t) [f]  (z ) \equiv f  (    \Phi_\vp^{-1} (t,z)   )  \qquad \forall t\in \R , \ \forall z \in \R^{2d} , \ \forall f\in L^1(\R^{2d}) \ .
%In other words, $ \|	 S_\vp (t) [f]	\|_{L^1}  =   \|	f	\|_{L^1}$.
%Of course $(S_\vp (t))_{t \in \R} $ is the semi-group associated to the PDE
%	 (\partial_t + p \nabla_x +  F_\vp  (t)\nabla_p) f (t,x,p)  = 0 \ , \qquad f(0) = f_0
%where $F_\vp(t) \equiv  - \nabla V * |\vp (t)|^2$\ .
%Our next step is to use Duhamel's formula  and the semigroups
%$ S^\hbar \equiv S_{\vp^\hbar }$
%$\widetilde S ^\hbar \equiv S_{\widetilde \vp ^\hbar}$
%to compare $f^\hbar$ and $\widetilde f^\hbar$.
%Namely, we start by writing
%	f^\hbar (t)
%	& =
%	 S^\hbar (t) f_0^\hbar
%	+
%	\hbar
%	\int_0^t S^\hbar (t-s ) R_\hbar   [  f^\hbar (s) ]  \d s     ,\\
%	\widetilde f ^\hbar (t)
%	& =
%	 \widetilde S ^\hbar (t) f_0^\hbar
%	 ,
%and conclude that for any test function $h(x,p)$ there holds
%	\<  h , f^\hbar (t) - \widetilde f ^\hbar (t)\>
%	 &  =
%	 \<  h ,   	 S^\hbar (t) f_0^\hbar  - \widetilde S ^\hbar (t) f_0^\hbar      \>
%	 +
%	 \hbar
%	 	\int_0^t  \< h , S^\hbar (t-s ) R_\hbar   [  f^\hbar (s) ]  \>  \d s   \\
%	& =
%		 \<  h\circ \Phi^\hbar  (t) - h \circ  \widetilde \Phi^\hbar (t) ,   	 f_0^\hbar        \>
%	 	+
%	 	\hbar
%	 	\int_0^t  \< h\circ \Phi^\hbar (t-s) ,  R_\hbar   [  f^\hbar (s) ]  \>  \d s   \  .
%where $\Phi^\hbar \equiv \Phi_{ \vp^\hbar		}$
%$\widetilde \Phi ^\hbar \equiv \Phi_{ \widetilde \vp ^\hbar		}$.
%Let us analyze the two terms in the above equation, separately.
%	\item \textbf{The first term.}
%Assuming for simplicity that $h(x,p)$ is real-valued, we find
%|		 \langle   h\circ \Phi^\hbar  (t) - h \circ  \widetilde \Phi^\hbar (t) ,   	 f_0^\hbar        \rangle			|
% &  \leq
%  \int_{\R^{d}  \times \R^d }
%  | f_0^\hbar (x,p)|
%  |    h( \Phi^\hbar (t,x,p)		) 	 - h( \widetilde \Phi^\hbar (t,x,p)		) 			|
%  \d x \d p 		\ , \\
%  & \leq
%\|	 h	\|_{Lip}
%  \int_{\R^{d}  \times \R^d }
%  | f_0^\hbar (x,p)|
%  |     \Phi^\hbar (t,x,p)	  - \widetilde \Phi^\hbar (t,x,p)		 		|
%  \d x \d p 		\ , \\
%  &  \leq
%C \exp(C |t|)
%\|	 h	\|_{Lip}
%\|	 f_0^\hbar	\|_{L^1 }
%  \int_0^{|t|}
%  \| \vp^\hbar(s) - \widetilde \vp^\hbar (s)		\|_{L^2} \d s
%	where in the last line we have used continuity of the map $\vp \mapsto \Phi_\vp(t,z)$ (see next Section for details).
%	\item \textbf{The second term.}

\appendix
\section{Well-posedness of the PDEs}
\label{section analysis PDEs}
In this section, we state basic well-posedness results for the Hartree-Hartree, and Vlasov-Hartee equations
that we have introduced in Section \ref{section introduction}.
For notational simplicity, denote by
\begin{equation}
\mathfrak{h}  \equiv L^2(\R^d)
\end{equation}
the one-particle Hilbert space.

\subsection{The Hartree-Hartree Equation}
In what follows, we
consider the
Hartree-Hartree equation
that couples the fermionic reduced density matrix, and the bosonic field.
For notational simplicity,
we assume
microscopic scaling regime--of course,
every result in this section is independent of the
scaling regime under consideration.
That is,  we consider the equation
\begin{align}
\begin{cases}
\label{hartree hartree 1}
& 	i      \partial_t \omega
[   -     \Delta
(V   * \rho_B )  ,
\omega
]			\\
& 	i     \partial_t  \vp
-    \Delta  \vp
(V  *  \rho_F)
\vp   \  , \\
& 	(\omega(0) , \vp(0) = (\omega_0 , \vp_0 )
\end{cases}
\end{align}
for some initial datum
$(ω_0 , _0 )
∈ℒ^1( )
Here we will be employing the notation
for the bosonic and fermionic position densities,
for $(t,x) ∈×^d : $
\begin{equation}
\rho_B(t,x) \equiv
|   \vp (t,x)	|^2
\qquad
\t{and}
\qquad
\rho_F (t,x) \equiv
\omega(t;x,x) \ .
\end{equation}

\vspace{1mm }
Here, we only consider bounded potentials.
The analysis of mean-field equations
with such interactions
is   classical.
Hence, we state
the main results in the next Theorem and omit the proofs.
For instance, we refer the reader to
e.g [12]
whose proof can be adapted to the problem at hand.

\begin{theorem}
\label{theorem wp 1}
Let $(\omega_0 , \vp_0 ) \in
\mathscr{L}^1( \h  ) \times \h
with $\omega_0^* = \omega_0$;
and assume the interaction potential  is bounded $V\in L^\infty(\R^d,\R) $.
the following statementes hold true
\begin{enumerate}[leftmargin=0.7cm]
\item
{\rm{(Global well-posedness)}}
There exists a unique global solution
$(\omega, \vp) \in C ( \R, \mathscr{L}^1( \h ) \times \h ) $
to the Hartree-Hartree equation \eqref{hartree hartree 1} in mild form:
\begin{align}
\label{hartree hartree mild form}
\omega(t) & =
e^{- it \Delta}
\omega_0
e^{+  it \Delta}
\int_0^t
e^{- i (t-s ) \Delta}
\big[
V* \rho_B(s) , \omega(s)
\big]
e^{ i (t-s ) \Delta}
\d s  \\
\vp(t)
&  =
e^{- i t \Delta}
\vp_0
\int_0^t
e^{- i (t-s ) \Delta}
\big(
V * \rho_F(s)
\big) \vp (s)
\d s \ .
\end{align}
Furthermore, there is continuity with respect to the initial data.

\item
{\rm (Unitary evolution)}
Let $(\omega,\vp)$ be the mild  solution of \eqref{hartree hartree mild form},
consider the time-dependent,
mean-field Hamiltonians
on $H^2(\R^d)$
\begin{align}
h_F (t) =  -\Delta + V*\rho_B(t)
\qquad
\t{and}
\qquad
h_B(t)  = - \Delta + V*\rho_F (t) \ .
\end{align}
consider the
two-parameter unitary evolution groups on
\h
that the mean-field Hamiltonians generate
\begin{align}
\begin{cases}
& 	i \partial_t U_F (t,s) = h_F(t) U_F (t,s) \\
& U_F (t,t) = \1
\end{cases}
\quad
\t{and}
\quad
\begin{cases}
& 	i \partial_t U_B (t,s) = h_B (t) U_B (t,s) \\
& U_B (t,t) = \1
\end{cases} \ .
\end{align}
Then, for all $t \in \R$
\begin{equation}
\omega(t) = U^*_F (t,0) \omega_0 U_F (t,0)
\qquad
\t{and}
\qquad
\vp(t)
(t,0) \vp_0 \ .
\end{equation}
In particular,
$\|	 \omega(t)\|_{\tr} = \|	 \omega_0	\|_{\tr}$,
$\omega(t)^* = \omega(t)$, and
$\| \vp(t) 	\|_{L^2 } = \| \vp_0	\|_{L^2}$
for all $t \in \R$.
Additionally, if $\omega_0 \geq 0 $,
then $\omega (t) \geq 0$
and, similarly,
if $\omega_0^2= \omega_0$,
$\omega(t)^2 = \omega(t)$.

\vspace{1mm}

\item
{\rm (Strong solutions)}
Consider the spectral decomposition
$\omega_0 = \sum_{k=0}^\infty
\lambda_k
\ket{\phi_k }
\bra{\phi_k},
and assume that
$ \sum_{k=0}^\infty  | \lambda_k|   \|  \phi_k\|_{H^2}^2    $
$\|	\vp_0\|_{H^2}$
are finite.
the  solution map  is continuously differentiable
$ (\omega , \vp) \in
C^1(\R; \mathscr{L}^1( \h  )\times \h )  $,
\eqref{hartree hartree 1}
holds in the strong sense.

\end{enumerate}

\end{theorem}

\subsection{The Vlasov-Hartee Equation}
In this section, we analyze the coupled Vlasov-Hartee equation that we have introduced in Section \ref{section introduction}.
Namely, the coupled PDEs
\begin{align}
\label{vlasov hartree 2}
\begin{cases}
& 		(\partial_t  + p\cdot  \nabla_x + F_\vp  (t) \cdot    \nabla_p   ) f  = 0 	\\
&		  i \partial_t \vp
- \frac{1 }{2} \Delta   \vp +   V_f(t)  \vp 	\\
&  (f,\vp) ( 0 )  = (f_ 0 , \vp_0)
\in L^1_+(\R^{2d})\times L^2(\R^d)
\end{cases}
\end{align}
$(f_0,_0)∈L_+^1(^2d) ×L^2(^d)$
is some initial  datum.
Here and in the sequel, we shall be using the following notation
for the boson-driven  force field,
and fermion-field potential energy,
for $(t,x) ∈×^d$:
\begin{equation}
F_\vp (t,x)
\equiv
- \int_{\R^d }
\nabla V(x-y ) | \vp(t,y)|^2 \d y
\qquad
\t{and}
\qquad
V_f	(t,x)
\equiv
\int_{	\R^{2d}	}  V(x-y ) f (t,y,p ) \d p  \d y
\end{equation}

\vspace{2mm}
The notion of solution we use is the following.
Let  $  I = (a,b) ⊂$ be an open interval
containing $0$.
We say that a bounded measurable map
$(f , )  :  [a,b] →L^1(^2d)×L^2(^d)$
is a \textit{weak-mild} solution on $I $ of the
Vlasov-Hartee equation \eqref{vlasov hartree 2}
with intial data $(f_0, _0)$,
if the following three conditions are satisfied:
\begin{itemize}[leftmargin=.8cm, label=$\square$]
\item For all $ h  \in  \mathscr S (\R^{2d}) $
the map $t \mapsto  \<  f(t) , h\>$ is differentiable on $I $
\item For all $ h \in \mathscr S (\R^{2d})$
and $ t \in I $
it holds true  that
\begin{align}
& 		 \frac{d}{dt}
\<  f(t), h \>
= \<	 f(t) , (p \cdot  \nabla_x +
F_\vp(t) \cdot   \nabla_p  ) h 	  	 	\> \\
\vp(t)  =   e^{- i t \Delta/2 } \vp_0
- i
\int_0^t e^{ - i (t-s ) \Delta/2 }
V_f (s)
\vp(s) \d s
\end{align}
\item $(f, \vp) (0)  =  (f_0 , \vp_0). $
\end{itemize}
We say that a solution is local-in-time if $I ≠$,
and global-in-time if $ I = $.

\vspace{2mm}

Similarly as for the Hartree-Hartree equation,
since we consider here only  regular potentials
we do not present a proof of the following well-posedness theorem.
Rather, we refer the reader to
\cite[Appendix A]{BenedikterPortaSaffirioSchlein2016}
for a similar result whose proof can be adapted to our problem.

\begin{theorem}
\label{theorem wp 2}
Let $(f_0, \vp_0) \in
L_+^1(\R^d) \times L^2(\R^d)
and assume that
$\nabla V \in \mathrm{Lip}(\R^d; \R)$.
Then, the following statements hold true .
\begin{enumerate}[leftmargin=0.7cm]
\item
{\rm (Global well-posedness)}
There exists
\in L^\infty
(\R ;
L_+^1(\R^d) \times L^2(\R^d)
unique global
weak-mild solution
to the   Vlasov-Hartee equation \eqref{vlasov hartree 2}.
Furthermore, there is continuity with respect to the initial data.

\vspace{1mm}
\item
{\rm (Characteristics/Unitary Representation)}
be the global weak-mild solution
to the Vlasov-Hartee equation.
We denote by
$\Phi_\vp(t) : \R^{2d}\rightarrow \R^{2d}$
volume-preserving diffeomorphism,
corresponding to  	the solution map of the ODE
\begin{align}
\label{ODE}
\begin{cases}
& 			 d/dt  \, x (t) = p (t) 		\\
& d/dt \, p (t)	 = F_\vp(t,  x(t) )	\\
&  (x,v) (0) = (x_0 , p_0 )
\end{cases}  \  ;
\end{align}
and we denote by
$( U_f (t,s)  )_{t,s\in \R}$
the two-parameter
family of unitary transformations defined through
\begin{equation}
\begin{cases}
& i \partial_t U_f (t,s)
\big(
- (1/2)\Delta
V_f (t)
\big)
U_f(t,s) \\
& U_f (t,t) = \1
\end{cases} \ .
\end{equation}
Then, for all $t \in \R$
there holds
\begin{align}
f(t) = f_0 \big(
\Phi_\vp^{-1} (t)
\big) \\
\vp(t) = U_f(t,0) \vp_0  \ .
\end{align}
In particular, $\|	f(t)	\|_{L^1} = \|f_0	\|_{L^1}$,
$f(t) \geq0$
$\|	 \vp(t)\|_{L^2} = \| \vp_0	\|_{L^2}, $
for all $ t \in \R $.

\vspace{1mm}

\item
{\rm (Strong solutions)}
Assume additionally that
$\|	   (1 + |x| + |p|)   \< \nabla_{x,p }\> 	 f_0\|_{L^1}$
$\|	 \vp_0	\|_{H^2}$
are finite.
the solution map
is differentiable
$(f, \vp ) \in
C^1(\R; 		L^1_+(\R^{2d})\times L^2(\R^d 		)
and the Vlasov-Hartee equation \eqref{vlasov hartree 2}
holds in the strong sense.
\end{enumerate}

\end{theorem}

\begin{remark}
\label{remark grow}
Let $(f,\vp)$ be the weak-mild solution
of the Vlasov-Hartree equation extracted from Theorem \ref{theorem wp 2},
and let $\Phi_\vp(t) $
be the solution map of the associated ODE.
Then, a Gr\"onwall argument shows
that there exists a map
$\kappa : \R \mapsto (0,\infty)$
such that the following bound is satisfied
\begin{equation}
| \Phi_\vp(t,z) |
\leq \kappa (t) |z| \ , \qquad \forall t\in \R , \ \forall z \in \R^{2d} \ .
\end{equation}
\end{remark}

%			Before we jump into the proof
%			of existence and uniqueness of weak/mild solutions,
%			we need a preliminary result, on representation of these solutions.
%			First, we proceed to define
%			the study the associated
%			ODE problem.
%			\begin{lemma}
%				Let us assume that the interaction potential is such that
%				$\nabla V \in \mathrm{Lip} (\R^d; \R^d )$.
%				\begin{enumerate}[leftmargin=* ]
%					\item
%					Let
%					$
%					\vp \in L^\infty (\R  ;  L^2)
%					$.
%					Then, for all $z_0 = (x_0 , v_0) \in \R^{2d}$
%					there exists a unique solution
%					$ z  = (x,v) \in C^1(	\R ; \R^{2d}  	)$
%					to the ODE
%					\begin{align}
%						\label{ODE}
%						\begin{cases}
%							& 			 d/dt  \, x (t) = v(t) 		\\
%							& d/dt \, v(t)	 = F [ \vp(t) ] ( x(t) )	\\
%							&  (x,v) (0) = (x_0 , v_0 )
%						\end{cases}  \ .
%					\end{align}
%					We denote  by $\Phi_\vp ( t, z) $
%					the solution to the   ODE \eqref{ODE}, with initial data $z \in \R^{2d}$, at time $t\in\R$.
%					For all $t \in \R$,
%					the map $\Phi_\vp(t) : \R^{2d} \rightarrow \R^{2d}$
%					is a volume-preserving diffeomorphism.
%					\item
%					There  is continiuity with respect to the initial data.
%					Namely, for all $   \vp \in C_b(\R  ;  L^2)$   there exists $C>0$ such that
%					for all $z_1,z_2 \in \R^{2d}$ and $t \in \R$, there holds
%					\begin{equation}
%						|	 \Phi_\vp ( t , z_1 ) 		-		 \Phi_\vp ( t,z_2) 	|
%						\leq  C \exp(C|t|) |   z_1 - z_2  |   \ .
%					\end{equation}
%					\item
%					There is continuity with respect to the boson field.
%					Namely, there exists a constant $C>0$
%					such that
%					for all $z \in  \R^{2d}$, $t \in \R$
%					and $\vp_1, \vp_2 \in   C_b(\R  ;  L^2)$,
%					there holds
%					\begin{equation}
%						|	 \Phi_{\vp_1 } ( t , z ) 		-		 \Phi_{\vp_2} ( t,z ) 	|
%						\leq  C    \exp(C |t|)
%						\int_0^{|t|} \|	   \vp_1 (s)		- \vp_2(s)	\|_{L^2}  \d s \  .
%					\end{equation}
%				\end{enumerate}
%			\end{lemma}
%			We can now relate the solutions of the ODE \eqref{ODE},
%			with the weak/mild solutions of the coupled Hartree system.
%			\begin{lemma}
%				The map
%				$(f,\vp) \in C( [a,b] ;  L_+^1 \times L^2	)$
%				is a weak/mild solution to the coupled Hartree system, with initial data $(f_0 , \vp_0)$,
%				if and only if
%				for all   $ t  \in  [a,b ] $
%				the following identities  holds true
%				\begin{align}
%					&
%					f ( t; x, v	 	)
%					=  f_0 (  \Phi_{	 \vp }^{-1 }  (t ; x,v) 	)
%					\ , \qquad  \t{ a.e } (x,v )\in \R^{2d }
%					\\
%					&
%					\vp(t)  =   e^{- i t \Delta/2 } \vp_0
%					- i
%					\int_0^t e^{ - i (t-s ) \Delta/2 }
%					\nu [f(s) ]
%					\vp(s) \d s
%				\end{align}
%			\end{lemma}
%			\begin{proposition}[Local Well-posedness]
%				\label{prop local wp}
%				For all initial data $(f_0 , \vp_0) \in L^1_+  \times L^2 $,
%				there exists $T =  T (  V, \|	 f_0	\|_{L^1} , \|	 \vp_0\|_{L^2}  )>0$
%				such that
%				there is  a unique
%				weak/mild solution
%				$  (f, \vp) \in  L^\infty   (  [0,T] ;  L_+^1  \times L^2  )$
%				to the coupled Hartee system, with initial data $(f_0 ,\vp_0)$.
%			\end{proposition}
%			\begin{proof}
%				We use a fixed point argument, and only   give a sketch of the proof.
%				We introduce the map
%				\begin{equation}
%					\M = (\M_1 , \M_2) :
%					X_T \equiv 	L^\infty ( [0,T] ;  L_+^1 \times L^2	)\rightarrow
%					X_T = 	L^\infty  ( [0,T] ;  L_+^1 \times L^2	)
%				\end{equation}
%				as follows.
%				For $ (f,\vp)		 \in 	X_T $ we let
%				\begin{align}
%					&
%					\M_1[f,\vp] (t ; x, v)
%					\equiv
%					f_0 (  \Phi_{	 \vp }^{-1 }  (t ; x,v) 	)
%					\ , \qquad  \t{ a.e } (x,v )\in \R^{2d }
%					\\
%					&
%					\M_2 [f,\vp] (t)
%					\equiv
%					e^{- i t \Delta/2 } \vp_0
%					- i
%					\int_0^t e^{ - i (t-s ) \Delta/2 }
%					\nu [f(s) ]
%					\vp(s) \d s   \ .
%				\end{align}
%				$\M_1$ is bounded because $\Phi_\vp (t)$
%				is a volume-preserving diffeomorphism.
%				$\M_2$ is bounded because
%				$ exp(- i t \Delta)$ is a unitary transformation in $L^2$,
%				and $ \| \nu [f(s)]  \|_{L^\infty} \leq   \|V\|_{L^\infty } \|	 f (s)\|_{L^1 } .$
%				Hence, $\M (X_T) \subset X_T$.
%				Furhter, let $R_1 \equiv \|	f_0	\|_{L^1}$
%				and $R_2 \equiv \|	 \vp_0	\|_{L^2 }$
%				and set
%				\begin{equation}
%					X_T(R_1, R_2)
%					\equiv
%					\{	 (f,\vp) \in X_T :   \|		f\|_{L^\infty L^1}	 \leq R_1 , \
%					\|	\vp	\|_{L^\infty L^2} \leq 2 R_2						 		\} \ .
%				\end{equation}
%				We may choose $T = T(V, R_1, R_2)$
%				such that $\M_T$ maps $X_T(R_1, R_2) $ into itself.
%				It suffices to show that it is a contraction  with
%				respect to an appropiate metric.
%				We define the metric
%				\begin{equation}
%					d \big[	 (f_1 , \vp_1)	 , (f_2 , \vp_2)	\big]
%					\equiv
%					\sup_{t\in[0,T]}
%					\sup_{   |h|_{\mathrm{Lip}} \leq 1 }
%					|    \<  f_1  (t) - f_2 (t) , h\> 	|
%					+
%					\sup_{t\in[0,T]}
%					\|	 \vp_1 (t)  - \vp_2 (t) 	 \|_{L^2 }
%				\end{equation}
%				Let $(f , \vp) , ( \widetilde f , \widetilde \vp )
%				\in X_T (R_1, R_2)$.
%				We find that for $ h \in \mathrm{Lip}$
%				\begin{align}
%					\nonumber
%					\sup_{t\in [0,T]}
%					|  \<  \M_1 [f,\vp ] (t)  -   \M_1 [\widetilde f , \widetilde \vp  ]  (t) , h 	\> |
%					& 	 \leq
%					\|	 h	\|_{\mathrm{Lip}	}
%					\|	 f_0	\|_{L^1 }
%					\sup_{t\in [0 ,T ]}
%					\sup_{z \in \R^{2d} }
%					|		\Phi_\vp (t,z )  -  \Phi_{\widetilde \vp} (t,z ) 	| \\
%					& \leq
%					\nonumber
%					\|	 h	\|_{\mathrm{Lip}	}
%					\|	 f_0	\|_{L^1 }
%					C e^{CT } T
%					\|	 \vp (t)   - \widetilde \vp(t) 	\|_{L^2}^2  \\
%					\nonumber
%					& \leq
%					\|	 h	\|_{\mathrm{Lip}	}
%					\|	 f_0	\|_{L^1 }
%					C e^{CT } T
%					( \| \vp	\|_{L^\infty L_2}  +
%					\|  \widetilde \vp	\|_{L^\infty L_2}   )
%					\| \vp - \widetilde \vp 	\|_{L^\infty L_2}    \\
%					& \leq
%					\|	 h	\|_{\mathrm{Lip}	}
%					4 C R_1 R_2
%					e^{CT } T
%					\ d[ (f, \vp)  , (\widetilde f,  \widetilde \vp)	] \ .
%				\end{align}
%				On the other hand,
%				since $V(x-y)$ is Lipschitz
%				\begin{align}
%					\nonumber
%					\|	\M_2[f ,\vp]   (t)		  -  \M_2  [\widetilde f , \widetilde \vp] (t)				\|_{L^2}
%					& 	\leq
%					\int_0^t
%					\|	 \nu[f(s)] \vp(s)     -  \nu[\widetilde f(s)]  \widetilde \vp (s) 	\|_{L^2 }
%					\d s  \\
%					\nonumber
%					& 	\leq
%					T
%					\sup_{s \in [0,T]}
%					\|	 (  \nu[f(s)]       -  \nu[\widetilde f(s)]   )    \vp (s) 	\|_{L^2 }  \\
%					\nonumber
%					& +
%					T  \sup_{s \in [0,T]}
%					\|	  \nu [ \widetilde f (s)	]	 (\vp (s)  -  \widetilde \vp (s) ) 	\|_{L^2 }  \\
%					\nonumber
%					& \leq
%					T   | V|_{\mathrm{Lip}}
%					\sup_{s \in [0,T]}
%					\sup_{   |h|_{\mathrm{Lip}} \leq 1 }
%					|    \<  f_1  (t) - f_2 (t) , h\> 	|	 \|	\vp	\|_{L^\infty L^2}	 	\\
%					\nonumber
%					&  +
%					T     |V|_{L^\infty}
%					\|	 \widetilde f 	\|_{L^\infty L^1}
%					\|	 \vp - \widetilde \vp	\|_{L^\infty L^2}	\\
%					& \leq
%					T C_V  R_1 R_2 \
%					\ d[ (f, \vp)  , (\widetilde f,  \widetilde \vp)	] \ .
%				\end{align}
%				Thus, it suffices to combine the two previous inequalities and take $T$ sufficiently  small,
%				depending only on $R_1$, $R_2$, and constants depending on $V$.
%			\end{proof}
%			\begin{proof}[Proof of Theorem \ref{theorem wp}]
%				Let $ (f_0, \vp_0)\in L^1_+(\R^6) \times L^2(\R^3 )$
%				and assume in addition that $\vp_0 \in H^2$.
%				Let $(f,\vp)$ be the local-in-time solution of the Vlasov-Hartee equation, as
%				given by Proposition \ref{prop local wp}, with maximum existence time $T>0$.
%				Since $\vp_0\in H^2$,   standard techniques from semi-linear Schr\"odinger equations    imply that
%				$\vp(t) $  is strongly differentiable with
%				$\vp\in C^1\big((0,T) , L^2\big)$.
%				In particular, the equation for the boson field now holds in the strong sense
%				\begin{equation}
%					i \partial_t \vp  = ( - (1/2) \Delta + \nu[f] )\vp    \ .
%				\end{equation}
%				Thus,   it holds that   on $   (0,T)$
%				\begin{equation}
%					\partial_t \|	 \vp   	\|_{L^2}^2 =
%					2	\mathrm{Im } \<  \vp    ( - (1/2) \Delta + \nu[f] )  \vp   \> = 0 \ .
%				\end{equation}
%				Hence, $ \|	 \vp (t)	\|_{L^2 } = \|	 \vp_0	\|_{L^2}$ for all $t \in [0,T]$.
%				Thus, using the density of $H^2 \subset L^2$
%				and the $L^2$ continuity with respect to the initial data,
%				one may conclude that the same conservation of $L^2$ norm holds,
%				for $\vp_0 \in L^2$.
%				On the other hand, for the
%				fermion density, since $\Phi_\vp (t)$ is always volume-preserving,
%				there holds $ \|	 f(t)	\|_{L^1} = \|		 f_0\|_{L^1}.	$
%				Since $T = T (V ,\|	f_0	\|_{L^1} , \| \vp_0	\|_{L^2}	)$
%				depends only on conserved quantities,
%				one may then patch local-in-time solutions
%				to obtain global well-posedness.
%			\end{proof}

\section{Calculation of the Infinitesimal Generator}
\label{appendix inf generator}
In this section, we give
more details of the calculation
of the infinitesimal generator
of the
fluctuation   dynamics $(t,s)$,
introduced in \eqref{unitary U}.
This is the time-dependent self-adjoint operator
$Ł(t)$ on $$ defined through the equations
\begin{equation}
i \hbar \partial_t \U(t,s)
\mathcal L (t) \U(t,s) \ ,
\qquad
\U(t,t)=\1
\end{equation}
where $ t,s∈$.
In what follows,
we let
$(ω,)$  be a strong  solution
of the Hartree-Hartree equation \eqref{hartree hartree},
so that the unitary maps
$_t$ and $_t$
defined in \eqref{definition W and R}
are differentiable with respect to $t∈$.
Note that the final result of this section is contained in Proposition \ref{prop appendix}
and requires no $H^2$ regularity
of the solution $(ω,)$.
an approximation argument
shows that the result also holds for
mild solutions--we leave the details to the reader.

\vspace{2mm}

Our purpose here is to given an explicit representation of $Ł(t)$ in terms
of creation- and annihilation-
As a first step in our calculation, we see that the unitarity of the maps easily imply that
$Ł(t)$ is the contribution of three terms. Namely, for all $ t ∈$ there holds
\begin{equation}
\label{generator 1}
\L(t)
\i \hbar \partial_t \calR_t^*  \calR_t \otimes \1
\1 \otimes \i \hbar \partial_t \W_t^* \W_t
( \calR_t \otimes \W_t 		)^*
\H
(\calR_t \otimes \W_t ) \ .
\end{equation}
We now proceed to  calculate each term separately.

\subsection{Calculation of $\partial_t \calR^*_t \calR_t $}
Let us first recall some facts and notations
that we have   introduced
in Section \ref{section second quantization}.
Namely, denoting by $ω_t ≡ω(t) = ∑_i=1^M |ϕ_i(t)⟩    ⟨ϕ_i(t)|$
the fermion component of the solution of the Hartree-Hartree equation \eqref{hartree hartree},
we let
$u_t ≡u(t)$ and $  v_t  ≡v(t)$
be the bounded operators on $L^2(^d)$
defined as
$u_t =    1 - ω_t
|  ϕ_i(t) ⟩
⟨ ϕ_i(t)  | .
The kernels of these operators define the  distributions
$	u_t,x (y)  ≡u_t (y,x)$
$				v_t,x (y) ≡v_t (y,x) $,
for all $t∈$ and $x,y∈^d$.
let us recall that we have introduced
in \eqref{hamiltonians} the  operator $h_F(t)$
as   the one-particle   Hamiltonian,
driving the fermion dynamics.
In particular,
one may verify that for all $t∈$ there holds
\begin{align}
\label{u and v equations}
\i \hbar \partial_t u_t   = [ h_F (t) , u_t  ]
\qquad
\t{and}
\qquad
\i \hbar \partial_t \overline{ v_t }
& =
h_F(t) \overline{v_t}
\overline{v_t} \overline{ h_F(t) }  \  .
\end{align}

\vspace{2mm}

In order to calculate the first term in the expansion \eqref{generator 1},
we start with the following auxiliary Lemma.
The  proof is a simplification of the argument contained in
\cite[Proposition 3.1]{BenedikterJaksicPortaSaffirioSchlein2016}
for   mixed states.
Note that a similar calculation has been carried out in
[6] for pure states,
which unfortunately is not
precise enough for our interests.

\begin{lemma}
\label{lemma generator R}
Assume that for all $j=1, \ldots, M $,
the orbitals are differentiable
in $L^2$,
$t\mapsto \phi_j(t)$.
Then, we find that for all $t\in \R$
there holds
\begin{equation}
\i \hbar \partial_t
\calR_t^* \, \calR_t
\mathcal{S}_F(t) \1
\int_{\R^{2d}} \calC  (t; x,y) a_x^* a_y \d x \d y
\frac{1}{2}
\Bigg(
\int_{ \R^{2d}}
\mathcal{D}(t ; x,y) a_x^* a_y^* \, dx dy
+ h.c
\Bigg)  \ .
\end{equation}
Here, $\calC (t)$ and $\calD (t) $ are  the operators on $L^2(\R^d)$
\begin{equation}
\calC (t)  =
(i \hbar \partial_t u_t) u_t
(i\hbar \partial_t \overline{v_t}) v_t
\qquad \t{and} \qquad
\calD (t)   =
(i \hbar \partial_t u_t) \overline{v_t}
(i\hbar \partial_t \overline{v_t}) \overline{u_t}       \ ,
\end{equation}
and $\calS_F(t)\in \R$ is the scalar term.
\begin{align}
\textstyle
\mathcal{S}_F(t)   =
- \sum_{j=1}^M
\<
\phi_j(t)
i \hbar \partial_t
\phi_j(t) \>_{L^2} \ .
\end{align}

\end{lemma}

\begin{proof}[Sketch of proof]
Let us drop the time labels
in order to simplify our notation:
$\calR_t = \calR$,
$u_t = u$
$v_t = v$.
We start with the following observation:
letting $f\in L^2(\R^d)$ be a test function,
we obtain the conjugation relations
\begin{equation}
\calR^* a^* (f) \calR
a(\overline{vf}) \ .
\end{equation}
and thus, taking the time derivative on both sides
and using $(d/dt \calR )^* \calR = - \calR^* (d/dt \calR)$
we obtain
\begin{equation}
[ (d/dt \, \calR^*) \calR , \calR^* a^*(f) \calR    ] =
a^*\Big(  \frac{du}{dt} f \big)
\frac{d  \overline v}{dt}
\overline{f}
\Big) \ .
\end{equation}
We plug in again  the conjugation relations back in the previous equation and multiply both sides with $i$
to to obtain
\begin{equation}
\label{lemma b2 eq 1}
( i \, d/dt \, \calR^*) \calR
a^*\Big(  i \frac{du}{dt} f \big)
i \frac{d  \overline v}{dt}
\overline{f}
\Big) \ .
\end{equation}
The previous equation
the self-adjoint
$( i \, d/dt \, \calR^*) \calR$
up to a scalar term,
and implies that
is quadratic in creation-
and annihilation operators.
there exists operators
$C$ and $D$
$L^2(\R^d )$
with operator kernels
and $D(x,y)$,
and a scalar $S_F \in \R $
such that
\begin{equation}
( i \, d/dt \, \calR^*) \calR
S_F  \1
\int_{\R^{2d}} C(x,y) a_x^* a_y \d x \d y
\frac{1}{2}
\Bigg(
\int_{\R^{2d}}
D(x,y) a_x^* a_y^* \d x \d y
+ h.c
\Bigg)  \ .
\end{equation}
Here,   we may assume without loss of generality that $\overline{C(y,x)} = C(x,y)$
$D(y,x) =  - D(x,y)$.
Let us first compute that operator contributions.
A lengthy but straightforward calculation using the CAR
allows us to compute the commutators
\begin{align}
\label{lemma b2 eq 2}
[ ( i \, d/dt \, \calR^*) \calR , a^*(uf)
a^*(C uf )
- a(C \overline{v} \overline{f})
+ a^*(D v f )
a (D \overline u \overline f) \  .
\end{align}
Thus, we compare the right hand sides of  \eqref{lemma b2 eq 1} and \eqref{lemma b2 eq 2}
to obtain, as operators on $L^2(\R^d)$, the following equations
\begin{align}
\begin{cases}
C u + D v &  = i  du/dt  \\
C \overline v + D \overline u & = i d \overline v / dt  \  \\
\end{cases}
\iff
(C \ D )
\begin{pmatrix}
u & \overline v \\
v & \overline u \\
\end{pmatrix}
= i d/dt ( u \ \overline v  ) \ .
\end{align}
In particular, the matrix  containing $u$ and $v$
is a unitary map on $L^2(\R^d) \oplus L^2(\R^d)$,
which we denote by $\nu$.
In fact, using the relations
$u^* = u $,
$v^* = \overline v$
it is easy to verify that  $\nu^* = \nu$.
Thus we invert the last equation and solve for the operators  $C$ and $D$
to find
\begin{equation}
C =  i  du/dt   \, u + i d \overline v /dt  \, v
\qquad
\t{and}
\qquad
D = i du/dt   \, \overline v + i d \overline v /dt \, \overline  u \ .
\end{equation}
Finally,  we multiply with $\hbar$
and identify
$\calC = \hbar C$
$\calD = \hbar D$.
As for the scalar contribution, we consider its vacuum expectation value to obtain
\begin{align}
=  \<  \Omega_F , i d\calR^* /dt \, \calR  \Omega_F\>_{\F_F}
= - \<     \calR \Omega_F , i  d\calR / dt \, \Omega_F    \>_{\F_F} \ .
\end{align}
Let us now calculate the right hand side.
To this end, we write in terms of the orbitals
$\omega(t) = \sum_{j=1}^M \ket{\phi_j(t)} \bra{\phi_j(t)}$
the vector
$\calR \Omega_F  = a^*( \phi_1) \cdots a^*(\phi_M)\Omega_F$.
A straightforward calculation
using the CAR and the orthogonality relations
$\<\phi_i, \phi_j\>_{L^2}= \delta_{i,j}$ now implies
\begin{equation}
= - \sum_{j=1}^M \< \phi_j (t),  i d \phi_j(t)/dt \>_{L^2} \ .
\end{equation}
This finishes the proof once we multiply with $\hbar$
and identify $\calS = \hbar S $.
\end{proof}

Let us now give a more explicit representation
of the term that we just calculated.
Namely, in view of \eqref{u and v equations}
it is easy to verify that
\begin{equation}
\mathcal C (t)
u_t h_F(t) u_t
\overline{v_t}   \overline{h_F(t)}   v_t
\qquad
\t{and}
\qquad
\mathcal D (t)
- u_t h_F(t) \overline{v_t}
\overline{v_t} \overline{h_F(t)} \overline{u_t}
\end{equation}
In particular, the second term may be simplified. Namely,  in view of the relations
$h_F (t)  = h_F^*(t) $,
$ u_t = u^*_t$ and $ v̅ _t  = v_t^*  $, one may check that
$ (vhu)^* = u h v̅	$.
Consequently, using the  anticommutation relation ${a_x^* ,  a_y^*} = 0 $,
one  finds that
\begin{equation}
\int_{ \R^d \times \R^d }  (  \overline{v h u}  )  (x,y) a_x^* a_y^* \d x \d y
\int_{ \R^d \times \R^d }
( u h \overline v)(x,y)
a_x^* a_y^* \d x \d y \ .
\end{equation}
We can then put the calculation from the above lemma in the following form,
which is the final result of our calculation
(we omit the time labels for convenience)
\begin{equation}
( 	i \hbar \partial_t \calR^*)   \calR
- \tr \big(
h_F \omega
\big)
\d \Gamma_F \Big[
u  h_F u
\overline{v}   \overline{h_F}   v
\Big]
\Bigg(
\int_{ \R^{2d} }
[ u h_F  \bar{v} ](x,y) a_x^* a_y^*
\d x \d y
+ h.c
\Bigg)
\end{equation}

\subsection{Calculation of $\partial_t \W^*_t \W_t $}
Let us now calculate the second contribution of \eqref{generator 1}.
The time derivatives of Weyl operators that are parametrized
by a   field $t↦α(t) ∈L^2(^d)$
can be regarded  as classical result,
and we record it in the following lemma.
For reference, see
\cite[Lemma 3.1]{GinibreVelo1974}.

\begin{lemma}
Assume that
$t \mapsto \vp_t \in L^2(\R^d)$
is   differentiable.
Then, for all $t\in \R$ there holds
\begin{equation}
i \hbar \partial_t
\W_t^* \W_t
N  \mathrm{Im}
\< \vp_t ,  \hbar  \partial_t \vp_t \>
\sqrt N
\Big(
b^* ({ i \hbar \partial_t \vp_t}   )
b  ( i \hbar \partial_t \vp_t)
\Big)  \ .
\end{equation}
\end{lemma}

We shall use the fact that $(t)$
solves the Hartree-Hartree equation \eqref{hartree hartree}.
Namely, we consider on $H^2(^d)$
the time-dependent bosonic   Hamiltonian
defined in \eqref{hamiltonians} and we conclude that for
all $t∈$
there holds
\begin{equation}
i \hbar \partial_t
\W_t^* \W_t
N \mathrm{Re}
\< \vp_t , h_B(t)  \vp_t \>
\sqrt N
\Big(
b^* ( h_B(t)  \vp_t   )
b  ( h_B(t)  \vp_t)
\Big)  \ .
\end{equation}

\subsection{Calculation of $\calR^*_t \W^*_t \H \W_t \calR_t $}

Our next task it to compute
( _t ⊗_t 		)^*
(̋_t ⊗_t )  .
To this end, we shall use extensively
the conjugation relations
for particle-hole transformations   (see Lemma \ref{lemma fermion properties})
\begin{align}
\calR_t^* a^*_x \calR_t
a^* ( u_{t,x})
a( \overline{  v_{t,x} } )
\qquad
\t{and}
\qquad
\calR_t^* a_x \calR_t
a ( u_{t,x})
a^*( \overline{  v_{t,x} } ) \ ,
\end{align}
for all $(t,x) ∈×^d$, and similarly for Weyl operators (see Lemma \ref{lemma properties bosons})
\begin{align}
\W_t^* b_x  \W_t
b_x +  	\sqrt  N \vp_t(x)
\qquad
\t{and}
\qquad
\W_t^* b^*_x  \W_t
b^*_x +  \sqrt N 	\overline{ 	\vp_t(x)  } \  .
\end{align}
where $(t,x) ∈×^d$.
A lengthy but otherwise straightforward calculation
using these conjugation relations,
together with $u^*=u$
yields the following result.
\begin{lemma}
\label{lemma generator conjugation}
The following holds for all $t \in \R$
\begin{align}
( \calR_t \otimes &  \W_t 		)^*					\nonumber
\H
(\calR_t \otimes \W_t )  \\
& 	=
\bigg(
N \|	  \nabla \vp_t	\|_{L^2}^2 		+ \t{Tr }( -\Delta\omega_t )
\nonumber
\lambda NM
\int_{ \R^{2d} }
\rho_F(t,x )
\rho_B(t,y)
\, dx dy
\bigg)  \, \1 \otimes \1
\\
&  +
\Bigg(						\nonumber
\d \Gamma_F
\big[
u_t h_F(t) u_t
\overline{v_t} h_F(t) v_t
\big]
\int_{ \R^{2d} }
( u_t h_F(t) \bar v _t)(x,y) a_x^* a_y^* dx \, dy
+ \t{h.c}
\Bigg)
\otimes \1 			 \\
\1 \otimes
\bigg(
\d \Gamma_B
\big[ h_B(t) \big]
\sqrt N \big(	b^* (h_B(t) \vp_t )		+ 	b(h_B(t) \vp_t )	\big)
\bigg)
\nonumber 			\\
&  +  \lambda \sqrt  N \L_{2,1} (t) + \lambda \L_{2,2} (t) \ .
\end{align}
Here, we denote
\begin{align}
\nonumber
\L_{2,1} (t)
& =
\int_{ \R^d \times \R^d  } V(x - y ) a^*(u_{t,x}) a( u_{t,x}) \otimes \big(		 \vp_t (y) b_y^* + h.c	\big)
\ \d x \d y 	\\
\nonumber
& -
\int_{ \R^d \times \R^d  } V(x - y ) a^*( \overline{v _{t,x}}) a(\overline{v_{t,x}}) \otimes
\big(		 \vp_t(y) b_y^* + h.c	\big) \ \d x \d y 	\\
\nonumber
& +
\int_{ \R^d \times \R^d  } V(x - y ) a^*(  {u_{t,x}}) a^*(\overline{v_{t,x}})
\otimes
\big(		 \vp_t(y) b_y^* + h.c	\big)  \ \d x \d y 	\\
& +
\int_{ \R^d \times \R^d  } V(x - y ) a(  \overline{v _{t,x}}    ) a( u_{t,x} ) \otimes
\big(		 \vp_t(y) b_y^* + h.c	\big)
\ \d x \d y
\end{align}
\begin{align}
\nonumber
\L_{2,2} (t)
& =
\int_{ \R^d \times \R^d  } V(x - y )
a^*(  u_{t,x}) a(    u_{t,x} )
\otimes b_y^*b_y
\ \d x \d y 	\\
\nonumber
& -
\int_{ \R^d \times \R^d  } V(x - y )
a^*( \overline{v _{t,x}})
\otimes b_y^*  b_y
\ \d x \d y 	\\
\nonumber
& +
\int_{ \R^d \times \R^d  } V(x - y )
a^*(  {u_{t,x}})
a^*(\overline{  v_{t,x} })
\otimes b_y^*b_y
\ \d x \d y 	\\
& +
\int_{ \R^d \times \R^d  } V(x - y ) a(  \overline{v _{t,x}}    ) a( u_{t,x}  ) \otimes b_y^*b_y
\ \d x \d y 	  \ .
\end{align}
\end{lemma}

\subsection{Putting everything together}

We put together the last three lemmas to find the following explicit representaiton of the generator $Ł(t)$.

\begin{proposition}
\label{prop appendix}
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
\begin{equation}
\L(t)
S(t) \1 \otimes \1
\d \Gamma [h_F(t) ] \otimes \1
\1 \otimes \d \Gamma [h_B(t)]
\lambda \sqrt N \L_{2,1}(t)
\lambda \L_{2,2}(t)  \ .
\end{equation}
the scalar term is
$S(t) = -  \lambda NM \int_{\R^{2d}}V(x-y) \rho_F(t,x) \rho_B(t,y) dxdy$,
$h_F(t)$ and $h_B(t)$ are the 1-particle Hamiltonians defined in
\eqref{hamiltonians}, and
the operators $\L_{2,1}(t) $ and $\L_{2,2}(t)$ are defined in Lemma \ref{lemma generator conjugation}.
\end{proposition}

%\section{Comments on the Initial Data}
%Let us consider here an external potential
%$U: \R^d \rightarrow \R$
%that we regard as a trap.
%Namely, we assume that its negative and positive parts satisfy
%U_- \in L^{1 + d/2}	(\R^d )
%+ L_\ve (\R^d)
%U_+ \in L_{\t{loc}}^1(\R^d) ,
%\ \lim_{|x| \rightarrow \infty}
%U_+ (x) = + \infty \ .
%We consider
%on the $N$-particle space
%$\mathfrak{ h }_N \equiv  L^2_a (\R^{dN})$
%the Hamiltonian
%	H_N
%	\equiv
%	\sum_{ i =1 }^N
%	(- \Delta_{x_i})
%	+
%	\sum_{ i = 1 }^N U (x_i )
%which describes fermions with positions $x_i\in \R^d$, moving under the presence of the trap $U (x_i)$. We shall denote its ground state energy by
%	E( N ) \equiv
%	\inf  \sigma (H_N ) \ , \qquad N \geq 1 \  ,
%and set the semi-classical scaling $\hbar =1/  N^{1/d}. $
%\textbf{Thomas-Fermi theory}.
%Heuristically, in the large $N $ limit one is able to introduce a mean-field description
%by means of the Thomas-Fermi energy.
%More precisely, the leading order term of $E(N)$ is expected to be
%described in terms
%of the functional
%\label{thomas fermi}
%	\mathcal E
%	(\rho)
%	\equiv
%	\frac{d}{d+2}
%	C_{TF}
%	\int_{	\R^d	} \rho(x)^{1 + 2/d} \d x
%	+
%	\int_{	\R^d	} U(x) \rho(x) \d x \  ,
%where $\rho(x)$ denotes the  number  density of the particles at the point $x\in \R^d$, and $C_{TF} = 4\pi^2 (d/ |\mathbb S _{d-1}|)^{2/d}$ stands for the Thomas-Fermi constant.
%The variational problem now reads
%	e_{TF}
%	\equiv
%	\inf
%	 \Big\{
%	 \mathcal E (\rho )
%	 :
%	 0 \leq \rho  \in L^1\cap L^{1+2/d}(\R^d),   \|	 \rho	\|_{L^1}=1
%	 \Big\}  \ .
%Under these conditions, it is known that \cite[Theorem 1.1]{Fournais2018} that
%the Thomas-Fermi energy functional actually describes the many-body quantum problem. Namely, it holds true that
%	\lim_{N\rightarrow \infty } \frac{E(N)}{N}
%	 =
%	 e_{TF} \ .
%\textbf{Convergence of states.}
%What can be said about states?
%Let us first introduce some notation.
%given  a normalized state $\Psi \in \mathfrak{ h }_N$,
%we  define its one-particle reduced density matrix
%as the trace-class operator with kernel
%	\gamma^{(1)}_{\Psi} ( x, x')
%	 \equiv
%	 N
%	  \int_{\R^{d(N-1)}}
%	 \Psi (x , x_1 , \ldots, x_{N_1} ) \overline{\Psi(x' , x_1 , \ldots, x_{N-1})}
%	 \d x_1 \cdots \d x_{N-1} \   \ .
%In the semi-classical limit, one studies the $\hbar\downarrow 0 $
%limit of the Wigner transform, defined as
%the function on phase space
%	 {f}_{\Psi}
%	(x,p)
%	\equiv
%	\int_{	\R^d	}
%	\gamma^{(1)}_\Psi ( x + y/2 , x - y/2 )  e^{- i  \frac{y\cdot p }{\hbar }} \d y  \ .
%Note that we have chosen the normalizations:
%$\t{Tr} \gamma_\Psi = N $, and $\int_{	\R^{2d}	} \mathscr{W}_\Psi(x,p) \d x \d p =1 $.
%Now, let us consider a sequence of approximate minimizers
%of the ground state problem.
%That is, let $ \{  \Psi \}_{N=1} $ be a sequence of normalized states such that
%	\< \Psi_N , H_N \Psi_N \>_{\mathfrak{ h }_N} = E(N) \big(1 + o(1) \big) ,
%	 \qquad
%	 \t{as}
%	 N \rightarrow \infty \ .
%Then--up to extraction of  a subsequence--it is known
%\cite[Theorem 1.2]{Fournais2018} that the following limits hold true
%\label{convergence minimizers}
%	\forall  \vp \in   C_c^\infty (\R^{2d})
%	\lim_{N\rightarrow \infty }
%	 \< \vp , f_{\Psi_N} \>
%	 =
%	 	 \< \vp ,  f_\rho  \>
%	 	 \quad
%	 	 \t{where}
%	 	 \quad
%	 	 f_\rho(x,p) \equiv \mathds{1}(|p|^2 \leq C_{TF} \rho(x)^{2/d})
% where $\rho$ is the \textit{unique} minimzer of the (striclty convex) functional defined in Eq. \eqref{thomas fermi}.
%In the original result, one actually has convergence
%  in the space  of test functions
%satisfying $D^\alpha \vp \in L^\infty (\R^{dk})$ for any multi-index with
%$\max \alpha_i  \leq 1. $
%However, for our purposes, convergence in \eqref{convergence minimizers} is everything we need.
%\textbf{Upgrading convergence}.
%Let us now present an abstract argument
%that shows that one can   upgrade  the conclusion \eqref{convergence minimizers}
%to strong convergence  in  certain Banach spaces of functions.
%In what follows, using the same notation as above, we shall denote for simplicity
%	f_N \equiv f_{\Psi_N} , \  N \geq 1  , \qquad \t{and} \qquad  f \equiv f_\rho \ .
%We will also introduce the associated \textit{Husimi} measure
%\label{husimi def}
%	m_N \equiv f_N * \mathscr{G}_{\hbar } , \qquad   \hbar = \frac{1}{N^{1/d}}, N \geq 1
%where $\mathscr{G}_\hbar (x,p) \equiv  \hbar^{-2d} \mathscr{G}_1 ( \hbar^{-1} (x,p))$
%is the approximation of the identity
%with Gaussian profile
%$\mathscr{G}_1 (x,p) \equiv  \pi^{-d/4} \exp(  x^2+p^2 )$.
%In particular, the authors in [29] prove as well that
%		\label{convergence minimizers husimi}
%		\forall  \vp \in   L^1(\R^d)  +L^\infty (\R^d)
%		\quad
%		\t{there holds}
%		\quad
%		\lim_{N\rightarrow \infty }
%		\< \vp , m_{N} \>
%		=
%		\< \vp ,  f_\rho  \> \ .
%As is well-known, the prescription given by \eqref{husimi def} turns Wigner transforms
%into measures. That is
%	0 \leq  m_N(x,p) \leq 1 \qquad \t{and}
%\int_{	\R^{2d}	} m_N(x,p) d x dp =1 \ .
%In what follows, for $s \geq 0 $,
% we introduce the Sobolev space $W_s $ of functions $m\in L_{loc}^1 (\R^{2d})$
% that have the finite norm
%|m|_s \equiv \sup_{ (\xi, \eta) \in \R^{2d}}  \<  \xi , \eta \>^{-s} |\hat m (\xi, \eta)| \  .
%The main advantage of introducing the space $W_s$
%is that one has the following formula at our disposal
%for the Wigner transforms $f_N:$
%	\hat f_N(\xi , \eta)
%	= \frac{1}{N}
%	\t{Tr}  \Big[
%	\O_{\xi ,\eta}  \gamma_N
%	\Big]
%	 \ , \qquad (\xi, \eta) \in \R^{2d} \  ,
%$\O_{\xi , \eta  } \equiv  \exp(  i x \cdot \xi + i p \cdot \eta  )$
%with $p \equiv - i \hbar \nabla_x$.
%In particular, while the $L^1$ norm of $f_N$ cannot be controled
%by trace norms of $\omega_N$, one now simply has
%| 		\hat f_N(\xi , \eta) | \leq  \frac{1}{N} \|	 \O_{\xi , \eta  }	\|   \t{Tr} \omega_N \leq 1  \ , \qquad
%(\xi, \eta) \in \R^{2d} \ .
%Thus, $ | 	 f_N	|_{s} \leq 1$ for all $N \geq 1$ and $s \geq 0 $.
%These  uniform bounds will replace the lack of $L^1$ boundeness in our analysis.
%\begin{proposition}[Convergence of states at zero temperature]
%The following statements hold true  (up to extraction of a subsequence)
%	\item For all $ p \in [1,2]$, $m_N$ converges  to $f$ strongly in $L^p$.
%	\item For all $s \geq 1$, $f_N$ converges to $f$ strongly in $W_s$.
%	\item
%If for all $N\geq1$,   $\gamma_N \equiv \gamma^{(1)}_{\Psi_N}$
%is an orthogonal projection, then
%	 $f_N$  converges to  $f$, strongly in $L^2$.
%	 Further, the following holds true for the $p=1 $ case
%	 \begin{equation}
%	 	\label{mass}
%	 	 \limsup_{N\rightarrow \infty } \|	 f - f_N	\|_{L^1}
%	 	  =
%	 	  	 	 \limsup_{N\rightarrow \infty }  \| f_N	\|_{L^1} -1 \ .
%	 \end{equation}
%If, additionally,  the following condition   is satisfied
%	\begin{equation}
%\label{sobolev type}
%		\sup_{N \geq 1 }
%		\sup_{z'\in \R^{2d}}
%		\frac{1}{|z'|} \int_{	\R^{2d}	}
%		|f_{N}  ( z  + \hbar z'   )  - f_N (z)  | \d z  = 0    \ ,
%	\end{equation}
%then  for all $ p \in [1,2]$,  $f_N$ converges to   $f$ strongly  in $L^p$.
%	\end{proposition}
%The condition given by \eqref{sobolev type}
%can be understood    in terms of integral of derivatives.
%Namely, one may easily replace it  with  the condition
%$\sup_{N \geq 1} \|	 \nabla f_N	\|_{L^1} <\infty$.
%Indeed, the Taylor formula
%$ f (z + \hbar z' ) - f(z) = \hbar z ' \int_0^1 \nabla f(z + t \hbar z' ) \d t $
%and a change of variables
%$z \mapsto z - t \hbar z' $
%	\begin{equation}
%	\label{sobolev type 2}
%	\sup_{z'\in \R^{2d}}
%	\frac{1}{|z'|} \int_{	\R^{2d}	}
%	|f_{N}  ( z  + \hbar z'   )  - f_N (z)  | \d z  \leq
%	\hbar \|	 \nabla f_N	\|_{L^1}     .
%The fact that we restrict ourselves to orthogonal projections
%$\gamma_N = \gamma_N^2$
%is a consequence of Lieb's principle.
%Namely, the minimizers of the ground state problem  for $E(N)$
%are \textit{always} given by Slater determinants.
%In view of Fatou's Lemma,
%$\liminf_{N\rightarrow \infty} \|	f_N\|_{L^1} \geq \|	 f	\|_{L^1} =1 $.
%Thus,  the  above Proposition     proves that
%	a sequence of orthogonal projections corresponding to minimzers  have their Wigner transforms converging in $L^1$ if and only if
%	\begin{equation}
%		 \lim_{N \rightarrow \infty}
%		 \|	f_N\|_{L^1}  =1 \ .
%	\end{equation}
%\begin{proof}[Proof of (1)]
%	Clearly, it is sufficient to establish the $p=1$ and $p=2$ cases, since everything in between follows from interpolation.
%	We treat these cases separately.
%	\vspace{3mm}
%	\underline{The $p=2$ case.}
%	Since the limiting measure $f $ is in $L^1(\R^{2d})$, we may
%	use $\vp  =  f$ in \eqref{convergence minimizers husimi}
%	to conclude that
%	\begin{equation}
%		\lim_{N \rightarrow \infty} \< f , m_N\>
%		=
%		 \<  f ,f  \> =1 \ .
%	\end{equation}
%Since the $L^2$ norm is generated by an inner product, we find that
%	\|	 m_N -f 	\|_{L^2}^2
%	 =
%	 \|	 f\|_{L^2}^2
%	  	 \|	  m_N\|_{L^2}^2
%	  	   -
%	  	   2 \< m_N , f\>
%  \end{equation}
%Since $\|	f\|_{L^2}=1$
%$\|	 m_N	\|_{L^2} \leq \| m_N	\|_{L^\infty} \|	m\|_{L^1} =1 $ we find that
%	\limsup_{N\rightarrow \infty }
%		\|	 m_N -f 	\|_{L^2}^2
%		\leq
%		2 - 2 \lim_{N \rightarrow \infty} 	  	   2 \< m_N , f\>
%		 =  0 \  .
%This finishes the proof of the $p=2$ case.
%\underline{The $p=1$ case}.
%We decompose the norm into two parts. One is controlled by the difference in $L^2$ norm, and the other one we use again \eqref{convergence minimizers husimi} with $\vp = 1 -f \in L^\infty(\R^{2d})$. That is, using that
%\|	 f - m_N	\|_{L^1}
%  =
%  \|	  f ( f - m_N	 ) \|_{L^1}
%  +
%  \|	(1-f )   m_N  	\|_{L^1}
%  \leq
%  \|	f	\|_{L^2} \|	 f - m_N	\|_{L^2}
%  + \< 1 -f ,  m_N \> \ .
%It suffices now to take the $N \rightarrow \infty$ limit, which finishes the $p=1$ case.
%\begin{proof}[Proof of (2)]
%Let $s \geq1$. 	We use the triangle inequality to find
%	\begin{equation}
% | f_N - f  |_s
% \leq
%  | f_N  - m_N|_s + | m_N - f|_s
% \leq
%  | f_N  - m_N|_s +  \|	m_N - f \|_{L^1}
%	\end{equation}
%where we used $  | \cdot |_s \leq \|	 \cdot 	\|_{L^1}$.
%The second term in the last displayed equation converges to zero, in view of part (1) of the Proposition.
%For the second term, we study its Fourier transform. Namely,
%	|  ( \hat f_N - \hat m_N )(\xi,\eta )|
%%	 =
%%	 |  ( \hat f_N  - \hat f_N \hat  {G} _\hbar )(\xi , \eta)  |
% =
%	| \hat f_N (\xi, \eta)|
%	| 1  -  \hat  {G}_1 ( \hbar (\xi , \eta)  )  |
%	\leq
%	C ( G_1)
%	  \hbar  \<  \xi , \eta \>
%		| \hat f_N (\xi, \eta)|
%		\leq
%			C ( G_1)
%		\hbar   \< \xi ,\eta \>
%where in the last line we used the uniform bound $\|  \hat f_N	\|_{L^\infty } \leq 1$.
%Here, the constant $C(G_1)$ is
%the Lipschitz constant for $\hat G_1$--recall that $G_1$ is a Schwartz function.
%The last inequality now implies that for any $s \geq 1$
%|  f_N  - m_N |_s
% \sup_{\xi ,\eta \in \R^d }
% \< \xi , \eta  \>^{-1}
% |  ( \hat f_N - \hat m_N )(\xi,\eta )|
% \leq C \hbar \ .
%This finishes the proof
%\begin{proof}[Proof of (3)]
%	First, we employ the fact that (up to scaling), the Wigner transformation is a unitary map between $L_{x,x'}^2$ and $L_{x,p}^2$.
%	In our setting, one is able to conclude that
%	\|	f_N	\|_{L^2}
%	=
%	 \hbar^{d/2}
%	 \|	 \gamma_{N} 	\|_{L^2}
%	 =
%	 \hbar^{d/2}
%	 \|	 \gamma_{N} 	\|_{HS}
%	 =
%	 \hbar^{d/2}
%	 \t{Tr}(\gamma^*_N \gamma_N)^{\frac{1}{2}}
%	 =
%	 \hbar^{d/2}
%	 \t{Tr}( \gamma_N)^{\frac{1}{2}}
%	 = (\hbar^d N )^{1/2} =1 \ .
%Consequently, $\|	f_N	\|_{L^2} =1 $ for all $N \geq1$.
%Thus, we can extract a subsequence such that $f_N$ converges weakly in $L^2$.
%Of course,  in view of \eqref{convergence minimizers}, this limit is given by $f = f_\rho$. Further, we note that because $f$ is a characteristic function, it holds true that $f = f^2$. Hence, we conclude
%	\|		 f  \|_{L^2} = \|	 f^2	\|_{L^1}^{1/2} = \|	f 	\|_{L_1}^{1/2} =1 \ ,
%because $\|	f\|_{L^1} =  \int_{	\R^d	} \rho(x) \d x = 1 $.
%Thus, $\lim_{N\rightarrow \infty} \|f_N	\|_{L^2}  =  \|	f\|_{L^2}$.
%Since $L^2$ is a uniformly convex space, we may conclude that convergence is strongly in $L^2$. This finishes the $p=2$ case.
%Let us now prove the formula \eqref{mass}.
%Again, we use the fact that $f=f^2$ because $f$ is a characteristic function.
%We estimate the $L^1$ norm difference as follows
%\label{L1 step 1}
%\|	f - f_N	\|_{L^1}
%\int_{	\R^{2d}	} f | f - f_N|
%\int_{	\R^{2d}	} (1-f) | f - f_N|
%\int_{	\R^{2d}	} f | f - f_N|
%\int_{	\R^{2d}	} (1-f) | f_N|
%The first term in \eqref{L1 step 1} can be estimates in view of the Cauchy-Schwarz inequality, and the $p=2$ case. Namely
%	\int_{	\R^{2d}	} f | f - f_N|
%	 \leq
%	  \|	f 	\|_{L^2}
%	  \|	f - f_N	\|_{L^2}  = 	  \|	f - f_N	\|_{L^2}  \ .
%The second one is more tricky.
%First, because of the $p=2$ case, we may assume that convergence holds pointwise almost everywhere.
%Hence, Fatou's Lemma implies  that
%	  1 = \int_{	\R^{2d}	} f |f|
%	  \leq \liminf_{N\rightarrow \infty}
%	  \int_{	\R^{2d}	} f |f_N | \ .
%On the other hand, we know that $ \|	f_N\|_{L^2}=1$ and thus
%the Cauchy-Schwarz inequality implies
%\limsup_{N \rightarrow \infty}
%	\int_{\R^{2d}} f |f_N|
%	 \leq
%	 \| f \|_{L^2}
%	 \limsup_{N \rightarrow \infty}
%	  \|	f_N	\|_{L^2} =1   \ .
%Thus, $\lim_{N \rightarrow \infty } \int_{\R^{2d}} f |f_N| =1$.
%Thus, going back to \eqref{L1 step 1} we see that
%	\limsup_{N \rightarrow \infty} \|	 f - f_N	\|_{L^1}
%	\leq
%	 	  	 	 \limsup_{N\rightarrow \infty }  \| f_N	\|_{L^1} -1 \ .
%The proof  of equality
%is finished one we combine the above inequality
%with the estimate
% $\|f_N\|_{L^1} -    1
%  \leq
%  \| f \|_{L^1 } +
%\| f - f_N\|_{L^1}
%= -1
%\| f - f_N\|_{L^1},$
%where we used $\int f = 1$.
%Let us now show that under the additional condition \eqref{sobolev type},
%strong convergence holds for $p=1$ as well.
%One may then use interpolation to conclude the same thing for $ p \in (1,2)$.
%We compare the distance with respect to the Husimi measure.
%Namely, the triangle inequality gives
%	 \|	  f_N 	\|_{L^1 }
%	 \leq
%	 	 \|	    m_N 	\|_{L^1 }
%	 	 +
%	 	 	 \|	  m_N - f	\|_{L^1 }
%	 	 	   =
%	 	 	   1  +	 	 	 \|	  m_N - f	\|_{L^1 }
%Let us now analyze the first term. We write
%		 	 f_N (z ) - m_N(z)
%%		 	  =
%%		 	  \int_{	\R^{2d}	}
%%		 	   \big(
%%		 	   f_N(z) - f_N(z')
%%		 	   \big)
%%		 	   \mathscr{G}_\hbar (z - z') \d z'
%		 	    =
%		 	     \int_{	\R^{2d}	}
%		 	    \big(
%		 	    f_N(z  ) - f_N( z  - \hbar  z' )
%		 	    \big)
%		 	    \mathscr{G}_1 (z') \d z'
% which then implies that
% \begin{align}
% 		 	 \|	  f_N - m_N 	\|_{L^1 }
% 		&  	 \leq
% 		 	 \int_{	\R^{2d}	 \times \R^{2d}}
%|z' |^{-1}
% 		 	 f_N(z  ) - f_N( z  - \hbar  z' )
% |
% 	|z'|	 	 \mathscr{G}_1 (z') \d z' \d  z  , \\
% 		 	 & \leq
%  \  		 	\|	  |z ' | \mathscr{G}_1	\|_{L^1} \
% 		 	\sup_{z'\in \R^{2d}}
% 		 	\frac{1}{|z'|} \int_{	\R^{2d}	}
% 		 	|f_{N}  ( z  + \hbar z'   )  - f_N (z)  | \d z  \ .
% \end{align}
%In view of the assumption of the Proposition, this finishes the proof.

%\underline{The $p=1$ case}.
%Again, we use the fact that $f=f^2$ because $f$ is a characteristic function.
%We estimate the $L^1$ norm difference as follows
%\label{L1 step 1}
%\|	f - f_N	\|_{L^1}
%\int_{	\R^{2d}	} f | f - f_N|
%\int_{	\R^{2d}	} (1-f) | f - f_N|
%\int_{	\R^{2d}	} f | f - f_N|
%\int_{	\R^{2d}	} (1-f) | f_N|
%The first term in \eqref{L1 step 1} can be estimates in view of the Cauchy-Schwarz inequality, and the $p=2$ case. Namely
%	\int_{	\R^{2d}	} f | f - f_N|
%	 \leq
%	  \|	f 	\|_{L^2}
%	  \|	f - f_N	\|_{L^2}  = 	  \|	f - f_N	\|_{L^2}  \ .
%The second one is more tricky.
%Heuristically, if one thinks of $1 -f  $
%as the projection outside of the set
%$\mathcal O  = \{ (x,p) \in \R^{2d} : |p |^2 \leq C_{TF} \rho(x)^{2/d} \}$,
%then the second term should be zero in the limit $N \rightarrow \infty$.
%Let us now check this is true.
%First, because of the $p=2$ case, we may assume that convergence holds pointwise almost everywhere.
%Hence, Fatou's Lemma implies  that
%	  1 = \int_{	\R^{2d}	} f |f|
%	  \leq \liminf_{N\rightarrow \infty}
%	  \int_{	\R^{2d}	} f |f_N | \ .
%On the other hand, we know that $ \|	f_N\|_{L^2}=1$ and thus
%the Cauchy-Schwarz inequality implies
%\limsup_{N \rightarrow \infty}
%	\int_{\R^{2d}} f |f_N|
%	 \leq
%	 \| f \|_{L^2}
%	 \limsup_{N \rightarrow \infty}
%	  \|	f_N	\|_{L^2} =1   \ .
%Thus, $\lim_{N \rightarrow \infty } \int_{\R^{2d}} f |f_N| =1$.

%\section{From Hartree to Vlasov: Proofs}
%Let us sketch the proof of Theorem \ref{theorem convergence 2}  --  it follows very closey that of (BPSS2016) Theorem 2.5,
%but including the boson influence.
%Let us state the auxiliary lemmata that we will need. In order to do so, we introduce the following objects.
%\textit{Regularization of initial data}.
%Let $f_0^\hbar = W^\hbar [\omega_0^\hbar ]$
%be as in the statement of  Theorem \ref{theorem convergence 2},
%and consider
%for $(x,p) \in \R^6$
%$\ve \in (0,1)$
%the following mollifier
%	g_1 (x,p) : =  \frac{1}{ (2\pi)^3  }
%	e^{  -  \frac{1}{2} (x^2 + p^2 )}
%	\qquad
%	\t{and}
%	\qquad
%	g_\ve (x,p ) : = \ve^{ -3 } g_1( \ve^{-1/2} x ,   \ve^{-1/2}  p )
%together with the regularized initial data
%	f_0^{\ve, \hbar }
%	:=    f_0^\hbar  *  g_\ve
%		 \in \S_{x,p} ( \R^3 \times \R^3 ) \  .
%We consider the Weyl quantization of the regularized initial data  defined above, namely
%	\omega_0^{ \ve , \hbar } : = \mathrm{Op}^\hbar [ f_0^{\ve, \hbar }   ]  \ .
%It turns out that $\omega_0^{\ve, \hbar }$ is again trace-class, with
%$0 \leq 	\omega_0^{ \ve , \hbar }  \leq 1$
%\mathrm{Tr} 	\omega_0^{ \ve , \hbar }  = N .
%We consider the solution to the coupled Hartree equation, corresponding to this regularized initial data.
%Namely we define.
%	We denote by
%	$  (  \omega^{\ve , \hbar }   ,    \vp^{\ve, \hbar }) $
%	the solution of the  coupled  Hartree equation \eqref{coupled hartree},
%	with  regularized initial data
%	$  (\omega_0^{\ve , \hbar } , \vp_0^\hbar)  .  $
% We will compare the fermionic densities $\omega^\hbar$ and $\omega^{\ve, \hbar}$
% with solutions to the Vlasov equation, but with modified initial data;
% either $f_0^\hbar$ or $f_0^{\ve, \hbar}$.
% We pull back this solutions to fermionic densities by means of the Weyl quantization.
% Notation-wise, this will be denoted by a tilde $\tilde  . $
%We denote by  $ (\widetilde f^\hbar ,  \widetilde  x^\hbar )$  the solution of the Vlasov-Newton \eqref{vlasov}
%equation with initial data
%$ ( f_0^\hbar  , x_0 , p_0)$,
%and   we denote  its Weyl quantization
%$ \widetilde \omega^{\hbar } : = \mathrm{Op}^\hbar [   \widetilde f^\hbar   ] $.
%	 We denote by  $ (\widetilde f^{\ve, \hbar} ,  \widetilde  x^{\ve, \hbar} )$
%	  the solution of the Vlasov-Newton \eqref{vlasov}
%	 equation with initial data
%	 $ ( f_0^{\ve, \hbar}  , x_0 , p_0)$,
%	 and   we denote  its Weyl quantization
%	 by
%	 $ \widetilde \omega^{\ve  , \hbar } : = \mathrm{Op}^\hbar [   \widetilde f^{\ve, \hbar  }   ] $.
%\noindent \textit{The 1-particle Hamiltonian}
%Let $\hbar > 0 $
%and consider the solution $(\omega^\hbar , \vp^\hbar)$ of the coupled Hartree equation.
%We denote the fermionic 1-particle Hamiltonian by
%	h_f^\hbar (t) : = - \frac{\hbar^2 }{2 } \Delta_x  +  \big( v * \rho_b^\hbar(t)  \big)  (x)
%	\qquad
%	 \t{ on }
%	 L^2 (\R^3 )
% with the bosonic density $\rho_b^\hbar (t, x) = |\vp^\hbar (t,x)|^2$.
% In what follows,
% we shall consider the unitary flow that it generates.
% Namely,
% we shall call the \textit{fermion flow}
% the unitary propagator solving
% \begin{equation}
%\label{fermion flow}
% 	\begin{cases}
%&  i \hbar \partial_t \U_\hbar  (t,s ) = h_f^\hbar (t) \U_\hbar  (t,s) \\
% 	& \U_\hbar   (t,t) = \1  \ .
% 	\end{cases}
% \end{equation}
%%\begin{lemma}[BPSS2016 Lemma 4.1]
%\begin{lemma}[BPSS2016 Lemma 4.2]
%	Assume that $ \| v \|_{W^{ 2 ,1 }}  < \infty$
%	and let $ \U_\hbar (t,s) $
%	be the fermion flow defined in \eqref{fermion flow}.
%	Then, there exists $C>0$ such that
%	for all $(\xi , \eta ) \in \R^3 \times \R^3 $
%	and $t,s \in \R $
%	there holds
%	\sup\Biggl\{
%	\frac{1}{ | k | }
%	\Big|
%	\mathrm{Tr}
%	\big[
%	e^{ i k \cdot x }, \U_\hbar^* (t,s) e^{i \xi \cdot x + i \eta \cdot p }
%	\U_\hbar (t,s)
%	\big]
%	\omega
%	\Big|
%	 \ : \
%	 k \in \R^d \ ,
%	 \mathrm{Tr}|\omega| \leq 1
%	\Biggr\}
%	\leq
%	\hbar ( | \xi| + |\eta |)
%	e^{ C  |t -s |   }
%where $ p = - i \hbar \nabla_x$ on $L^2 (\R^3 ). $
%\begin{lemma}[BPSS2016 Lemma 5.1]
%		Assume that $ \| v \|_{W^{ 2 ,1 }}  < \infty$
%	and let $ \U_\hbar (t,s) $
%	be the fermion flow defined in \eqref{fermion flow}.
%	Then,
%	under the assumptions of Theorem \ref{theorem convergence 2},
%	there exists $C>0$ such that
%	for all $(\xi , \eta ) \in \R^3 \times \R^3 $
%	and $t,s \in \R $
%	there holds
%	\begin{align}
%		\sup\Biggl\{
%		\Big|    \,
%		\mathrm{Tr}
%		\big[ x_i ,  \big[ x_ j
%	, \U_\hbar^* (t,s) e^{i \xi \cdot x + i \eta \cdot p }
%		\U_\hbar (t,s)
%		\big] \big ]
%		\omega
%	 \, 	\Big|
%		\ : \
%	i,j = 1 , \ldots, d , \
%	 & 	\mathrm{Tr}|\omega| \leq 1
%		\Biggr\} \\
%& \leq C \hbar^2
% ( | \xi| + |\eta |)^2
%		e^{ C  |t -s |   }			\nonumber
%	\end{align}
%	where $ p = - i \hbar \nabla_x$ on $L^2 (\R^3 ). $

\begin{thebibliography}{10}

[1]
R. Adami,
F. Golse,
A. Teta.
\textit{Rigorous derivation of the cubic NLS in dimension one. }
J. Stat. Phys. 127 (2007), no. 6, 1193–1220.

[2]
Z. Ammari,
F. Nier.
\textit{Mean field propagation of Wigner measures and BBGKY hierarchy for general bosonic states. }
J. Math. Pures Appl. (9) 95 (2011), no. 6, 585–626.

[3]
A. 	Athanassoulis,
T. Paul,
F. Pezzotti,
M. Pulvirenti.
\textit{			Strong semiclassical approximation of Wigner functions for the Hartree dynamics. }
Rend. Lincei Mat. Appl. 22, 525–552 (2011)

[4]
C. Bardos,
F. Golse,
A. D. Gottlieb,
N. J. Mauser.
Mean field dynamics of fermions and the time-dependent Hartree-Fock equation.
J. Math. Pures Appl. (9) 82 (2003), no. 6, 665–683.

[5]
N. 			Benedikter.
\textit{Effective Dynamics of Interacting Fermions from Semiclassical Theory to the Random Phase Approximation. }

[6]
N, 			Benedikter,
Desio. D.
\textit{Two Comments on the Derivation of the Time-Dependent Hartree-Fock Equation.}

[7]
N. 		Benedikter,
V. 	 Jakšić,
M. 			  Porta,
C.		  Saffirio,
B. 			   Schlein,
\textit{Mean-Field Evolution of Fermionic Mixed States. }
Comm. Pure Appl. Math., 69: 2250-2303 (2016).

[8]
N. Benedikter,
M. Porta,
C. Saffirio,
B. Schlein.
\textit{From the Hartree Dynamics to the Vlasov Equation. }
Arch Rational Mech Anal 221, 273–334 (2016).

[9]
N. Benedikter,
M. Porta,
B. Schlein.
\textit{			Mean–Field Evolution of Fermionic Systems. }
Commun. Math. Phys. 331, 1087–1131 (2014).

[10]
N. Benedikter,
M. Porta,
B. Schlein.
\textit{Mean-field dynamics of fermions with relativistic dispersion. }
J. Math. Phys. 1 February 2014; 55 (2): 021901.

[11]
N. Benedikter,
M. Porta,
B. Schlein.
\textit{Effective Evolution Equations from Quantum Dynamics.}
vol.7. Springer, Berlin (2016).

[12]
A. Bove,
G. D. Prato,
G. Fano.
\textit{An existence proof for the Hartree-Fock time-dependent problem with bounded two-body interaction.}
Communications in Mathematical Physics 37 (1974): 183-191.

[13]
E. Cárdenas.
\textit{Norm convergence of confined fermionic systems at zero temperature}.
In preparation.

[14]
T. Chen,
C. Hainzl,
N. Pavlovi\'c,
R. Seiringer.
\textit{Unconditional Uniqueness for the Cubic Gross-Pitaevskii Hierarchy via Quantum de Finetti.} Comm. Pure Appl. Math (2015), \textbf{68}: 1845-1884. https://doi.org/10.1002/cpa.21552

[15]
L. Chen,
J. Oon Lee.
\textit{Rate of convergence in nonlinear Hartree dynamics with factorized initial data.}
J. Math. Phys. 1 May 2011; 52 (5): 052108

[16]
L. Chen,
J. Oon Lee,
B. Schlein.
\textit{Rate of convergence towards Hartree dynamics.}
J. Stat. Phys. 144 (2011), no. 4, 872–903.

[17]
J. Chong,
L. Lafleche,
C. Saffirio.
\textit{From many-body quantum dynamics to the Hartree-Fock and Vlasov equations with singular potentials}. 2021. Preprint: https://arxiv.org/abs/2103.10946

[18] %This paper is the "first" to demonstrate a doubly degenerate BF mixture. The instability due to energy scales is pointed out here.
B. DeSalvo,
K. Patel,
J. Johansen,
C. Chin.
\textit{Observation of a Degenerate Fermi Gas Trapped by a Bose-Einstein Condensate.} Phys. Rev. Lett. \textbf{119}, 233401 (2017,12), https://link.aps.org/doi/10.1103/PhysRevLett.119.233401

[19]
C. Dietze,
J. Lee.
\textit{Uniform in Time Convergence to Bose-Einstein Condensation for a Weakly Interacting Bose Gas with an External Potential. }

[20] %Here there are heavy bosons and light fermions. Good reference for physical motivation about caring for interspecies interactions.
B. J. DeSalvo,
K. Patel,
G. Cai,
C. Chin.
\textit{Observation of fermion-mediated interactions between bosonic atoms.}
Nature \textbf{568}, 61–64 (2019). https://doi.org/10.1038/s41586-019-1055-0

[21]
M. Duda,
XY. Chen,
A. Schindewolf,
R. Bause,
J. Milczewski,
R. Schmidt,
I. Bloch,
XY. Luo.
\textit{Transition from a polaronic condensate to a degenerate Fermi gas of heteronuclear molecules.} Nat. Phys. \textbf{19}, 720–725 (2023). https://doi.org/10.1038/s41567-023-01948-1

[22]
A. Elgart,
L. Erd\"os,
B. Schlein,
H.-T. Yau.
\textit{Nonlinear Hartree equation as the mean field limit of weakly coupled fermions. }
J. Math. Pures Appl. (9) 83 (2004), no. 10, 1241–1273.

[23]
L. Erd\"os,
B. Schlein,
H.-T. Yau.
\textit{Derivation of the cubic non-linear Schro dinger equation from quantum dynamics of many-body systems. }
Inv. Math. 167 (2006), 515– 614.

[24]
L. Erd\"os,
B. Schlein,
H.-T. Yau.
\textit{Rigorous derivation of the Gross-Pitaevskii equation. }
Phys. Rev. Lett. 98 (2007), no. 4, 040404.

[25]
L. Erd\"os,
B. Schlein,
H.-T. Yau.
\textit{Rigorous derivation of the Gross-Pitaevskii equa- tion with a large interaction potential. }
J. Amer. Math. Soc. 22 (2009), 1099–1156.

[26]
L. Erd\"os,
B. Schlein,
H.-T. Yau.
\textit{Derivation of the Gross-Pitaevskii equation for the dynamics of Bose-Einstein condensate, }
Ann. of Math. (2) 172 (2010), no. 1, 291–370.

[27]
L. Erd\"os,
H.-T. Yau.
\textit{Derivation of the nonlinear Schro dinger equation from a many body Coulomb system. }
Adv. Theor. Math. Phys. 5 (2001), no. 6, 1169–1205.

[28] %dually superfluid bose-fermi mixture.
I. Ferrier,
M. Delehaye,
S. Laurent,
A. T. Grier,
M. Pierce,
B. S. Rem,
F. Chevy,
C. Salomon.
\textit{A mixture of Bose and Fermi superfluids.} Science, \textbf{345},1035-1038(2014).DOI:10.1126/science.1255380

[29]
S. Fournais,
M. Lewin,
J. P. Solovej.
\textit{The semi-classical limit of large fermionic systems}.
Calc. Var. (2018) 57:105

[30]
L. Fresta,
M. Porta,
B. Schlein.
\textit{Effective Dynamics of Extended Fermi Gases in the High-Density Regime. }
Commun. Math. Phys. 401, 1701–1751 (2023).

[31]
J. Fr\"ohlich,
A. Knowles.
\textit{A microscopic derivation of the time-dependent Hartree- Fock equation with Coulomb two-body interaction, }
J. Stat. Phys. 145 (2011), no. 1, 23–50.

[32]
J. Fr\"ohlich,
A. Knowles,
S. Schwarz.
\textit{On the mean-field limit of bosons with Coulomb two-body interaction. }
Comm. Math. Phys. 288 (2009), no. 3, 1023–1059.

[33]
I. 			Gasser,
R. 			Illner,
P.A. 			Markowich,
C. 			 Schmeiser.
\textit{			Semiclassical, {   	$t \rightarrow  \infty$}
asymptotics and dispersive effects for HF systems. }
Math. Model. Numer. Anal. 32, 699–713 (1998)

[34]
J. Ginibre,
G. Velo.
\textit{The classical field limit of scattering theory for non-relativistic many-boson systems. I. }
Commun. Math. Phys. 66, 37–76 (1979).

[35]
J. Ginibre,
G. Velo.
\textit{The classical field limit of scattering theory for non-relativistic many-boson systems. II.}
Commun. Math. Phys. 68 (1979), 45–68.

[36]
F. Golse,
C. Mouhot,
T. Paul.
\textit{On the Mean Field and Classical Limits of Quantum Mechanics. }
Communications in Mathematical Physics, 343(1):165–205, Apr. 2016.

[37]
F. Golse,
T. Paul.
\textit{			The Schr odinger Equation in the Mean-Field and Semiclassical Regime.}
Arch. Rational Mech. Anal. 223 (2017), 57–94.

[38]
F. Golse,
T. Paul.
\textit{Semiclassical Evolution With Low Regularity. }
Journal de Mathématiques Pures et Appliquées, 151:257–311, July 2021.

[39]
F. Golse,
T. Paul.
\textit{Optimal transport pseudo metrics for quantum and classical densities.}
Journal of Functional Analysis, 282(9):109417, May 2022.

[40]
M. Grillakis,
M. Machedon,
D. Margetis.
\textit{Second-order corrections to mean field evolution of weakly interacting bosons. I. }
Comm. Math. Phys. 294 (2010), no. 1, 273– 301.

[41]
M. Grillakis,
M. Machedon,
D. Margetis,
\textit{Second-order corrections to mean field evolution of weakly interacting bosons.
II. }
Adv. Math. 228 (2011), no. 3, 1788–1815.

[42]
M. Grillakis,
M. Machedon.
\textit{Pair excitations and the mean field approximation of interacting bosons, I. }
ArXiv:1208.3763 [math-ph] (2012).

[43]
K. Hepp.
\textit{The classical limit for quantum mechanical correlation functions.}
Commun. Math. Phys. 35 (1974), 265–277.

[44]
K. Kirkpatrick,
B. Schlein,
G. Staffilani.
\textit{Derivation of the two dimensional nonlinear Schro  dinger equation from many body quantum dynamics.}
Amer. J. Math. 133 (2011), no. 1, 91-130.

[45]
A. Knowles,
P. Pickl.
\textit{Mean-field dynamics: singular potentials and rate of convergence. }
Comm. Math. Phys. 298 (2010), no. 1, 101–138.

[46]
L. Lafleche.
\textit{			Quantum optimal transport and weak topologies.}

[47]
S. Fournais,
M. Lewin,
J. P. Solovej.
\textit{The semi-classical limit of large fermionic systems}.
Calc. Var. (2018) 57:105 https://doi.org/10.1007/s00526-018-1374-2.

[48]
P. L.  Lions,
T. Paul.
\textit{Sur les mesures de Wigner.}
Rev. Mat. Iberoam. 9, 553–618(1993)

[49]
P.A. Markowich,
N.J. 			Mauser,
\textit{The classical limit of a self-consistent quantum Vlasov equation. }
Math. Models Methods Appl. Sci. 3(1), 109–124 (1993)

[50]
P. Mocz,
L. Lancaster,
A. Fialkov,
F. Becerra,
P.-H. Chavanis.
\textit{Schrödinger-Poisson–Vlasov-Poisson correspondence}
Phys. Rev. D 97, 083519 – Published 25 April 2018.

[51] %This is the original theoretical reference for Bose-Fermi Mixtures
K. Mølmer.