\section{Introduction}

In this paper, we consider the dynamics of a system of $N$ bosons living in $\mathbb{R}^d$ in the so-called \emph{mean-field regime} and we are interested in the time-dependent fluctuations around the Hartree state solution.

The $N$-body system is described by a time-dependent wave function $\Psi_N(t,x_1,...,x_N)$ in the symmetric space $\gH^N=\bigotimes_{\text{sym}}^N L^2(\mathbb{R}^d)$ and its evolution is governed by Schr\"odinger's equation
\begin{eqnarray} \label{eq:Schrodinger-dynamics}
\left\{ \begin{gathered}
 i\, \dot{\Psi}_{N}(t) = H_N \Psi_{N}(t),  \hfill \\
  \Psi_{N}(0)=\Psi_{N,0},\hfill \\
\end{gathered}  \right.\end{eqnarray}
or, equivalently,
$$\Psi_N(t)=e^{-it H_N}\Psi_N(0).$$
The many-particle Hamiltonian $H_N$ is given by
\begin{equation}
H_N= \sum\limits_{j = 1}^N (-\Delta)_{x_j} + \frac{1}{N-1} \sum\limits_{1 \le j < k \le N} {w(x_j-x_k)},
\label{eq:def_H_N}
\end{equation}
where the Laplacian $-\Delta_{x_j}$ describes the kinetic energy of the $j$-th particle and $w:\mathbb{R}^d \to \mathbb{R}$ is a measurable, even function describing the interactions between the particles.

The fact that we are considering the mean-field regime is apparent in the factor $1/(N-1)$ which makes the one-body term and the two-body interaction term of order $N$ in the Hamiltonian $H_N$. We could as well take $1/N$ instead of $1/(N-1)$ but the latter choice simplifies some expressions.

We shall always assume that the interaction potential $w$ satisfies the operator inequality
\begin{eqnarray} \label{eq:asp-w} w^2 \le C(1-\Delta)
\end{eqnarray}
on $L^2(\mathbb{R}^d)$, for some constant $C>0$. The condition (\ref{eq:asp-w}) in particular ensures that $w$ is relatively form bounded by the Laplacian with bound as small as we want:
$$ |w| \le \sqrt{C(1-\Delta)} \le \varepsilon (1-\Delta)+ C \varepsilon^{-1}$$
for every $\varepsilon>0$. The reader may keep in mind the example of $w$ being the Coulomb potential $|x|^{-1}$ in $\mathbb{R}^3$. We may also consider bosons living in a bounded domain of $\mathbb{R}^d$ (with appropriate boundary conditions), or replace the Laplacian $-\Delta$ by the Schr\"odinger operator $-\Delta+V(x)$ with an external potential $V(x)$, or the pseudo-relativistic counterpart $\sqrt{1-\Delta}$. In these cases our results still apply without significant changes (see Remark \ref{rm:ext}).

Under the assumption (\ref{eq:asp-w}), the Hamiltonian $H_N$ is bounded from below. Therefore, it can be defined as a self-adjoint operator (still denoted by $H_N$) by Friedrichs' method ~\cite{ReeSim2}, and this gives a proper meaning to the many-particle Schr\"odinger equation~\eqref{eq:Schrodinger-dynamics}, by Stone's theorem.

Bosons have the ability to undergo {\em Bose-Einstein condensation}, which means that a macroscopic number of the particles live in a common quantum state $u(t)\in L^2(\mathbb{R}^d)$. It is a very important fact that (complete) Bose-Einstein condensation is stable under the Schr\"odinger flow in the limit $N\to\ii$. More precisely, if the initial datum is a pure Hartree state
$$\Psi_{N,0}(x_1,...,x_N)=(u_0)^{\otimes N}(x_1,...,x_N)=u_0(x_1)\cdots u_0(x_N),$$
then the many-particle wave function is well described by a Hartree state
$$\Psi_{N}(t)\approx u(t)^{\otimes N}$$
for all times $t\ge 0$ in the limit $N\to\ii$, in a sense to be made precise below (and, in particular, \emph{not} in norm!). Heuristically, the mean-field potential experienced by each particle can be approximated by $|u(t)|^2\ast w$. This observation leads to the nonlinear Hartree equation for the condensate wave function $u(t)$:
\begin{eqnarray} \label{eq:Hartree-equation}
\left\{ \begin{gathered}
 i\, \dot u(t) =  \big(-\Delta +|u(t)|^2*w -\mu(t)\big) u(t),  \hfill \\
  u(0)=u_{0}.\hfill \\
\end{gathered}  \right.
\end{eqnarray}
Here the gauge parameter $\mu(t)\in \mathbb{R}$ can be freely chosen and different choices lead to different phases of $u(t)$. In this paper it will be convenient to choose
\begin{eqnarray} \label{eq:def-mu}
\mu(t):=\frac12\iint_{\mathbb{R}^d\times\mathbb{R}^d}|u(t,x)|^2w(x-y)|u(t,y)|^2\,dx\,dy
\end{eqnarray}
which implies the compatibility of the energies:
\begin{align*}
N \langle u(t), (-\Delta +|u(t)|^2*w/2) u(t) \rangle \approx
\langle \Psi_{N}(t), H_N \Psi_{N}(t)\rangle = \langle i \dot{\Psi}_{N}(t), \Psi_{N}(t)\rangle \\
\approx N  \langle i \dot{u}(t), u(t) \rangle = N  \langle {u}(t), (-\Delta+|u(t)|^2*w -\mu(t))  u(t) \rangle.
\end{align*}
By the usual argument based on Duhamel's formula, one can show that for every initial datum $u_0\in H^1(\mathbb{R}^d)$, the Hartree equation (\ref{eq:Hartree-equation}) admits a unique solution
$$u(t)\in C^0([0,\infty),H^1(\mathbb{R}^d))\cap C^1([0,\infty),H^{-1}(\mathbb{R}^d))$$
with the gauge constant $\mu(t)$ given in (\ref{eq:def-mu}). Moreover, the norm $\|u(t)\|_{L^2}$ and the energy $\left \langle u(t), \left( -\Delta + |u(t)|^2*w/2\right) u(t) \right \rangle$ are constant in time.

In fact, the approximation $\Psi_{N}(t)\approx u(t)^{\otimes N}$ holds in the topology of {\em reduced density matrices}. For every $k\geq1$, the $k$-particle density matrix $\gamma_{\Psi_{N}(t)}^{(k)}$ of $\Psi_{N}(t)$ is the trace class operator on the $k$-particle space $\mathfrak{H}^k=\bigotimes_{\rm sym}^kL^2(\mathbb{R}^d)$ with kernel
\[\begin{gathered}
\gamma_{\Psi_{N}(t)}^{(k)}(x_1,...,x_k; y_1,...,y_k)\hfill \\
   \quad \quad = \int \Psi_{N}(t,x_1, ... , x_N) \overline{\Psi_{N}(t,y_1,...,y_k, x_{k+1}, ... , x_N)}\, {\rm d} x_{k+1}... {\rm d} x_{N}.\hfill \\
\end{gathered} \]
It is known (see, for example, ~\cite{Spohn-80,ErdYau-01,RodSch-09,KnoPic-10,CheLeeSch-11}) that if $\Psi_N(0)=u(0)^{\otimes N}$, then the $k$-particle density matrix of $\Psi_{N}(t)$ converges to that of the Hartree state $u(t)^{\otimes N}$, that is,
\begin{eqnarray} \label{eq:cv-density-matrices}
\lim_{N\to \infty} \Tr_{\mathfrak{H}^k} \Big| \gamma^{(k)}_{\Psi_{N}(t)}- |u(t) \rangle \langle u(t)|^{\otimes k} \Big| =0,
\end{eqnarray}
for every $k\in \mathbb{N}$ and every $t\geq0$.

We mention that the convergence (\ref{eq:cv-density-matrices}) does not imply the corresponding approximation in the norm topology of $\mathfrak{H}^N$ for the wave function $\Psi_{N}(t)$. We will see later that even if $\Psi_N(0)=u(0)^{\otimes N}$, then for every $t>0$ fixed, $\Psi_{N}(t)$ {\em never} stays close to the Hartree state $u(t)^{\otimes N}$ in the norm of $\mathfrak{H}^N$, except in the non-interacting case $w\equiv0$.

The aim of our paper is to derive an effective equation which gives the exact behavior of the wave function $\Psi_{N}(t)$ in the norm topology as $N\to\ii$. More precisely, we shall determine continuous mappings $t\mapsto \phi_k(t)\in \bigotimes_{\rm sym}^k\mathfrak{H}_+(t)$, where $\mathfrak{H}_+(t)=\{u(t)\}^\bot \subset \mathfrak{H}$, such that for all times $t \ge 0$,
\begin{equation}
\boxed{\lim_{N\to\ii}\norm{\Psi_N(t)-\sum_{k=0}^N u(t)^{\otimes (N-k)}\otimes_s\phi_k(t)}_{\mathfrak{H}^N}=0.}
\label{eq:fluctuations-intro}
\end{equation}
The precise statement is given in Theorem \ref{thm:fluctuations} below. The evolution of the family of functions $(\phi_k(t))_{k=0}^\infty$ is governed by an effective quadratic Hamiltonian in Fock space $\mathcal{F}=\bigoplus_{n\ge 0} \mathfrak{H}^n$, which is derived from Bogoliubov's approximation, as will be explained in detail in Section \ref{sec:Bogoliubov}. The convergence (\ref{eq:fluctuations-intro}) is much more precise than (\ref{eq:cv-density-matrices}) and it implies (\ref{eq:cv-density-matrices}) easily (see Corollary~\ref{cor:Hartree} below).

The behavior of the many-particle wave function in the mean-field regime has been studied in several situations but, as far as we know, never close to a Hartree state $u(t)^{\otimes N}$ in the $N$-body space $\mathfrak{H}^N$ like here. The only exception is \cite{BenKirSch-11}, where the fluctuations around the Hartree evolution of factorized initial data are shown to satisfy a central limit theorem. Most previous works focus instead on the description of the fluctuations around a \emph{coherent state in Fock space}, see for example ~\cite{Hepp-74,GinVel-79,GinVel-79b,GriMacMar-10,GriMacMar-11,Chen-12,GriMac-12,BenOliSch-12}. In the mean-field limit, the evolution of a coherent state is also governed by the {\em same} nonlinear Hartree equation~\eqref{eq:Hartree-equation} in the topology of density matrices, but the Bogoliubov Hamiltonian describing the fluctuations is {\em different} from the Hartree case, as we shall explain in Section~\ref{sec:coherent}.

Note that, in the time-independent setting, a convergence similar to (\ref{eq:fluctuations-intro}) was recently established for eigenvectors of $H_N$ by Lewin, Nam, Serfaty and Solovej in \cite{LewNamSerSol-13}. In fact, some tools from~\cite{LewNamSerSol-13} will be used in this work. This includes, in particular, the unitary operator $U_N(t)$ defined later and a quantitative estimate in truncated Fock spaces. The convergence (\ref{eq:fluctuations-intro}) will then follows from some energy estimates on the Bogoliubov dynamics and on the Schr\"odinger dynamics.

The paper is organized as follows. In Section \ref{sec:main-result} we state precisely our main result. In Section \ref{sec:coherent} we quickly discuss the fluctuations around coherent states for comparison. Then we derive the Bogoliubov Hamiltonian and study its dynamics in Section \ref{sec:Bogoliubov}. The main result is proved in Section \ref{sec:reformulation}.

\bigskip

\noindent\textbf{Acknowledgement.} We acknowledge financial support from the European Research Council under the European Community's Seventh Framework Programme FP7/2007-2013: Grant Agreement MNIQS 258023 (M.L.~and P.T.N.) and Grant Agreement MAQD 240518 (B.S.).

\section{Main result} \label{sec:main-result}

\noindent{\bf Notation.} Any Hilbert space we consider has an inner product which is conjugate linear in the first variable and linear in the second. We always denote by $C$ a (large) positive constant which depends only on $w$.
The symmetric tensor product $\Psi_k\otimes_s\Psi_\ell\in\mathfrak{H}^{k+\ell}$ of two symmetric functions $\Psi_k\in \mathfrak{H}^k$ and $\Psi_\ell\in\mathfrak{H}^\ell$ is defined by
\begin{multline*}
\Psi_k\otimes_s\Psi_\ell(x_1,...,x_{k+\ell})\\
=\frac{1}{\sqrt{k!\ell!(k+\ell)!}}\sum_{\sigma\in\mathfrak{S}_{k+\ell}}\Psi_k(x_{\sigma(1)},...,x_{\sigma(k)})\Psi_\ell(x_{\sigma(k+1)},...,x_{\sigma(k+\ell)}).
\end{multline*}

\begin{theorem}[Dynamics of the fluctuations around a Hartree state]\label{thm:fluctuations}
Assume that $w^2 \le C(1-\Delta)$ on $ \mathfrak{H}=L^2(\mathbb{R}^d)$. Let $u_0\in H^1(\mathbb{R}^d)$ with $\|u_0\|_{L^2}=1$ and let $u(t)$ be the unique solution to the Hartree equation~\eqref{eq:Hartree-equation}. Let $\mathfrak{H}_+(t):=\{u(t)\}^\perp \subset \mathfrak{H}$, and $\mathfrak{H}_+^n(t)=\bigotimes_{\rm sym}^n\mathfrak{H}_+(t)$ be the $n$-fold symmetric tensor product.

Consider a sequence of functions $\phi_{n,0}\in \mathfrak{H}_+^n(0)$ such that, for any $\varepsilon > 0$,
\begin{equation}
\sum_{n=0}^\ii \norm{\phi_{n,0}}_{\mathfrak{H}^n}^2=1 \quad \text{ and} \quad \sum_{n=0}^\ii \, n \, \langle \phi_{n,0} , (1-\Delta_{x_1}) \phi_{n,0} \rangle < \ii.
\label{eq:condition_initial_datum}
\end{equation}
If $\Psi_N(t)$ is the solution of the Schr\"odinger equation~\eqref{eq:Schrodinger-dynamics} with initial datum
\begin{equation}
\Psi_N(0)=\sum_{n=0}^N u_0^{\otimes (N-n)}\otimes_s\phi_{n,0},
\label{eq:form_initial_datum}
\end{equation}
then we have for all times $t\ge 0$,
\begin{equation}
\lim_{N\to\ii}\norm{\Psi_N(t)-\sum_{n=0}^N u(t)^{\otimes (N-n)}\otimes_s\phi_n(t)}_{\mathfrak{H}^N}=0,
\label{eq:fluctuations}
\end{equation}
with $t\mapsto \phi_n(t)\in\mathfrak{H}_+^n(t)$ a continuous map for all $n\geq0$. More precisely, $\Phi(t):=(\phi_n(t))_{n\geq0}$ is the unique solution to the effective Bogoliubov equation in Fock space $\mathcal{F}=\bigoplus_{n\geq 0}\mathfrak{H}^n:$
\begin{eqnarray} \label{eq:Bogoliubov-dynamics}
\left
\{ \begin{gathered}
  i\, \dot{\Phi}(t) = \mathbb{H}(t)  \Phi(t) \hfill \\
  \Phi(0)=\Phi_0:=(\phi_{n,0})_{n\geq0}, \hfill \\
\end{gathered}
\right.
\end{eqnarray}
where $\mathbb{H}(t)$ is a quadratic Hamiltonian in Fock space, written in second quantized form
\begin{align} \label{eq:Bogoliubov-Hamiltonian}
\mathbb{H}(t)=&\int_{\mathbb{R}^d}a^\dagger(x) \big(-\Delta+|u(t)|^2\ast w-\mu(t) + K_1(t) \big)a(x)\,dx\\
&+\frac12\iint_{\mathbb{R}^d\times\mathbb{R}^d}\Big(K_2(t,x,y)a^\dagger(x)a^\dagger(y)+\overline{K_2(t,x,y)}a(x)a(y)\Big)dx\,dy.\nonumber
\end{align}
Here $K_1(t)$ is the operator $K_1(t)={Q(t)}\tilde K_1(t)Q(t)$ which acts on $L^2(\mathbb{R}^d)$, where $Q(t)=1-|u(t)\rangle\langle u(t)|$ and where $\tilde{K}_1$ is the Hilbert-Schmidt operator on $L^2(\mathbb{R}^d)$ with kernel $\tilde{K}_1(t,x,y)=u(t,x)w(x-y)\overline{u(t,y)}$. On the other hand, $K_2(t)=(Q(t)\otimes Q(t)) \tilde K_2(t)$ is the projection onto $\mathfrak{H}_+(t)\otimes\mathfrak{H}_+(t)$ of the symmetric two-body function $\tilde K_2(t,x,y)=u(t,x)w(x-y)u(t,y)$.
\end{theorem}

Here $a^\dagger(x)$ and $a(x)$ are the usual creation and annihilation operators (see Section \ref{sec:coherent} for definitions). The Bogoliubov Hamiltonian $\mathbb{H}(t)$ and the effective equation \eqref{eq:Bogoliubov-dynamics} will be discussed in detail in Section~\ref{sec:Bogoliubov} below. It is not obvious from the definition of $\mathbb{H}(t)$ that
$$\Phi(t)\in \mathcal{F}_+(t):=\bigotimes_{n\ge 0} \mathfrak{H}_+^n(t) \subset \mathcal{F}$$
for all times, but this fact will be proved in the sequel. Note that the Bogoliubov Hamiltonian $\mathbb{H}(t)$ does {\em not} preserve the subspace $\mathcal{F}_+(t)$ due to the term involving the Hartree mean-field Hamiltonian
$$h(t):=-\Delta+|u(t)|^2\ast w-\mu(t)$$
which itself does not preserve $\mathfrak{H}_+(t)$. In fact, for all times $t$, the vector $\mathbb{H}(t) \Phi(t)$ belongs to the space tangent to the manifold $\{ \tilde \Phi \in \mathcal{F}_+(t),\ \|\tilde \Phi\|_\mathcal{F}=1\}$ at the point $\Phi(t)$, which is the correct condition to have $\Phi(t)\in\mathcal{F}_+(t)$.

It is possible to rewrite~\eqref{eq:Bogoliubov-dynamics} as an infinite system of coupled linear equations:
\begin{equation*}
 \left\{\begin{array}{rcl}
i\dot \phi_0(t)&=& \displaystyle\sqrt2\iint_{\mathbb{R}^d\times\mathbb{R}^d}\overline{K_2(t,x_1,x_2)}\phi_2(t,x_1,x_2)\,dx_1\,dx_2\\[0.3cm]
i\dot \phi_1(t,x_1)&=& (h(t)+K_1(t))\phi_1(t,x_1)\\[0.3cm]
&&~+\displaystyle\sqrt6\iint_{\mathbb{R}^d\times\mathbb{R}^d}\overline{K_2(t,x_2,x_3)}\phi_3(t,x_1,x_2,x_3)dx_2\,dx_3\\[0.3cm]
i\dot \phi_2(t,x_1,x_2)&=& \big((h(t)+K_1(t))_{x_1}+(h(t)+K_1(t))_{x_2}\big)\phi_2(t,x_1,x_2)\\[0.3cm]
&&~+ \displaystyle \sqrt2\, K_2(t,x_1,x_2)\phi_0(t)\\[0.3cm]
&&~+\displaystyle\sqrt{12}\iint_{\mathbb{R}^d\times\mathbb{R}^d}\overline{K_2(t,x_3,x_4)}\phi_4(t,x_1,x_2,x_3,x_4)dx_3\,dx_4\\[0.3cm]
&\vdots&
\end{array}
\right.
\end{equation*}
In words, each $\phi_n(t)$ evolves with the particle-preserving Hamilton operator $\sum_{j=1}^n (h(t)+K_1(t))_{x_j}$ and it is coupled to its second left and second right neighbors $\phi_{n-2}(t)$ and $\phi_{n+2}(t)$, through the function $K_2(t,x,y)$. Note that the even functions $\phi_{2j}(t)$ and odd functions $\phi_{2j+1}(t)$ are not coupled.

The reader should think of the simplest case where $\phi_{0,0}=1$ and $\phi_{n,0}=0$ for all $n\geq1$, which means that the initial datum is the Hartree state
$$\Psi_N(0)=u_0^{\otimes N}$$
for all $N$. In this case the condition~\eqref{eq:condition_initial_datum} boils down to $\|u_0\|_{L^2}=1$ and $u_0\in H^1(\mathbb{R}^d)$. Then~\eqref{eq:fluctuations} gives the exact behavior of the wave function $\Psi_N(t)$ for large $N$, in terms of the ($N$-independent) functions $\phi_n(t)$ which describe the fluctuations around it. Note that in this simple case, it can be verified that the odd functions $\phi_{2j+1}(t)\equiv0$ for all $j$ but the even functions $\phi_{2j}(t)$ are all different from 0 for all $j$ and $t>0$ small enough, provided that $K_2(0)\neq0$. For instance the derivative of $\phi_2$ does not vanish at time zero,
$$\dot\phi_2(0)=- \frac{i}2 K_2(0)\neq 0,$$
and therefore $\phi_2(t)\neq0$ for small times. We deduce from~\eqref{eq:fluctuations} that $\Psi_N(t)$ does \emph{not} stay close to the pure Hartree state $u(t)^{\otimes N}$ in norm in the limit $N\to\ii$. In fact, in the expansion in~\eqref{eq:fluctuations} there are infinitely many terms in the series which do not vanish if $K_2(0)\neq0$, and this is always the case except when $w\equiv 0$.

Finally, we remark that the estimate~\eqref{eq:fluctuations} is much stronger than the convergence of the density matrices. This is because any of the terms $u(t)^{\otimes (N-n)}\otimes_s\phi_n(t)$ in the series has most of its particles condensed into the Hartree state $u(t)$ and only $n$ excited particles outside of the condensate. The functions $\phi_n(t)$ do not participate to the first order of the density matrices in the limit $N\to\ii$. We quickly explain this now for clarity.

\begin{corollary}[Convergence of density matrices]\label{cor:Hartree}
Under the assumptions of Theorem~\ref{thm:fluctuations}, we have
\begin{equation}\label{eq:tr-norm} \lim_{N\to\ii}\tr\left|\gamma^{(k)}_{\Psi_N(t)}-|u(t)\rangle\langle u(t)|^{\otimes k} \right|=0 \end{equation}
for every fixed $k\geq1$ and $t\in\mathbb{R}^+$.
\end{corollary}

\begin{proof
Recall that for any normalized $\Psi,\Psi'\in\mathfrak{H}^N$
\begin{equation}
\tr\left|\gamma^{(1)}_{\Psi}-\gamma^{(1)}_{\Psi'}\right|\leq 2\norm{\Psi-\Psi'},
\label{eq:estimate_DM}
\end{equation}
which follows from the estimate
\begin{align*}
\left|\tr\big(A(\gamma^{(1)}_{\Psi}-\gamma^{(1)}_{\Psi'})\big)\right|&=\left|\pscal{\Psi,A\otimes 1_{\mathfrak{H}^{N-1}}\Psi}-\pscal{\Psi',A\otimes 1_{\mathfrak{H}^{N-1}}\Psi'}\right|\\
&\leq 2\norm{A}\,\norm{\Psi-\Psi'}
\end{align*}
for any bounded operator $A$. Combining~\eqref{eq:fluctuations} and~\eqref{eq:estimate_DM} we deduce that
\begin{eqnarray} \label{eq:density-cv-1}
\lim_{N\to\ii}\tr\left|\gamma^{(1)}_{\Psi_N(t)}-\gamma^{(1)}_{\sum_{n=0}^N u(t)^{\otimes (N-n)}\otimes_s\phi_n(t)}\right|=0.
\end{eqnarray}
Now we compute, with $P(t)=|u(t)\rangle\langle u(t)|$,
\begin{align*}
&\pscal{u(t),\gamma^{(1)}_{\sum_{n=0}^N u(t)^{\otimes (N-n)}\otimes_s\phi_n(t)} u(t)}\\
&\qquad=\sum_{n,m=0}^N\pscal{u(t)^{\otimes (N-n)}\otimes_s\phi_n(t),\left(P(t)\otimes 1_{\mathfrak{H}^{N-1}}\right)u(t)^{\otimes (N-m)}\otimes_s\phi_m(t)}\\
&\qquad=\sum_{n=0}^N\left(1-\frac{n}{N}\right)\norm{\phi_n(t)}^2_{\mathfrak{H}^n}
\end{align*}
where we have used the fact that $\phi_n(t)\in\mathfrak{H}_+^N(t)$ for all times. The last term tends to 1 by the dominated convergence theorem for series (actually it will be proved later that $\sum n\|\phi_n(t)\|^2$ stays finite for all times, which even gives the convergence rate $C(t)N^{-1}$).
It follows that
\begin{eqnarray} \label{eq:density-cv-2}
 \lim_{N\to \infty} \tr\, \left| \gamma^{(1)}_{\sum_{n=0}^N u(t)^{\otimes (N-n)}\otimes_s\phi_n(t)} - |u(t)\rangle\langle u(t)| \right| = 0. \end{eqnarray}
From (\ref{eq:density-cv-1}) and (\ref{eq:density-cv-2}) we obtain the convergence of the one-particle density matrices. The result for higher density matrices follows from a similar argument (or follows directly from the convergence of the one-particle density matrices ~\cite[Corollary 2.4]{LewNamRou-13}).
\end{proof}

\begin{remark}[Convergence rates]\label{rmk:rate}
In this work, we do not focus on convergence rates. However, using the techniques that we will discuss in the next sections, one can also obtain a quantitative bound for the difference in (\ref{eq:fluctuations}). The optimal bound is of the order $N^{-1/2}$, for every fixed $t \in \mathbb{R}$, and can be proved proceeding similarly as we do below, with the difference that instead of bounding only the expectation of $\mathcal{N}$, in Lemma \ref{le:N+-bound} we would need an estimate for the expectation of higher power of $\mathcal{N}$ (correspondingly, one needs a stronger assumption on the expectation of powers of $\mathcal{N}$ in the initial data). Such estimates can be established similarly as in the proof of Lemma \ref{le:N+-bound}. Although a bound of the order $N^{-1/2}$ for the fluctuation in (\ref{eq:fluctuations}) immediately implies a bound of the same order for the reduced densities in (\ref{eq:tr-norm}), the optimal estimate for the latter is proportional to $N^{-1}$ (see \cite{CheLeeSch-11}).
\end{remark}

\begin{remark}[Extensions] \label{rm:ext} We can consider not only the typical Hamiltonian $H_N$ in (\ref{eq:def_H_N}) but also many other cases.
\begin{itemize}

\item[(i)] (External potential) We can add to the Laplacian $-\Delta$ in $H_N$ an external potential $V(x)$ which satisfies
$$V \in L_{\rm loc}^{d/2} (\mathbb{R}^d,\mathbb{R}) \quad {\text{and}}~~-V(x) \le \alpha (-\Delta)+C$$
for some $0\le \alpha<1$. Note that $V(x)$ is allowed to grow at infinity, which corresponds to a confined system. The potential $V(x)$ appears in the Bogoliubov Hamiltonian in a natural way. Our results in Theorem \ref{thm:fluctuations} still hold true provided that $w^2 \le C(1-\Delta)$ and that the second condition in (\ref{eq:condition_initial_datum}) on the initial datum is replaced by
$$
\sum_{n=0}^\ii \, n \, \langle \phi_{n,0} , (1-\Delta_{x_1}+|V(x_1)|) \phi_{n,0} \rangle < \ii.
$$

\item[(ii)] (Boson stars) We can also consider bosons living in $\mathbb{R}^3$ with the one-body pseudo-relativistic kinetic operator $\sqrt{1-\Delta}$ and the Newton interaction $w(x-y)=-\kappa |x-y|^{-1}$, $0<\kappa<2/\pi$. If $u_0\in H^1(\mathbb{R}^3)$, then the Hartree dynamics $u(t)$ is well-posed in $H^1(\mathbb{R}^3)$ (see \cite{Lenzmann-07}) and our results in Theorem \ref{thm:fluctuations} still hold true under the same $H^1$-condition (\ref{eq:condition_initial_datum}) on the initial datum.
\end{itemize}
\end{remark}

\section{Comparison with coherent states}\label{sec:coherent}

One common approach to derive Hartree's nonlinear equation from the many-body dynamics is to work in Fock space and to use coherent states. We quickly recall this method here and we compare it with our main result. The Fock space is
$$\mathcal{F}=\mathbb{C}\oplus\bigoplus_{n\geq1}\mathfrak{H}^n$$
and it is useful to introduce the creation and annihilation operators $a^\dagger(f)$ and $a(g)$ which are unbounded operators acting on $\mathcal{F}$. More precisely, $a^\dagger(f)$ is defined by
$$a^\dagger(f)\left(\phi_0\oplus\bigoplus_{n\geq1} \phi_n\right)=0\oplus \phi_0 f \oplus\bigoplus_{n\geq1} \big(f\otimes_s \phi_n\big)$$
and $a(f)$ is its adjoint. The operator-valued distributions $a^\dagger(x)$ and $a(x)$ used in (\ref{eq:Bogoliubov-Hamiltonian}) may be defined so that
$$ a^\dagger (f)=\int _{\mathbb{R}^d} f(x)a^\dagger(x) dx \quad {\rm and}\quad a(f)=\int _{\mathbb{R}^d} \overline{f(x)}a^\dagger(x) dx$$
for all $f\in \mathfrak{H}$. These operators satisfy the canonical commutation relations
\begin{align*} [a(f),g(g)] &=[a^\dagger(f),a^\dagger(g)]=0,\quad [a(f),a^\dagger(g)]=\langle f,g\rangle_{\mathfrak{H}} \\
[a(x),a(y)]&=[a^\dagger(x),a^\dagger(y)]=0,\quad [a(x),a^\dagger(y)]=\delta(x-y)
\end{align*}
for all $f,g\in \mathfrak{H}$ and for all $x,y\in \mathbb{R}^d$.

The \emph{Weyl unitary operator} is defined by
$$W(f):=\exp\big(a^\dagger(f)-a(f)\big)$$
and a \emph{coherent state} is a Weyl-rotation of the vacuum $\Omega=1\oplus0\cdots$:
$$W(f)\Omega=e^{-\norm{f}^2/2}\sum_{n\geq0}\frac{1}{\sqrt{n!}}f^{\otimes n}.$$
The average particle number of a coherent state is, therefore,
$$e^{-\norm{f}^2}\sum_{n\geq0}\frac{n}{n!}\norm{f}^{2 n}=\norm{f}^2.$$

Instead of starting the many-particle dynamics in the neighborhood of a Hartree state in $\mathfrak{H}^N$, one can also use a coherent state in Fock space. The analysis is simplified by the useful algebraic properties of coherent states. For example, coherent states are eigenvectors of annihilation operators, which follows from the action of the Weyl rotation $W(f)$ on the creation and annihilation operators
\begin{eqnarray} \label{eq:Weyl-action}
W(f)^* a^\dagger(g) W(f)=a^\dagger(g) + \langle f,g \rangle, \quad W(f)^* a(g) W(f)=a(g) + \langle g,f \rangle.
\end{eqnarray}

For a precise discussion on the dynamics in Fock space, let us introduce the second quantized forms of operators. If $A$ is a one-body operator on $\mathfrak{H}$ with kernel $A(x,y)$, then its second quantized form $d\Gamma(A)$ is an operator on $\mathcal{F}$ defined by
$$d\Gamma(A):= \bigoplus_{n\ge 1} \sum_{j=1}^n A_j= \iint_{\mathbb{R}^d\times \mathbb{R}^d} a^\dagger(x) A(x,y) a(y) \,dx\,dy,$$
which is frequently written as $\int a^\dagger(x) A a(x) \,dx$ for short. For example, $\mathcal{N}:=d\Gamma(\1_\mathfrak{H})$ is the number operator. Similarly, if $W$ is a symmetric operator on $\mathfrak{H}^2$ with kernel $W(x,y;x',y')$, then its second quantized form is given by
\begin{multline*}
\bigoplus_{n\ge 2} \bigg( \sum_{1\le j<k \le n} W_{j,k} \bigg)\\= \frac{1}{2}\iiiint W(x,y;x',y')a^\dagger(x)a^\dagger(y)a(x') a(y') \,dx\,dy\,dx'\,dy'.
\end{multline*}
In particular, the formal kernel of the multiplication operator $w(x-y)$ on $\mathfrak{H}^2$ is $w(x-y)\delta_{(x,y)}(x',y')$ and the second quantized form of $H_N$ is an operator on Fock space:
\begin{equation}
\mathbb{H}_N=d\Gamma(-\Delta) +\frac{1}{2(N-1)}\iint_{\mathbb{R}^d\times\mathbb{R}^d} w(x-y)a^\dagger (x)a^\dagger (y) a(x) a(y)\nonumber\\
\end{equation}
which coincides with $H_N$ in the sector $\mathfrak{H}^N$. Because of the constant in front of the interaction which depends on $N$ and not on the number operator $\mathcal{N}$, the restriction of $\mathbb{H}_N$ to another $k$-particle subspaces is not related to $H_k$.
Under certain assumptions on $w$, for example (\ref{eq:asp-w}), the dynamics generated by $\mathbb{H}_N$ is well-defined on the quadratic form domain of $d\Gamma(1-\Delta)$.

It was first proved by Hepp \cite{Hepp-74} for regular interactions, and then extended by Ginibre-Velo \cite{GinVel-79,GinVel-79b} for singular interactions (see also~\cite{RodSch-09,GriMacMar-10,GriMacMar-11}) that if the initial state (in Fock space) is close to a coherent state $W(\sqrt{N} u_0)\Omega$, then the dynamics generated by $\mathbb{H}_N$ stays close to $W(\sqrt{N} u(t))\Omega$ for all $t\ge 0$ in the sense of density matrices and it is also possible to identify the fluctuations. In our typical setting, this  result can be stated as follows.

\begin{theorem}[Dynamics of the fluctuations around a coherent state]\label{thm:fluctuations_coherent}
Assume that $w^2 \le C(1-\Delta)$. Let $u_0\in H^1(\mathbb{R}^d)$ with $\|u_0\|_{L^2}=1$ and let $u(t)$ be the unique solution to the nonlinear Hartree equation~\eqref{eq:Hartree-equation}.

Consider an initial datum $\Psi_{N,0}$ in Fock space which is such that
$$W(\sqrt{N} u_0)^*\Psi_{N,0} \to \Xi_0$$
strongly in $\mathcal{F}$ and weakly in $Q({\ensuremath{\rm d}\Gamma}(1-\Delta))$. Let $\Psi_N(t)=\exp(-it\mathbb{H}_N)\Psi_{N,0}$ be the solution of the Schr\"odinger equation in $\mathcal{F}$, with initial datum $\Psi_{N,0}$. Then we have for all times $t\ge 0$:
\begin{equation}
\lim_{N\to\ii}W(\sqrt{N} u(t))^*\Psi_N(t) = \Xi(t),
\label{eq:fluctuations_coherent}
\end{equation}
strongly in $\mathcal{F}$ and weakly in $Q({\ensuremath{\rm d}\Gamma}(1-\Delta))$, where $\Xi(t)$ is the unique solution to the effective equation in Fock spac
\begin{eqnarray} \label{eq:Bogoliubov-dynamics_coherent}
\left
\{ \begin{gathered}
  i\, \dot{\Xi}(t) = \tilde\mathbb{H}(t)  \Xi(t) \hfill \\
  \Xi(0)=\Xi_0. \hfill \\
\end{gathered}
\right.
\end{eqnarray}
Here $\tilde \mathbb{H}(t)$ is a quadratic Hamiltonian on $\mathcal{F}$,
 written in second quantized form
\begin{align*}
\tilde\mathbb{H}(t)=&\int_{\mathbb{R}^d}a^\dagger(x) \big( -\Delta +|u(t)|^2\ast w-\mu(t) + \tilde K_1(t) \big)a(x)\,dx\nonumber\\
&+\frac12\iint_{\mathbb{R}^d\times\mathbb{R}^d}\Big(\tilde K_2(t,x,y)a^\dagger(x)a^\dagger(y)+\overline{\tilde K_2(t,x,y)}a(x)a(y)\Big)dx\,dy,
\end{align*}
with $\tilde K_1(t,x,y)=u(t,x)w(x-y)\overline{u(t,y)}$, $\tilde K_2(t,x,y)=u(t,x)w(x-y)u(t,y)$.
\end{theorem}

Theorem \ref{thm:fluctuations_coherent} can be proved using the argument of the proof of Theorem \ref{thm:fluctuations} and we will omit the details. Although the result in the case of coherent states looks very similar to the Hartree case, the effective Bogoliubov Hamiltonians $\tilde\mathbb{H}(t)$ and $\mathbb{H}(t)$ are \emph{not} the same. The part involving $h(t)$ is identical, but the other terms involve the functions $\tilde K_j(t,x,y)$'s instead of the projected ones $K_j(t,x,y)$'s.

In fact, the coherent state approach can also be used to study fluctuations around a Hartree state, if the initial data is exactly a Hartree state. This can be seen as follows. Theorem~\ref{thm:fluctuations_coherent} implies that
$$\lim_{N\to\ii}\norm{e^{it\mathbb{H}_N}W(\sqrt{N}u_0)\Omega-W(\sqrt{N}u(t))\tilde{\mathbb{U}}(t,0)\Omega}_{\mathcal{F}}=0,$$
where $\tilde{\mathbb{U}}(t,0)$ is the unitary propagator associated with the time-dependent Bogoliubov Hamiltonian $\tilde\mathbb{H}(t)$ in Fock space. A more precise analysis shows that
\[ \norm{e^{it\mathbb{H}_N}W(\sqrt{N}u_0)\Omega-W(\sqrt{N}u(t))\tilde{\mathbb{U}}(t,0)\Omega}_{\mathcal{F}} \lesssim N^{-1/2} \]
for any fixed time $t$. Hence projecting onto the $N$-particle subspace, we find
\begin{equation}
\norm{\frac{(N/e)^{N/2}}{\sqrt{N!}}e^{itH_N}u_0^{\otimes N}-\1_{\mathfrak{H}^N}W(\sqrt{N}u(t))\tilde{\mathbb{U}}(t,0)\Omega}_{\mathfrak{H}^N} \lesssim N^{-1/2} \, .
\label{eq:projection_coherent}
\end{equation}
Multiplying with $\sqrt{N!} / (N/e)^{N/2} \simeq N^{1/4}$, we get
$$ \lim_{N \to \infty} \norm{e^{itH_N}u_0^{\otimes N}-(2\pi N)^{1/4}\1_{\mathfrak{H}^N}W(\sqrt{N}u(t))\tilde{\mathbb{U}}(t,0)\Omega}_{\mathfrak{H}^N}= 0.$$
The effective dynamics $(2\pi N)^{1/4}\1_{\mathfrak{H}^N}W(\sqrt{N}u(t))\tilde{\mathbb{U}}(t,0)\Omega$ agrees with our fluctuations in~\eqref{eq:fluctuations} to the leading order, even if it looks much more complicated. In fact, this analysis can be extended to initial $N$-particle wave functions
having the form $d_N \, \1_{\mathfrak{H}^N} W(\sqrt{N} u_0) \Psi$, assuming the normalization constant $d_N$ to be small enough (of order $N^{1/4}$) and $\Psi$ to have sufficiently small number of particles (of order one, independent of $N$). In other words, to get information about the evolution of an $N$-particle initial data using the coherent state method requires further assumptions on the initial datum, which guarantee sufficient closeness to an Hartree state (in an appropriate sense). On the other hand, our direct method in $\mathfrak{H}^N$ applies to any initial datum of the general form~\eqref{eq:form_initial_datum} satisfying~\eqref{eq:condition_initial_datum}, for which no convergence rate is known.

\section{The Bogoliubov dynamics}\label{sec:Bogoliubov}

\subsection{Describing fluctuations around a given $u(t)$} \label{subsec:UN}
Consider a time-dependent normalized vector $u(t)$ in $\mathfrak{H}$ (which does not necessarily satisfy the Hartree equation). As in~\cite[Sec. 2.3]{LewNamSerSol-13}, we can write any function $\Psi\in \mathfrak{H}^N$ as follows
\begin{equation}
\Psi=\psi_0\, u(t)^{\otimes N}+ \psi_1 \otimes_s u(t)^{\otimes (N-1)} +\psi_2 \otimes_s u(t)^{\otimes (N-2)}+\cdots + \psi_N
\label{eq:decomp_Psi}
\end{equation}
where $\psi_n\in \mathfrak{H}_+(t)^{n}:= \bigotimes_{\rm sym}^n \mathfrak{H}_+(t)$ with $\mathfrak{H}_+(t):= \{u(t)\}^{\bot}$. Following~\cite{LewNamSerSol-13}, we define the unitary
\begin{equation}
 \label{eq:def-unitary-UN}
\begin{array}{cccl}
U_{N}(t): & \mathfrak{H}^N & \to & \displaystyle\mathcal{F}_+(t)^{\leq N}=\bigoplus_{n=0}^N \mathfrak{H}_+(t)^n \\[0.3cm]
 & \Psi & \mapsto & \psi_0\oplus \psi_1 \oplus\cdots \oplus \psi_N.
\end{array}
\end{equation}
Note that we have the inclusions
\begin{equation}
\mathcal{F}_+(t)^{\leq N}\subset \mathcal{F}_+(t) = \bigoplus_{n=0}^\infty \mathfrak{H}_+(t)^n \subset \mathcal{F}.
\label{eq:inclusions}
\end{equation}
Therefore, we can always see $U_{N}(t)$ as a partial isometry from $\mathfrak{H}^N$ to $\mathcal{F}$, where $U_{N}(t)^*$ is extended by $0$ outside of $\mathcal{F}_+(t)^{\leq N}$.

As shown in ~\cite[Proposition 14]{LewNamSerSol-13}, the operators $U_{N}(t)$ and $U_{N}(t)^*$ can be equivalently written as
\begin{eqnarray} \label{eq:alt-def-UN}
U_{N}(t)= \bigoplus_{j=0}^N Q(t)^{\otimes j} \frac{a (u(t))^{N-j}}{\sqrt{(N-j)!}} ,\quad U_{N}(t)^* = \bigoplus _{j = 0}^N\frac{{{{a^\dagger (u(t))}^{N - j}}}}{{\sqrt {(N - j)!} }},
\end{eqnarray}
where $Q(t)=1-|u(t)\rangle\langle u(t)|$, and $a(u(t))$ and $a^\dagger(u(t))$ are the usual annihilation and creation operators on $\mathcal{F}$. Moreover, we have the following identities on $\mathcal{F}_+(t)^{\le N}$:
\begin{eqnarray*}
U_{N}(t) \, a^\dagger (u(t)) a(u(t)) \,U_{N}(t)^* &=& N- \mathcal{N}_+(t) ,\hfill\\
U_{N}(t)\, a^\dagger(f) a(u(t)) \,U_{N}(t)^* &=& a^\dagger(f) \sqrt{N-\mathcal{N}_+(t)},\hfill\\
U_{N}(t) \,a^\dagger (u(t)) a(f) \,U_{N}(t)^* &=& \sqrt{N-\mathcal{N}_+(t)}\, a(f),\hfill\\
U_{N}(t)\, a^\dagger(f) a(g) \,U_{N}(t)^* &=& a^\dagger(f) a(g),
\end{eqnarray*}
for all $f,g\in \mathfrak{H}_+$, where $\mathcal{N}_+(t)=\mathcal{N}-a^\dagger(u(t)) a(u(t))$ is the number operator on $\mathcal{F}_+(t)$.

In the analysis of the dynamics of $U_{N}(t) \Psi_{N}(t)$, not only $U_{N}(t)$ but also its time-derivative plays an important role. We have
\begin{lemma} Let $u(t)$ be an arbitrary (sufficiently regular) trajectory on $\mathfrak{H}$ satisfying $\|u(t)\|=\|u(0)\|$ for all $t\ge 0$. Then the time-derivative of $U_N(t)$ is
\begin{multline} \label{eq:derivative-UN}
i\dot{U}_{N}(t)=\Big( a^\dagger(u(t))a(Q(t)i \dot u(t))-\sqrt{N-\mathcal{N}_+(t)}\,a(Q(t)i\dot u(t))\\
-a^\dagger(Q(t)i\dot u(t))\sqrt{N-\mathcal{N}_+(t)}-\pscal{i\dot u(t),u(t)}(N-\mathcal{N}_+(t))\Big)U_{N}(t).
\end{multline}
\end{lemma}
Note that the last three terms in the right side of (\ref{eq:derivative-UN})
do not go outside of $\mathcal{F}_+(t)$, but the first does.

\begin{proof} We start by writing
\begin{equation}
U_{N}(t)= \bigoplus_{k=0}^N \; Q(t)^{\otimes k}\frac{a(u(t))^{N-k}}{\sqrt{(N-k)!}}=:\bigoplus_{k=0}^N \mathcal{L}_k(t).
\label{eq:formula-psi-k-a}
\end{equation}
Now taking the time-derivative, we get
\begin{align} \label{eq:derivative-Lk}
\dot\mathcal{L}_k(t) &=\frac{N-k}{\sqrt{(N-k)!}}Q(t)^{\otimes k} a(\dot{u}(t))a(u(t))^{N-k-1} \nonumber\\
&~+ \frac{1}{\sqrt{(N-k)!}}\;\Big(\dot{Q}(t)\otimes Q(t)\otimes\cdots\otimes Q(t)+\cdots\nonumber\\
& \qquad\qquad\qquad\qquad\quad + Q(t)\otimes\cdots\otimes Q(t)\otimes \dot{Q}(t)\Big)a(u(t))^{N-k}.
\end{align}
In the first term of the right side of (\ref{eq:derivative-Lk}), we have used that $a(\dot{u}(t))$ and $a(u(t))$ commute. It is clear that for any $v$, we have
$$Q(t)^{\otimes (n+1)} a^\dagger(v)_{|\mathfrak{H}^{n}}={a^\dagger(Q(t)v)Q(t)^{\otimes n}}_{|\mathfrak{H}^{n}}$$
and
$$Q(t)^{\otimes n} a(v)_{|\mathfrak{H}^{n+1}}={a(Q(t)v)Q(t)^{\otimes (n+1)}}_{|\mathfrak{H}^{n+1}}+\pscal{v,u(t)}Q(t)^{\otimes n} a(u(t))_{|\mathfrak{H}^{n+1}}.$$
So the first term of the right side of (\ref{eq:derivative-Lk}) can be rewritten as
\begin{align} \label{eq:derivative-Lk-a}
&\frac{N-k}{\sqrt{(N-k)!}}Q(t)^{\otimes k}\Big( a(\dot{u}(t))a(u(t))^{N-k-1}\Psi\Big)\nonumber\\
&=\frac{\sqrt{N-k}}{\sqrt{(N-k-1)!}} a(Q(t)\dot{u}(t)) Q(t)^{\otimes (k+1)}a(u(t))^{N-k-1}\nonumber\\
&\qquad+\pscal{\dot{u}(t),u(t)}\frac{N-k}{\sqrt{(N-k)!}}Q(t)^{\otimes k} a(u(t))^{N-k}\nonumber\\
&=\sqrt{N-k}\, a(Q(t)\dot{u}(t)) \mathcal{L}_{k+1}(t)+\pscal{\dot{u}(t),u(t)}(N-k)\mathcal{L}_k(t).
\end{align}
Note that
\begin{align*}
\dot{Q}(t)=-|u(t)\rangle\langle \dot{u}(t)|-|\dot{u}(t)\rangle\langle {u}(t)|=-|u(t)\rangle\langle \dot{u}(t)|Q(t)-Q(t)|\dot u(t)\rangle\langle {u}(t)|
\end{align*}
because
\begin{align*}
|u(t)\rangle\langle \dot{u}(t)|&=|u(t)\rangle\langle \dot{u}(t)|Q(t)+\pscal{\dot{u}(t),u(t)}P(t)\\
|\dot u(t)\rangle\langle {u}(t)|&=Q(t)|\dot u(t)\rangle\langle {u}(t)|+\pscal{u(t),\dot u(t)}P(t),
\end{align*}
where $P(t)=|u(t)\rangle \langle u(t)|$, and
$$\pscal{\dot{u}(t),u(t)}+\pscal{u(t),\dot u(t)}=\frac{d}{dt} ||u(t)||^2=0.$$
Since $a^\dagger(f)a(g)={\ensuremath{\rm d}\Gamma}(|f\rangle\langle g|)$, we find
\begin{align*}
&\dot{Q}(t)\otimes Q(t)\otimes\cdots\otimes Q(t)+\cdots+ Q(t)\otimes\cdots\otimes Q(t)\otimes \dot{Q}(t)\\
&\qquad=-\left(\sum_{j=1}^k|u(t)\rangle\langle Q(t)\dot{u}(t)|_j\right)Q(t)^{\otimes k}-Q(t)^{\otimes k}\left(\sum_{j=1}^k|Q(t)\dot u(t)\rangle\langle {u}(t)|_j\right)\\
&\qquad=-a^\dagger(u(t))a(Q(t)\dot u(t))Q(t)^{\otimes k}-Q(t)^{\otimes k}a^\dagger(Q(t)\dot u(t))a(u(t))\\
&\qquad=-a^\dagger(u(t))a(Q(t)\dot u(t))Q(t)^{\otimes k}-a^\dagger(Q(t)\dot u(t))Q(t)^{\otimes (k-1)}a(u(t)).
\end{align*}
So the second term of the right side of (\ref{eq:derivative-Lk}) can be rewritten as
\begin{flalign} \label{eq:derivative-Lk-b}
&\frac{1}{\sqrt{(N-k)!}}\;\Big(\dot{Q}(t)\otimes Q(t)\otimes\cdots\otimes Q(t)+\cdots \nonumber\\
&\qquad\qquad\qquad\qquad\qquad+ Q(t)\otimes\cdots\otimes Q(t)\otimes \dot{Q}(t)\Big)a(u(t))^{N-k}\nonumber\\
&=-a^\dagger(u(t))a(Q(t)\dot u(t))\mathcal{L}_k-\sqrt{N-k+1}\,a^\dagger(Q(t)\dot u (t))\mathcal{L}_{k-1}.
\end{flalign}

Substituting (\ref{eq:derivative-Lk-a}) and (\ref{eq:derivative-Lk-b}) into (\ref{eq:derivative-Lk}) and taking the sum over $k$, we deduce from (\ref{eq:formula-psi-k-a}) the formula
\begin{multline*}
\dot{U}_N(t)=\Big(-a^\dagger(u(t))a(Q(t)\dot u(t))+\sqrt{N-\mathcal{N}_+(t)}\,a(Q(t)\dot u(t))\\
-a^\dagger(Q(t)\dot u(t))\sqrt{N-\mathcal{N}_+(t)}+\pscal{\dot u(t),u(t)}(N-\mathcal{N}_+(t))\Big)U_N(t),
\end{multline*}
which is equivalent to (\ref{eq:derivative-UN}).
\end{proof}

\subsection{Derivation of the effective evolution} \label{subsec:Bogoliubov-dynamics} Let $\Psi_{N}(t)$ satisfy the Schr\"odinger equation (\ref{eq:Schrodinger-dynamics}). Let $u(t)$ be the solution to the nonlinear Hartree equation and let $U_{N}(t)$ be as in (\ref{eq:def-unitary-UN}). Let us consider the vector
$$\boxed{\Phi_{N}(t):= U_{N}(t) \Psi_{N}(t)}$$
in the Fock space $\mathcal{F}_+^{\le N}(t) \subset \mathcal{F}_+(t)\subset \mathcal{F}$. From (\ref{eq:Schrodinger-dynamics}) and (\ref{eq:derivative-UN}), we get
\begin{multline} \label{eq:derivative-Phi-Nt}
i\dot{\Phi}_{N}(t)=\Big( U_{N}(t) H_N U_{N}(t)^* + a^\dagger(u(t))a(Q(t)i\dot u(t))-\pscal{i\dot u(t),u(t)}(N-\mathcal{N}_+(t))\\ -\sqrt{N-\mathcal{N}_+(t)}\,a(Q(t)i\dot u(t))-a^\dagger(Q(t)i\dot u(t))\sqrt{N-\mathcal{N}_+(t)}  \Big) \Phi_{N}(t).
\end{multline}
A lengthy but straightforward computation (see \cite[Section 4]{LewNamSerSol-13}) shows that
\begin{flalign} \label{eq:UNHNUN-original}
& U_{N}(t)H_{N} U_{N}(t)^* = N e(t) + {\ensuremath{\rm d}\Gamma} \Big( Q(t) \big(h(t)+K_1(t)-e(t)\big) Q(t) \Big)\nonumber\\
&\quad+ \sqrt{N-\mathcal{N}_+(t)} a \big( Q(t)h(t) u(t)\big)+a ^*\big( Q(t) h(t)u(t)\big)\sqrt{N-\mathcal{N}_+(t)} \nonumber\\
&\quad+\frac12\iint_{\mathbb{R}^d\times\mathbb{R}^d}\bigg( K_2(t,x,y)a^\dagger(x)a^\dagger(y)+\overline{K_2(t,x,y)}a(x)a(y)\bigg)dx\,dy \nonumber\\
&\quad + R_{1,N}(t) + R_{2,N}(t)
\end{flalign}
where
\begin{align} \label{eq:RN1}
R&_{1,N} (t)= \frac{1}{2} {\ensuremath{\rm d}\Gamma}( Q(t) [w*|u(t)|^2+K_{1}(t)-\mu(t)]Q(t))\frac{1-\mathcal{N}_+(t)}{N-1} \nonumber\\
&- a \big(Q(t) [w*|u(t)|^2]u(t)\big)\frac{\mathcal{N}_+(t) \sqrt{N-\mathcal{N}_+(t)}}{N-1} \nonumber\\
&+  \frac{1}{2} \iint dx\, dy\, K_2(t,x,y)a^\dagger(x) a^\dagger(y)  \bigg( \frac{\sqrt{(N-\mathcal{N}_+(t))(N-\mathcal{N}_+(t)-1)}}{N-1} -1 \bigg)\nonumber \\
& + \frac{\sqrt{N-\mathcal{N}_+(t)}}{N-1} \iiiint  \, \big( \1 \otimes \overline{Q(t)} \, w \,
Q(t) \otimes Q(t) \big) (x,y; x', y') \nonumber\\ & \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \times u (t,x) \, a^\dagger(y) a(x')a(y')\,dx\, dy\, dx'\, dy' \nonumber \\
&+\text{h.c.}
\end{align}
and
\begin{align} \label{eq:RN2}
R_{2,N} (t) =  \frac{1}{2(N-1)} \iiiint & \big( \overline{Q(t)}\otimes \overline{Q(t)} \, w \, Q(t) \otimes Q(t) \big) (x,y; x',y') \nonumber \\ & \times a^\dagger(x) a^\dagger(y) a(x')a(y')\,dx\, dy\, dx'\, dy' .
\end{align}
Here $w$ is always understood as the multiplication operator by the function $w(x-y)$ in $\mathfrak{H}^2$ and recall that
\begin{align*}
P(t) & =1-Q(t)=|u(t) \rangle \langle u(t)|,\\
h(t) &=-\Delta+|u(t)|^2*w-\mu(t)\\
e(t) &=\langle u(t), h(t) u(t) \rangle = \langle u(t), (-\Delta + |u(t)|^2*w/2) u(t) \rangle.
\end{align*}
Now we substitute (\ref{eq:UNHNUN-original}) into (\ref{eq:derivative-Phi-Nt}) and use the Hartree equation $i\dot u(t) =h(t) u(t)$. Since
$$
 a^\dagger(u(t))a(Q(t) h(t) u(t))+d\Gamma( Q(t) h(t) Q(t))= {\ensuremath{\rm d}\Gamma}( h(t) ) - {\ensuremath{\rm d}\Gamma}( h(t) P(t)) ,
$$
and $\Phi_N(t)\in \mathcal{F}_+(t)$, we find that
\begin{eqnarray} \label{eq:equation_Phi_N_t}
\boxed{i\dot{\Phi}_{N}(t) =  \Big(  \mathbb{H}(t) + R_{1,N}(t)+R_{2,N}(t) \Big)\Phi_{N}(t)}
\end{eqnarray}
where $\mathbb{H}(t)$ is the Bogoliubov Hamiltonian given in (\ref{eq:Bogoliubov-Hamiltonian}):
\begin{align*}
\mathbb{H}(t)= d\Gamma(h(t)+K_1(t))+ \frac12\iint\Big(  K_2(t,x,y)&a^\dagger(x)a^\dagger(y)\\
&+\overline{K_2(t,x,y)}a(x)a(y)\Big)\, dx\,dy.
\end{align*}

Heuristically, if most particles live in the state $u(t)$, then $\mathcal{N}_+(t)/N$ is negligible and hence  $R_{1,N}(t)+R_{2,N}(t)$ can be ignored. The main goal of our paper is to show that ${\Phi}_{N}(t)$ actually converges to the exact solution of the Bogoliubov equation
\begin{eqnarray} \label{eq:formal-Bogoliubov-equation}
\boxed{i\dot{\Phi}(t) = \mathbb{H}(t) \Phi(t).}
\end{eqnarray}

The convergence $\Phi_N(t)\to \Phi(t)$ will be reformulated and justified in Section \ref{sec:reformulation}. In the rest of this section, we shall study the well-posedness of equation (\ref{eq:formal-Bogoliubov-equation}).

\subsection{Well-posedness of the Bogoliubov dynamics}\label{sec:Dynamics-quadratic-form} Now we state and prove a result concerning the well-posedness of the Bogoliubov dynamics (\ref{eq:formal-Bogoliubov-equation}). The main difficulty is of course the fact that $\mathbb{H}(t)$ is time-dependent.

\begin{theorem}[Bogoliubov dynamics] \label{thm:Bogoliubov-dynamics} Assume that $w^2 \le C(1-\Delta)$. Let $u_0\in H^1(\mathbb{R}^d)$ with $\|u_0\|=1$ and let $u(t)$ be the unique solution to the Hartree equation~\eqref{eq:Hartree-equation}. For every state $\Phi_0$ in the quadratic form domain $Q(d\Gamma(1-\Delta))$, there exists a unique solution to the Bogoliubov equation
\begin{eqnarray} \label{eq:Bog-eq}
\left\{ \begin{gathered}
   i \dot\Phi(t)=\mathbb{H}(t) \Phi(t),\hfill \\
   \Phi(0)=\Phi_0\hfill \\
\end{gathered}  \right.
\end{eqnarray}
such that $\Phi\in C^0([0,\ii),\mathcal{F})\cap L^\ii_{\rm loc}(0,\ii;Q(d\Gamma(1-\Delta)))$.
Moreover, $\Phi(t)\in \mathcal{F}_+(t)$ and
$$\langle \Phi(t), d\Gamma(1-\Delta) \Phi(t) \rangle \le C e^{Ct} \langle \Phi_0, d\Gamma(1-\Delta) \Phi_0 \rangle .$$
\end{theorem}

The proof of Theorem \ref{thm:Bogoliubov-dynamics} is based on the following abstract result, which is inspired by works of Kisy\'nski \cite{Kis-64} and Simon~\cite{Simon-71}.

\begin{theorem}[Dynamics generated by quadratic forms]\label{thm:dynamics-unbounded} Let $A\ge 1$ be a self-adjoint operator and $\{H(t)\}_{t\in [0,1)}$ a family of symmetric quadratic forms on a Hilbert space $\mathcal{H}$. Assume that
\begin{itemize}
\item[(a)] There exists a positive operator $B$ which commutes with $A$ and
$$ CA \ge H(t) \ge C^{-1}A-CB \quad{and}\quad i[H(t),B] \le CA.$$

\item[(b)] The time-derivative of $H(t)$ exists and it is bounded in $Q(A)$:
$$-C \langle x, A x\rangle \le \frac{d}{dt} \langle x, H(t) x \rangle \le C\langle x, A x\rangle\quad{for~all~}x\in Q(A).$$
\end{itemize}
Then for every $x_0\in Q(A)$, there exists a unique $x(t)\in L^\ii(0,1;Q(A))$ such that $t\mapsto x(t)$ is continuous in $\mathcal{H}$ and
$$
\left\{ \begin{gathered}
   i\dot{x}(t)= H(t)x(t)\quad {in}~~Q(A)^*,\hfill \\
   x(0)=x_0.\hfill \\
\end{gathered}  \right.
$$
Moreover, for all $t\in [0,1)$ we have $\|x(t)\|=\|x_0\|$ and
$$\langle x(t), A x(t) \rangle \le Ce^{Ct} (\langle x_0, A x_0 \rangle +C).$$
\end{theorem}

In Theorem \ref{thm:dynamics-unbounded}, the condition (b) is similar to the assumption on the boundedness of the commutator $i[H(t),A]$, which is necessary to control the energy of the dynamics generated by a time-dependent Hamiltonian. A proof of Theorem \ref{thm:dynamics-unbounded} can be found in Appendix \ref{app:dynamics}.

In order to apply Theorem \ref{thm:dynamics-unbounded}, we need the following bounds on the Bogoliubov Hamiltonian.

\begin{lemma}[Bogoliubov Hamiltonian] \label{le:Bog-Hamiltonian} Assume that $w^2 \le C(1-\Delta)$. Let $u_0\in H^1(\mathbb{R}^d)$ with $\|u_0\|=1$ and let $u(t)$ be the unique solution to the Hartree equation~\eqref{eq:Hartree-equation}. Then the Bogoliubov Hamiltonian $\mathbb{H}(t)$ defined in (\ref{eq:Bogoliubov-Hamiltonian}) satisfies
\begin{align} \label{eq:a-Bog-1}
C d\Gamma(1-\Delta) \ge \mathbb{H}(t) \ge {\ensuremath{\rm d}\Gamma} (-\Delta-C).
\end{align}
on Fock space $\mathcal{F}$. Moreover, \[
\pm \dot{\mathbb{H}}(t) \le C d\Gamma(1-\Delta) \quad \text{and }  \pm i[\mathbb{H} (t),\mathcal{N} ]\le C (\mathcal{N}+1). \]
\end{lemma}

Assuming Lemma \ref{le:Bog-Hamiltonian} for the moment, we can give the

\begin{proof}[Proof of Theorem \ref{thm:Bogoliubov-dynamics}] Thanks to Lemma \ref{le:Bog-Hamiltonian}, we can apply Theorem \ref{thm:dynamics-unbounded} with $H(t)=\mathbb{H}(t)$, $A=C{\ensuremath{\rm d}\Gamma}(1-\Delta)$ and $B=\mathcal{N}+1$. Thus we obtain all desired results on $\Phi(t)$ except the fact that $\Phi(t)\in \mathcal{F}_+(t)$. To prove that $\Phi(t)\in \mathcal{F}_+(t)$, we compute
\begin{eqnarray*}
 \frac{d}{dt} \| a(u(t))\Phi(t)\|^2 &=& 2 \Re \left \langle  a(u(t))\Phi(t), \frac{d}{dt} \big(a(u(t)) \Phi(t) \big) \right\rangle \hfill\\
&=& 2 {\rm Im} \left \langle  a(u(t))\Phi(t), \big( -a(i\dot u(t)) + a(u(t)) \mathbb{H}(t) \big) \Phi(t) \right\rangle.
\end{eqnarray*}
By using the Hartree equation $i\dot{u}_t=h(t) u(t)$ and the commutator relation
\begin{eqnarray*}
[a(u(t)),\mathbb{H}(t)] = [a(u(t)),{\ensuremath{\rm d}\Gamma}(h(t)) ]= a(h(t) u(t))
\end{eqnarray*}
we get
\begin{eqnarray*}
\frac{d}{dt} \| a(u(t))\Phi(t)\|^2 &=& 2{\rm Im} \left \langle  a(u(t))\Phi(t), \mathbb{H}(t) a(u(t)) \Phi(t) \right\rangle =0.
\end{eqnarray*}
Since $a(u_0)\Phi_0=0$, we conclude that $a(u(t))\Phi(t)=0$, and hence $\Phi(t)\in \mathcal{F}_+(t)$ for all $t\ge 0$.
\end{proof}

To prove Lemma \ref{le:Bog-Hamiltonian}, we need the following technical lemma.

\begin{lemma}[Kernel estimates] \label{le:kernel} For $j=1,2$, the kernels $K_j(t,x,y)$ are bounded uniformly in time:
$$\|K_j(t,\cdot,\cdot)\|_{L^2(\mathbb{R}^d \times \mathbb{R}^d)}\le C, \quad \|\dot{K}_j(t,\cdot,\cdot)\|_{H^{-1}(\mathbb{R}^d \times \mathbb{R}^d)} \le C.$$
\end{lemma}

Recall that $K_1(t)={Q(t)} \tilde K_1(t) Q(t)$ with $\tilde{K}_1(t,x,y)=u(t,x)w(x-y)\overline{u(t,y)}$ and $K_2(t)=(Q(t) \otimes Q(t)) \tilde K_2(t)$ with $\tilde K_2(t,x,y)=u(t,x)w(x-y)u(t,y)$.

\begin{proof} From the assumption $w^2 \le C(1-\Delta)$ and from the uniform bound $\|u(t)\|_{H^1(\mathbb{R}^d)} \le C$, it follows that $\| w^2*|u(t)|^2\|_{L^\infty} \le C$. Thus
\begin{align*} \int_{\mathbb{R}^d\times \mathbb{R}^d} |u(t,x)|^2 w(x-y)^2 |u(t,y)|^2 dx\,dy = \langle u(t), (w^2*|u(t)|^2)u(t) \rangle \le C
\end{align*}
for $j=1,2$. The desired bounds on $K_j$'s follow immediately.

Next, the time-derivatives $\dot{K}_j(t,x,y)$'s can be decomposed into many terms using the Hartree equation $i\dot{u}(t)=(-\Delta +|u(t)|^2*w-\mu(t))u(t)$. Let us consider the term $(\Delta_x u)(t,x)w(x-y) u(t,y)$ in detail.  Let $v\in H^1(\mathbb{R}^d \times \mathbb{R}^d)$. As in \cite[Lemma 6.2]{CheLeeSch-11}, by using an integration by part and the fact that $\nabla_x (w(x-y))=-\nabla_y (w(x-y))$, we get
\begin{align*} & \iint_{\mathbb{R}^d\times \mathbb{R}^d} (-\Delta_x u)(t,x)u(t,y) w(x-y) v(x,y) \, dx \, dy \nonumber\\
=&\iint_{\mathbb{R}^d\times \mathbb{R}^d} (\nabla_x u)(t,x)u(t,y) w(x-y) [\nabla_x v(x,y)+ \nabla_y v(x,y)]\, dx \, dy \nonumber\\
&+ \iint_{\mathbb{R}^d\times \mathbb{R}^d} (\nabla_x u)(t,x) (\nabla_y u)(t,y) w(x-y) v(x,y) \, dx \, dy.  \nonumber
\end{align*}
From the estimates $\| w^2*|u(t)|^2 \|_{L^\infty} \le C$, $\| w(x-y)v(x,y)\|_{L^2(\mathbb{R}^d \times \mathbb{R}^d)} \le C\|v\|_{H^1(\mathbb{R}^d \times \mathbb{R}^d)}$ and the Cauchy-Schwarz inequality we can conclude that
\begin{align*} \Big| \iint (-\Delta_x u)(t,x)u(t,y) w(x-y) v(x,y) \, dx \, dy \Big| \le  C \| v \|_{H^1}.
\end{align*}
Thus $(-\Delta_x u)(t,x)u(t,y) w(x-y)$ is bounded in $H^{-1}(\mathbb{R}^d \times \mathbb{R}^d)$. Since the other terms can be treated similarly, we obtain the desired bounds on $\dot{K}_j(t,x,y)$'s.
\end{proof}

Now we are able to give the

\begin{proof}[Proof of Lemma \ref{le:Bog-Hamiltonian}] First, we prove that,
\begin{eqnarray} \label{eq:a-Bog-1b}
C d\Gamma(1-\Delta) \ge \mathbb{H}(t) \ge {\ensuremath{\rm d}\Gamma} (1-\Delta)-C (\mathcal{N} +1)
.\end{eqnarray}
The term $d\Gamma(h(t))$ can be treated easily using the uniform bound $\| |u(t)|^2*w\|_{L^\infty} \le C$ and $\|K_1(t)\|_{L^2\to L^2}\le C$. Now let us consider the pairing term
$$\frac12\iint_{\mathbb{R}^d\times\mathbb{R}^d}\bigg(K_2(t,x,y) a^\dagger(x)a^\dagger(y) + \overline{K_2(t,x,y)} a(x)a(y) \bigg) \,dx \, dy .$$
For every state $\psi=(\psi_n)_{n\ge 0}$ in Fock space $\mathcal{F}=\bigoplus_{n\ge 0} \mathfrak{H}^n $ we can estimate
\begin{align} \label{eq:paring-K2-start}
& \bigg| \Big\langle \psi , \iint_{\mathbb{R}^d\times\mathbb{R}^d} dx dy \, \overline{{K_2}(t,x,y)} a(x)a(y) \, \psi \Big\rangle \bigg| \nonumber\\
& = \bigg| \sum_{n\ge 0} \sqrt{(n+1)(n+2)} \iint_{\mathbb{R}^d\times\mathbb{R}^d}  \,dx \,dy \int_{\mathbb{R}^{dn}} d{\bf z}  \overline{{K_2}(t,x,y)} \overline{\psi_{n}({\bf z})}\psi_{n+2}(x,y,{\bf z}) \bigg| \nonumber\\
&\le \sum_{n\ge 0} \sqrt{(n+1)(n+2)} \|K_2(t,.,.) \|_{L^2(\mathbb{R}^d \times \mathbb{R}^d)} \| \psi_n \| \| \psi_{n+2} \| \nonumber\\
& \le \frac{1}{2}\|K_2(t,.,.) \|_{L^2} \Big( \sum_{n\ge 0} (n+1)\| \psi_n \|^2 +  \sum_{n\ge 0} (n+2) \| \psi_{n+2} \|^2 \Big) \nonumber\\
&\le  \|K_2(t,.,.) \|_{L^2} \big\langle \psi, (\mathcal{N}+2) \psi \big \rangle.
\end{align}
Here we have introduced the variable ${\bf z}\in \mathbb{R}^{dn}$ and used the Cauchy-Schwarz inequality. From the kernel bound $\|K_2(t,.,.) \|_{L^2} \le C$ in Lemma \ref{le:kernel}, we can control the pairing term as
\begin{align} \label{eq:paring}
\pm \frac12\iint_{\mathbb{R}^d\times\mathbb{R}^d}\Big(K_2(t,x,y) a^\dagger(x)a^\dagger(y) +\text{h.c.}\Big) \, dx \, dy \leq C (\mathcal{N}+1).
\end{align}
Thus (\ref{eq:a-Bog-1}) follows immediately.

The time-derivative $\dot{\mathbb{H}}(t)$ can be estimated in the same way, using the kernel bound $\|\dot{K}_j(t,.,.) \|_{H^{-1}} \le C$ instead of $\|K_j(t,.,.) \|_{L^2} \le C$. More precisely,  similarly to (\ref{eq:paring-K2-start})   we have now
\begin{align*
& \bigg| \Big\langle \psi , \iint_{\mathbb{R}^d\times\mathbb{R}^d} dx \, dy \, \overline{{\dot{K_2}}(t,x,y)} a(x)a(y) \, \psi \Big\rangle \bigg| \nonumber\\
&\hspace{3cm}  \le \|\dot{K}_2(t,.,.) \|_{H^{-1}} \Big( \big\langle \psi, (\mathcal{N}+2) \psi \big \rangle + \big\langle \psi, d\Gamma(1-\Delta) \psi \big \rangle \Big)
\end{align*}
for every state $\psi\in Q(d\Gamma(1-\Delta))$. By similar bounds, we conclude that
\begin{eqnarray*}
\pm \dot{\mathbb{H}}(t) \le C d\Gamma(1-\Delta).
\end{eqnarray*}

To prove the bound for the commutator $[\mathbb{H} (t), \mathcal{N}]$ we observe that
\begin{align} \label{eq:comm-pairing}i [\mathbb{H}(t),\mathcal{N} ] &= \frac{i}{2}\iint_{\mathbb{R}^d\times\mathbb{R}^d}\bigg(K_2(t,x,y)[a^\dagger(x)a^\dagger(y),\mathcal{N} ] \nonumber\\
&\qquad\qquad\qquad\qquad+\overline{K_2(t,x,y)}[a(x)a(y),\mathcal{N} ]\bigg)dx\,dy \\
&= - i \iint_{\mathbb{R}^d\times\mathbb{R}^d} K_2(t,x,y) a^\dagger(x)a^\dagger(y) dx \, dy + \text{h.c.} \nonumber \end{align}
Similarly to (\ref{eq:paring-K2-start}), we find, for an arbitrary $\psi \in \mathcal{F}$,
\[ \begin{split} \Big| \iint_{\mathbb{R}^d \times \mathbb{R}^d} K_2 (t,x,y) \,  \langle \psi,  \,a^\dagger(x) &a^\dagger(y)  \psi \rangle \, dx \, dy \Big| \\ & \leq\; C \| K_2 (t,.,.) \|_{L^2} \| (\mathcal{N}+1)^{1/2} \psi \|^2 .\end{split} \]
Since  $\| K_2 (t,.,.) \|_{L^2} \leq C$ by Lemma \ref{le:kernel}, we can conclude that $\pm [i\mathbb{H} (t),\mathcal{N}] \leq C (\mathcal{N}+1)$.
\end{proof}

\section{Convergence of the fluctuations}\label{sec:reformulation}

\subsection{Statement involving $U_N(t)$}

Now we reformulate the convergence in Theorem \ref{thm:fluctuations}. Let $w^2\le C(1-\Delta)$. Let $\Psi_{N}(t)=e^{itH_N}\Psi_{N,0}$ be the Schr\"odinger evolution with initial wave function $\Psi_{N,0}\in \bigotimes_{\rm sym}^N H^1(\mathbb{R}^d)$. Let $u(t)$ be the unique Hartree solution with initial datum $u_0\in H^1(\mathbb{R}^d)$. Let $\Phi(t)$ satisfy the Bogoliubov equation (\ref{eq:Bog-eq}) with initial datum $\Phi_0 \in Q({\ensuremath{\rm d}\Gamma}(1-\Delta))\subset \mathcal{F}$ as in Theorem \ref{thm:Bogoliubov-dynamics} . Let $U_{N}(t)$ be defined from $u(t)$ as in (\ref{eq:def-unitary-UN}). Our main result in this section is

\begin{theorem} [Convergence to Bogoliubov dynamics] \label{thm:fluctuations_U_N}  If
\begin{eqnarray} \label{eq:cv-initial-datum}
U_{N,0} \Psi_{N,0} \to \Phi_0
\end{eqnarray}
strongly in $\mathcal{F}$ and weakly in $Q({\ensuremath{\rm d}\Gamma}(1-\Delta))$ as $N\to \infty$, then
$$\boxed{\lim_{N\to \infty} U_{N}(t) \Psi_{N}(t)= \Phi(t)}$$
strongly in $\mathcal{F}$ and weakly in $Q({\ensuremath{\rm d}\Gamma}(1-\Delta))$, for every fixed $t\geq0$.
\end{theorem}

\begin{remark} The condition (\ref{eq:cv-initial-datum}) is clearly satisfied if we start with a Hartree state at time 0, $\Psi_{N,0}=u_0^{\otimes N}$ with $u_0\in H^1(\mathbb{R}^d)$. In this case we have $U_{N,0}\Psi_{N,0}=|0 \rangle = \Phi_0$ for all $N$.
\end{remark}

In Theorem~\ref{thm:fluctuations_U_N}, $\Phi(t)$ is a state on $\mathcal{F}_+(t)\subset \mathcal{F}$ which describes the fluctuations around the Hartree state $u(t)$ through the unitary $U_N(t)$. Note that
\begin{align*}
\lim_{N\to \infty} \| U_N(t) \Psi_N(t) - \Phi(t)\|_{\mathcal{F}} &=\lim_{N\to \infty} \| U_N(t) \Psi_N(t) - \1_+^{\le N}\Phi_N(t)\|_{\mathcal{F}}\\
&=\lim_{N\to \infty} \|  \Psi_N(t) - U_N(t)^* \1_+^{\le N} \Phi(t)\|_{\mathfrak{H}^N} \\
&=\lim_{N\to \infty}\norm{\Psi_N(t)-\sum_{n=0}^N u(t)^{\otimes (N-n)}\otimes_s\phi_n(t)}_{\mathfrak{H}^N}
\end{align*}
where $\1_{+}^{\le N}$ is the projection onto the truncated Fock space $\mathcal{F}_+(t)^{\le N}$. Therefore, Theorem \ref{thm:fluctuations_U_N} implies Theorem \ref{thm:fluctuations} immediately. Theorem \ref{thm:fluctuations_U_N} will be proved in the rest of this section.

\subsection{Bound on number of particles outside of the condensate} The goal of this subsection is to estimate the expectation of the number of particles operator, and of its powers, in the evolved many body state $\phi_N (t) = U_N (t) \Psi_{N} (t)$. These estimates will play a crucial role in the proof of Theorem \ref{thm:fluctuations_U_N}. To control the growth of the expectation of the number of particles operator, we need the following bounds on the remainder $R_{N,1} (t)$ defined in (\ref{eq:RN1}).

\begin{lemma}[Remainder Hamiltonian] \label{le:remainder-bound} Let us denote $\1_{+}^{\le M}$ the projection onto the truncated Fock space $\mathcal{F}_+(t)^{\le M}$. For all $1\le M\le N$, we have
\begin{align} \label{eq:bound-ENt}
\pm \1_{+}^{\le M} R_{1,N} (t) \1_{+}^{\le M} \le C \sqrt{\frac{M}{N}}  (\mathcal{N} +1).
\end{align}
Moreover
\begin{equation}\label{eq:bound-coR} \pm \1_{+}^{\le N} i\left[ R_{1,N} (t),\mathcal{N} \right] \1_{+}^{\le N} \le C \, (\mathcal{N} +1). \end{equation}
\end{lemma}

\begin{proof} The first estimate (\ref{eq:bound-ENt}) was already proved in \cite[Proposition 15]{LewNamSerSol-13}.

Now we consider the commutator $i[R_{1,N} (t),\mathcal{N} ]$. A straightforward computation similar to \eqref{eq:comm-pairing} shows that
\begin{align} \label{eq:comm-RN}
& i[R_{1,N}(t),\mathcal{N}] \nonumber \\
 &= - i a \big(Q(t) [w*|u(t)|^2]u(t)\big) \frac{\mathcal{N}_+(t) \, \sqrt{N-\mathcal{N}_+(t)}}{N-1} \nonumber\\
& +2i\bigg( \frac{\sqrt{(N-\mathcal{N}_+(t))(N-\mathcal{N}_+(t)-1)}}{N-1} -1 \bigg)\!
\iint \overline{K_2(t,x,y)}a(x) a(y)\,dx \,dy \,  \nonumber\\
& + i\frac{\sqrt{N-\mathcal{N}_+(t)}}{N-1}  \iiiint  \big( \1 \otimes \overline{Q(t)} \, w \, Q(t) \otimes Q(t) \big) (x,y; x',y')
\nonumber\\ &\hspace{5cm} \times
 u (t,x) a^\dagger(y) a(x')a(y')\,dx\, dy\, dx'\, dy' \nonumber\\
&+ \text{h.c.}
\end{align}
Recall that, in the last term, $w$ is understood as the multiplication operator by the function $w(x-y)$.

We shall show that all terms in (\ref{eq:comm-RN}), when projected onto $\mathcal{F}_+(t)^{\le N}$, can be bounded by $C(\mathcal{N}+1)$. Let $v(t) = i Q(t) [w*|u(t)|^2] u(t)$. Since $\| v (t) \|_{L^2} \leq \| w * |u (t)|^2 \|_{L^\infty} \| u (t) \|_{L^2} \leq C$, the first term in (\ref{eq:comm-RN}) can be bounded, after projecting onto $\mathcal{F}_+(t)^{\le N}$, by
\[ \pm \left\{ a (v(t)) \frac{\mathcal{N}_+(t) \, \sqrt{N-\mathcal{N}_+(t)}}{N-1}  + \text{h.c.}\right\}  \leq C (\mathcal{N}+1) \]
The second term on the right side of (\ref{eq:comm-RN}) can be treated similarly to the pairing term of $\mathbb{H}(t)$ in (\ref{eq:comm-pairing}). As for the third term on the right side of (\ref{eq:comm-RN}), we write
\begin{align}\label{eq:cub-1}
&\frac{\sqrt{N-\mathcal{N}_+}}{N-1} \iint\!\!\! \iiint \,dx \,dy \,dy_1 \,dx' \,dy' Q_t (y,y_1) w (x-y_1) Q_t (x;x') Q_t (y_1;y') \nonumber\\ &\hspace{7cm} \times u (t,x) a^\dagger (y) a (x') a (y') \nonumber\\
&=\frac{\sqrt{N-\mathcal{N}_+}}{N-1} \iiint \,dx \,dy \,dy'  w (x-y) Q_t (y;y') \, u (t,x) a^\dagger (Q_{t,y}) a (Q_{t,x}) a (y') \nonumber\\
&=\frac{\sqrt{N-\mathcal{N}_+}}{N-1} \iint \,dx \,dy\,  w (x-y)  \, u (t,x) a^\dagger (Q_{t,y}) a (Q_{t,x}) a (y) \\
&\hspace{.4cm} - \frac{\sqrt{N-\mathcal{N}_+}}{N-1} \iint \,dx \,dy\,    w (x-y) \, u (t,y) \, u (t,x) \, a^\dagger (Q_{t,y}) a (Q_{t,x}) a (u (t)) \nonumber
\end{align}
where we introduced the notation $Q_t (x,x')$ to denote the kernel of the operator $Q(t) = 1-|u (t) \rangle \langle u(t)|$, and where $Q_{t,x} (z) = Q_t (z,x)$. To bound the first term on the right side of (\ref{eq:cub-1}), we notice that, for arbitrary $\psi \in \mathcal{F}_+^{\le N}$,
\[ \begin{split}
\Big|  \iint & dx \, dy   \, w (x-y) \, u (t,x) \, \left\langle \psi,  \sqrt{N-\mathcal{N}_+} \, a^\dagger (Q_{t,y}) a (Q_{t,x}) a (y) \psi \right\rangle \Big| \\
&\leq \iint dx\, dy \, |w (x-y)| \, |u (t,x)| \left\| a(Q_{t,y}) \sqrt{N-\mathcal{N}_+} \psi \right\| \, \left\| a (y) a (Q_{t,x}) \psi \right\| .
\end{split} \]
By the Cauchy-Schwarz inequality we find
\begin{align*}
\Big|  \iint dx\, &dy   \, w (x-y) \, u (t,x) \, \left\langle \psi,  \sqrt{N-\mathcal{N}_+} \, a^\dagger (Q_{t,y}) a (Q_{t,x}) a (y) \psi \right\rangle \Big| \\
&\leq \iint dx dy \, |w (x-y)|^2 \, |u (t,x)|^2 \left\| a(Q_{t,y}) \sqrt{N-\mathcal{N}_+} \psi \right\|^2 \\
&\hspace{.4cm} + \iint dx\, dy \, \left\| a (y) a (Q_{t,x}) \psi \right\|^2 \\
&\leq C \int dy  \left\| a(Q_{t,y}) \sqrt{N-\mathcal{N}_+} \psi \right\|^2 + \int dx \, \left\| a (Q_{t,x}) (\mathcal{N}+1)^{1/2} \psi \right\|^2 \\
&\leq CN \langle \psi, d\Gamma (Q (t)) \psi \rangle \leq C N \langle \psi , \mathcal{N} \psi \rangle .
\end{align*}
As for the second term on the right side of (\ref{eq:cub-1}), we have, for an arbitrary $\psi \in \mathcal{F}_+^{\le N}$,
\[ \begin{split}
\Big| \iint dx\, dy\, &w(x-y) \, u(t,x) u (t,y) \left\langle \psi, \sqrt{N-\mathcal{N}_+} \, a^\dagger (Q_{t,y}) a(Q_{t,x}) a(u_t) \psi \right\rangle \Big| \\
&\leq \iint dx\, dy \, |w (x-y)| |u (t,x)| |u (t,y)| \, \left\| a(Q_{t,y}) \sqrt{N-\mathcal{N}_+} \psi \right\| \\ &\hspace{7cm} \times  \left\| a(Q_{t,x}) a(u_t) \psi \right\| .
 \end{split} \]
By the Cauchy-Schwarz inequality again, we find
\[ \begin{split}
\Big| \int dx \,dy\, &w(x-y) \, u(t,x) u (t,y) \left\langle \psi, \sqrt{N-\mathcal{N}_+} \, a^\dagger (Q_{t,y}) a(Q_{t,x}) a(u_t) \psi \right\rangle \Big| \\
&\leq \int dx dy \, |w (x-y)|^2 |u (t,x)|^2 \, \left\| a(Q_{t,y}) \sqrt{N-\mathcal{N}_+} \psi \right\|^2 \\ &\hspace{.4cm} +\int dx dy \, |u (t,y)|^2 \,  \left\| a(Q_{t,x}) a(u_t) \psi \right\|^2 \\
&\leq C N \, \langle \psi, \mathcal{N} \psi \rangle . \end{split} \]
Inserting in (\ref{eq:cub-1}), we conclude that the third term on the right side of (\ref{eq:comm-RN}), after projecting onto $\mathcal{F}_+^{\le N}$, is also bounded by $C (\mathcal{N}+1)$.
\end{proof}

As an easy consequence of Lemma \ref{le:remainder-bound}, we have
\begin{lemma}[Number of particles outside of the condensate] \label{le:N+-bound} There exists a constant $C$ such that
$$\langle \Phi_N(t), (\mathcal{N}+1)  \Phi_N(t) \rangle \le C e^{Ct} \langle \Phi_N(0), (\mathcal{N}+1) \Phi_N(0) \rangle.$$
\end{lemma}
\begin{proof} From the equation $i\dot{\Phi}_N(t)=(\mathbb{H}(t)+R_N(t)) \Phi_N(t)$ and the commutator bounds in Lemma \ref{le:Bog-Hamiltonian} and Lemma \ref{le:remainder-bound}, we have
\begin{align*} \frac{d}{dt} \langle \Phi_N(t), (\mathcal{N} +1) \Phi_N(t) \rangle &= \big\langle \Phi_N(t), i[\mathbb{H}(t)+R_N(t),\mathcal{N}+1] \Phi_N(t) \big\rangle \\
&= \big\langle \Phi_N(t), i[R_{1,N}(t),\mathcal{N}] \Phi_N(t) \big\rangle \\
&\le C \langle \Phi_N(t), (\mathcal{N} +1) \, \Phi_N(t) \rangle.
\end{align*}
The desired bound then follows from Gronwall's inequality.
\end{proof}

\subsection{Proof of the main result}

Now we are able to prove Theorem \ref{thm:fluctuations_U_N} which itself implies Theorem~\ref{thm:fluctuations}.

\begin{proof} We have
\begin{align}\label{eq:diff-1}
\frac{d}{dt} &\| \Phi_{N}(t)-\Phi(t)\|^2 = 2\Re \Big\langle \Phi_{N}(t)-\Phi(t), \dot{\Phi}_{N}(t)-\dot \Phi(t)\Big\rangle \nonumber\\
&= -2 {\rm Im} \Big\langle \Phi_{N}(t)-\Phi(t), (\mathbb{H}(t) +R_N(t)){\Phi}_{N}(t)-\mathbb{H}(t) {\Phi}(t)\Big\rangle \nonumber\\
&= -2 {\rm Im} \Big\langle \Phi_{N}(t), R_N(t) {\Phi}(t)\Big\rangle \nonumber\\
&=-2  {\rm Im} \Big\langle \Phi_{N}(t), R_{1,N} (t) {\Phi}(t)\Big\rangle -2  {\rm Im} \Big\langle \Phi_{N}(t), R_{2,N} (t) {\Phi}(t)\Big\rangle
\end{align}
where $R_{1,N}(t)$ and $R_{2,N}(t)$ are defined in (\ref{eq:RN1}) and (\ref{eq:RN2}), respectively. We shall show that both terms on the right side of (\ref{eq:diff-1}) tend to $0$ as $N\to \infty$. To this end, we denote by $\1_+^{\leq M}$ the orthogonal projection onto the truncated Fock space $\mathcal{F}_+(t)^{\leq M}$ and
$$\Phi(t)^{\le M}:=\1_+^{\le M} \Phi(t), \quad \Phi(t)^{> M}:=(\1-\1_+^{\le M}) \Phi(t) =\Phi(t) -\Phi(t)^{\le M}.$$
{\bf Estimate on the first remainder term.} We start by estimating the first term on the right side of (\ref{eq:diff-1}). Since $R_{1,N}$ can change the number of particles at most by two, for any $M \in \mathbb{N}$ we have
\[ \begin{split}
\langle \Phi_N (t) , R_{1,N} (t) \Phi (t) \rangle = \, &\langle \Phi_N (t), R_{1,N} (t) \Phi (t)^{\leq M} \rangle + \langle \Phi_N (t), R_{1,N} (t) \Phi (t)^{> M} \rangle \\
= \; &\langle \Phi_N (t), \1_+^{\leq (M +2)} R_{1,N} (t) \1_+^{\leq (M+2)} \Phi (t)^{\leq M} \rangle \\ &+ \langle \Phi_N (t), \1^{\leq (N+2)} \, R_{1,N} (t)  \1^{\leq (N+2)} \Phi (t)^{> M} \rangle.
\end{split} \]
By Lemma \ref{le:remainder-bound}, we can choose $C>0$ so large that
\[ \begin{split} \1_+^{\leq (M +2)} R_{1,N} (t) \1_+^{\leq (M+2)} + C \sqrt{\frac{M}{N}} (\mathcal{N}+1) &\geq 0 \\ \1_+^{\leq (N+2)} R_{1,N} (t) \1_+^{\leq (N+2)} + C (\mathcal{N}+1) &\geq 0.
\end{split} \]
Then
\[ \begin{split}
\langle \Phi_N &(t) , R_{1,N} (t) \Phi (t) \rangle \\ = \; &\left\langle \Phi_N (t), \left[ \1_+^{\leq (M +2)} R_{1,N} (t) \1_+^{\leq (M+2)}+ C \sqrt{\frac{M}{N}} \, (\mathcal{N}+1) \right] \Phi (t)^{\leq M} \right\rangle \\ &- C
\sqrt{\frac{M}{N}} \, \left\langle \Phi_N (t), (\mathcal{N}+1) \Phi (t)^{\le M} \right\rangle \\
&+ \left\langle \Phi_N (t), \left[ \1^{\leq (N+2)} \, R_{1,N} (t)  \1^{\leq (N+2)} +C (\mathcal{N}+1) \right] \, \Phi (t)^{> M} \right\rangle \\ &- C \left\langle \Phi_N (t), (\mathcal{N}+1) \Phi (t)^{> M} \right\rangle.
\end{split} \]
By the Cauchy-Schwarz inequality, we estimate
\[ \begin{split}
\big| \langle &\Phi_N (t) , R_{1,N} (t) \Phi (t) \rangle \big| \\ \leq \; &\left\langle \Phi_N (t), \left[ \1_+^{\leq (M +2)} R_{1,N} (t) \1_+^{\leq (M+2)}+ C \sqrt{M/N} \, (\mathcal{N}+1) \right] \Phi_N (t) \right\rangle^{1/2} \\
&\hspace{.1cm} \times
\left\langle \Phi^{\le M} (t), \left[ \1_+^{\leq (M +2)} R_{1,N} (t) \1_+^{\leq (M+2)}+ C \sqrt{M/N} \, (\mathcal{N}+1) \right] \Phi^{\le M} (t) \right\rangle^{1/2} \\
&+  \left\langle \Phi_N (t), \left[ \1^{\leq (N+2)} \, R_{1,N} (t)  \1^{\leq (N+2)} +C (\mathcal{N}+1) \right] \, \Phi_N (t) \right\rangle^{1/2} \\
&\hspace{.1cm} \times
\left\langle \Phi (t)^{> M}, \left[ \1^{\leq (N+2)} \, R_{1,N} (t)  \1^{\leq (N+2)} +C (\mathcal{N}+1) \right] \, \Phi (t)^{> M} \right\rangle^{1/2} \\
&+C \sqrt{M/N} \, \left\langle \Phi_N (t), (\mathcal{N}+1) \Phi_N (t) \right\rangle^{1/2} \left\langle \Phi (t)^{\leq M}, (\mathcal{N}+1) \Phi (t)^{\leq M} \right\rangle^{1/2} \\
&+C \, \left\langle \Phi_N (t), (\mathcal{N}+1) \Phi_N (t) \right\rangle^{1/2} \left\langle \Phi (t)^{> M}, (\mathcal{N}+1) \Phi (t)^{> M} \right\rangle^{1/2} .
\end{split} \]
Using again Lemma \ref{le:remainder-bound} we find
\[ \begin{split}
\big| \langle \Phi_N &(t) , R_{1,N} (t) \Phi (t) \rangle \big| \\ \leq \; &C \sqrt{M/N} \, \left\langle \Phi_N (t), (\mathcal{N}+1) \Phi_N (t) \right\rangle^{1/2}  \left\langle \Phi (t)^{\le M}, (\mathcal{N}+1) \Phi^{\le M} (t) \right\rangle^{1/2} \\ &+C \left\langle \Phi_N (t), (\mathcal{N}+1) \Phi_N (t) \right\rangle^{1/2}  \left\langle \Phi (t)^{> M}, (\mathcal{N}+1) \Phi^{> M} (t) \right\rangle^{1/2}. \end{split} \]
From Lemma \ref{le:N+-bound}, we conclude that
\[ \begin{split}
\big| \langle \Phi_N &(t) , R_{1,N} (t) \Phi (t) \rangle \big| \\ \leq \; &C e^{Ct} \sqrt{M/N} \, \left\langle \Phi_N (0), (\mathcal{N}+1) \Phi_N (0) \right\rangle^{1/2}  \left\langle \Phi (0), (\mathcal{N}+1) \Phi (0) \right\rangle^{1/2} \\ &+C e^{Ct} \left\langle \Phi_N (0), (\mathcal{N}+1) \Phi_N (0) \right\rangle^{1/2}  \left\langle \Phi (t)^{> M}, (\mathcal{N}+1) \Phi^{> M} (t) \right\rangle^{1/2}. \end{split} \]
Letting first $N \to \infty$ and at the end $M \to \infty$, we obtain that
\begin{eqnarray} \label{eq:last-RN1}
\lim_{N\to \infty}\langle \Phi_N (t) , R_{1,N} (t) \Phi (t) \rangle= 0.
\end{eqnarray}
{\bf Estimate on the second remainder term.} Next, we consider the second term on the right side of (\ref{eq:diff-1}). Using the shorthand notation $A = (\overline{Q(t)} \otimes \overline{Q(t)}) \, w \, (Q(t) \otimes Q(t))$, we have
\[ \begin{split} \Big\langle \Phi_{N}(t), R_{2,N} (t) {\Phi}(t)\Big\rangle = \; &\frac{1}{N-1}\sum_{n\leq N}  \sum_{i<j}^n \left\langle \Phi^{(n)}_N (t), A_{ij} \, \Phi^{(n)} (t) \right\rangle \\ = \; &\frac{1}{N-1}\sum_{n \leq M} \frac{n(n-1)}{2} \, \left\langle \Phi^{(n)}_N (t), A_{12} \, \Phi^{(n)} (t) \right\rangle \\ &+ \frac{1}{N-1} \sum_{M \leq n \leq N} \frac{n(n-1)}{2} \, \left\langle \Phi^{(n)}_N (t), A_{12} \, \Phi^{(n)} (t) \right\rangle
\end{split}\]
where $A_{ij}$ denotes the two-particle operator $A$ acting on particles $i$ and $j$, and where we used the permutation symmetry of $\Phi_N (t)$ and $\Phi (t)$. By the Cauchy-Schwarz inequality, we find
\begin{align}\label{eq:quart-bd}
\Big|  \Big\langle \Phi_{N}(t), R_{2,N} (t) {\Phi}(t)\Big\rangle \Big| \leq  &\frac{M^{3/2}}{N} \left( \sum_{n =0}^\infty n \| A_{12} \Phi^{(n)} (t) \|_2^2 \right)^{1/2} \\ & + \langle \Phi_N (t), \mathcal{N} \Phi_N (t) \rangle^{1/2}  \left( \sum_{n =M}^\infty n \| A_{12} \Phi^{(n)} (t) \|_2^2 \right)^{1/2}.\nonumber
\end{align}
Next, we observe that
\[ \begin{split} \| A_{12} \Phi^{(n)} (t) \|^2 &= \Big\langle \Phi^{(n)} (t), A_{12}^* A_{12} \Phi^{(n)} (t) \Big\rangle \\ &= \Big\langle \Phi^{(n)} (t),  Q_1 (t) Q_2 (t) \, w (x_1 - x_2) Q_1 (t) Q_2 (t) \\ &\hspace{4cm} \times w (x_1 - x_2) Q_1 (t) Q_2 (t) \Phi^{(n)} (t) \Big\rangle \\
&\leq \Big\langle  \Phi^{(n)} (t),  Q_1 (t) Q_2 (t) \, w^2 (x_1 - x_2) Q_1 (t) Q_2 (t) \Phi^{(n)} (t) \Big\rangle
\end{split} \]
where $Q_1 (t)$ and $Q_2 (t)$ denote the orthogonal projection $Q (t)$ acting on the first and the second particles, respectively. Using the assumption $w^2 \leq C(1-\Delta)$, we find
\[ \begin{split}
 \| A_{12} \Phi^{(n)} (t) \|^2
&\leq C \Big\langle  \Phi^{(n)} (t),  Q_1 (t) Q_2 (t) \, (1-\Delta_{x_1}) Q_1 (t) Q_2 (t) \Phi^{(n)} (t) \Big\rangle \\
&\leq C \Big\langle \Phi^{(n)} (t),  Q_1 (t)  \, (1-\Delta_{x_1}) Q_1 (t) \Phi^{(n)} (t) \Big\rangle.
\end{split} \]
We write $Q_1 (t) = \1 - P_1 (t)$, with $P (t) = |u(t) \rangle \langle u(t)|$ and $P_1 (t)$ denotes the orthogonal projection $P(t)$ acting on the first particle. We easily find
\[ \begin{split}
\| A_{12} \Phi^{(n)} (t) \|^2  &\leq C \Big\langle \Phi^{(n)} (t), (1-\Delta_{x_1})  \Phi^{(n)} (t) \Big\rangle + \| u (t) \|_{H^1}^2 \| \Phi^{(n)} (t) \|_2^2 \\ &\leq C \Big\langle \Phi^{(n)} (t), (1-\Delta_{x_1})  \Phi^{(n)} (t) \Big\rangle.
\end{split} \]
Inserting the latter estimate into (\ref{eq:quart-bd}) and using Theorem \ref{thm:Bogoliubov-dynamics} and Lemma \ref{le:N+-bound}, we obtain
\begin{align*}
&\Big|  \Big\langle \Phi_{N}(t), R_{2,N} (t) {\Phi}(t)\Big\rangle \Big| \nonumber\\
\leq  &\frac{M^{3/2}}{N} \Big\langle \Phi (t) , d\Gamma (1-\Delta) \Phi (t) \Big\rangle^{1/2} \nonumber\\ &+  \Big\langle \Phi_N (t), \mathcal{N} \Phi_N (t) \Big\rangle^{1/2} \Big\langle \Phi (t)^{>M} , d\Gamma (1-\Delta) \Phi (t)^{>M} \Big\rangle^{1/2}\nonumber \\
\leq \; & C e^{Ct} \frac{M^{3/2}}{N}   \Big\langle \Phi (0) , d\Gamma (1-\Delta) \Phi (0) \Big\rangle^{1/2} \nonumber\\ &+ C e^{Ct} \Big\langle \Phi_N (0), \mathcal{N} \Phi_N (0) \Big\rangle^{1/2} \Big\langle \Phi (t)^{>M} , d\Gamma (1-\Delta) \Phi (t)^{>M} \Big\rangle^{1/2}.
\end{align*}
Under the assumption (\ref{eq:cv-initial-datum}), we have
$$\Big\langle \Phi (0) , d\Gamma (1-\Delta) \Phi (0) \Big\rangle \le C_0 \quad{\rm and}~~\Big\langle \Phi_N (0), \mathcal{N} \Phi_N (0) \Big\rangle \le C_0$$
where $C_0$ is dependent on the initial data, but independent of $N$. Therefore,
\begin{align*}
&\Big|  \Big\langle \Phi_{N}(t), R_{2,N} (t) {\Phi}(t)\Big\rangle \Big| \nonumber\\
\leq \; & C e^{Ct} \frac{M^{3/2}}{N} +C e^{Ct}   \Big\langle \Phi (t)^{>M} , d\Gamma (1-\Delta) \Phi (t)^{>M} \Big\rangle^{1/2} .
\end{align*}
For every $t>0$ fixed, since $\langle \Phi (t), d\Gamma (1-\Delta) \Phi (t) \rangle \le C$ by Theorem \ref{thm:Bogoliubov-dynamics}, it follows that
\begin{eqnarray*}
\lim_{M\to \infty} \langle \Phi^{>M} (t), d\Gamma (1-\Delta) \Phi^{>M} (t) \rangle= 0.
\end{eqnarray*}
Thus letting $N \to \infty$ and afterwards $M \to \infty$, we can conclude that
\begin{eqnarray} \label{eq:last-RN2}
\lim_{N\to \infty} \Big\langle \Phi_{N}(t), R_{2,N} (t) {\Phi}(t)\Big\rangle =0 .
\end{eqnarray}
{\bf Conclusion.} Substituting (\ref{eq:last-RN1}) and (\ref{eq:last-RN2}) into (\ref{eq:diff-1}), we obtain
$$ \lim_{N\to \infty} \frac{d}{dt} \| \Phi_{N}(t)-\Phi(t)\|^2 =0.$$
Since $\lim_{N\to \infty} \| \Phi_{N}(0)-\Phi(0)\|^2=0$ by the assumption (\ref{eq:cv-initial-datum}), we conclude that
$$\lim_{N\to \infty} \| \Phi_{N}(t)-\Phi(t)\|^2=0$$
for every $t\ge 0$ fixed. The proof is complete.
\end{proof}