# Dynamics of mean-field bosons at positive temperature

Marco Caporaletti  Andreas Deuchert  Benjamin Schlein

(August 27, 2024)

###### Abstract

We study the time evolution of an initially trapped weakly interacting Bose
gas at positive temperature, after the trapping potential has been switched
off. It has been recently shown in [24] that the one-particle density matrix
of Gibbs states of the interacting trapped gas is given, to leading order in
$N$, as $N\to\infty$, by the one of the ideal gas, with the condensate wave
function replaced by the minimizer of the Hartree energy functional. We show
that this structure is stable with respect to the many-body evolution in the
following sense: the dynamics can be approximated in terms of the time-
dependent Hartree equation for the condensate wave function and in terms of
the free evolution for the thermally excited particles. The main technical
novelty of our work is the use of the Hartree–Fock–Bogoliubov equations to
define a fluctuation dynamics.

###### Contents

  1. 1 Introduction and main results
    1. 1.1 Background and summary
    2. 1.2 Notation
    3. 1.3 The model
    4. 1.4 Equilibrium states of trapped Bose gases
    5. 1.5 Construction of the initial data
    6. 1.6 Main result
    7. 1.7 Proof strategy and organization of the article
  2. 2 Properties of the 1-pdm of the ideal gas and of Hartree minimizers
  3. 3 The Hartree–Fock–Bogoliubov equations
    1. 3.1 Well-posedness
    2. 3.2 Approximation of the HFB dynamics
  4. 4 Fluctuation dynamics
    1. 4.1 Weyl operators and Bogoliubov transformations
    2. 4.2 Construction of the fluctuation dynamics
    3. 4.3 Bound by the expectation of the number of particles
    4. 4.4 Proof of Theorem 1.2
  5. 5 Study of the fluctuation dynamics
    1. 5.1 Generator of the fluctuation dynamics
    2. 5.2 The truncated fluctuation dynamics
    3. 5.3 Weak bounds on the original dynamics
    4. 5.4 Comparison of the two fluctuation dynamics
    5. 5.5 Proof of Proposition 4.3
  6. References

##  1 Introduction and main results

###  1.1 Background and summary

The time evolution of bosonic many-particle systems in the mean-field limit
has received a considerable amount of attention in recent years. The main
theme in most works is the study of the dynamics of Bose gases that are
initially prepared in a ground state (equilibrium state at zero temperature)
of a trapped Hamiltonian after the trap has been switched off. It is well-
known, see e.g. [51, 66] and references therein, that under appropriate
conditions on the interaction potential, ground states of interacting mean-
field systems display complete Bose–Einstein condensation (BEC), that is, the
largest eigenvalue $\lambda_{N}$ of the one-particle density matrix (1-pdm) of
the ground state wave function satisfies $\lambda_{N}/N\to 1$ as $N\to\infty$.
It turns out that this structure is stable w.r.t. the dynamics generated by
the time-dependent many-particle Schrödinger equation. More precisely, it can
be shown that the system displays complete BEC also at time $t>0$, where the
condensate wave function (the eigenfunction related to the largest eigenvalue
of the 1-pdm) is given by the solution to the time-dependent Hartree equation.

The first results in this direction were obtained in [44] and then in [37, 38,
69]. The case of Coulomb interactions has been studied more recently in [5,
30] and in [25] with a relativistic dispersion relation. In [31, 36] the
convergence towards the Hartree equation has been interpreted as an Egorov
type theorem, [2] focuses on the propagation of Wigner measures, and bounds on
the rate of convergence towards the Hartree dynamics have been obtained in
[65, 49, 21, 1]. Once the convergence towards the time-dependent Hartree
equation is established, it is natural to ask whether it is possible to make
statements about the fluctuations around the Hartree dynamics. As shown in
[44, 39, 40, 6, 16, 54], for mean-field systems it is possible to approximate
the dynamics of particles outside the condensate by the time evolution
generated by a Bogoliubov Hamiltonian. This, in particular, allows one to
obtain a norm approximation for the wave function, and not only for the 1-pdm.
A systematic expansion for correlation functions in terms of the solution of
the Hartree equation and a Bogoliubov two-point function that approximates the
original dynamics to arbitrary precision has been obtained in [12].

Another interesting scaling limit for the Bose gas is the Gross–Pitaevskii
(GP) limit, which can be used to describe modern experiments with cold alkali
gases. In contrast to the mean-field limit, where particle collisions are
frequent but weak, particles interact rarely but strongly in the GP limit.
These strong collisions induce microscopic correlations between the particles
that are the main reason why the GP limit is more challenging from a
mathematical point of view than the mean-field limit. It has been shown in
[56, 55] that ground states of trapped gases display complete BEC also in the
GP limit. Moreover, in the series of works [26, 27, 28, 29] it has been
established that there is complete BEC also at later times and that the
condensate wave function is given by the solution to the time-dependent GP
equation. For other results in this direction, we refer to [48, 47, 19, 20,
17, 18, 61, 8, 45, 15, 41, 50, 42, 43, 58, 59, 11, 14, 22]. More references
concerning the dynamics of BECs in the above scaling limits can be found in
the lecture notes [9] and in the review article [60].

Here we are interested in the dynamics of mean-field bosons, initially
prepared in an equilibrium state at positive temperature (Gibbs state). Just
above the critical temperature for BEC, correlation functions of the many-body
Gibbs state have been shown to converge, in the mean-field limit, to
correlation functions of the invariant measure associated with the nonlinear
Hartree equation, see [52, 32, 53, 34, 35]; observe that, in dimensions
$d=2,3$, this requires an appropriate renormalization of the interaction. In
[33], also time-dependent correlations have been proven to converge. In the
present paper, we are primarily interested in initial data prepared at
equilibrium below the critical temperature, where the BEC coexists with a
thermal cloud and both phases contain a macroscopic number of particles. In
[24], it has been recently shown that at temperatures lower than (but
comparable with) the critical temperature of the corresponding ideal (non-
interacting) gas, the 1-pdm of the Gibbs state of a trapped mean-field system
can be approximated by the 1-pdm of the ideal gas, where the condensate wave
function has been replaced by the minimizer of the Hartree energy functional.
Our goal here is to investigate the dynamics of such positive temperature
states, after the trapping potential has been switched off. In our main
theorem we show that the 1-pdm of the evolved state is given, to leading
order, by the one of the ideal gas, propagated with the non-interacting time
evolution, with the condensate wave function replaced by the solution of the
time-dependent Hartree equation. In particular, the number of particles inside
the condensate and outside the condensate is preserved by the time evolution,
to leading order.

###  1.2 Notation

We denote by $\langle\cdot,\cdot\rangle$ the inner product in
$L^{2}(\mathbb{R}^{d})$, and by $\mathcal{B}(L^{2}(\mathbb{R}^{d}))$ the space
of bounded operators on $L^{2}(\mathbb{R}^{d})$. The $L^{2}$ norm of a
function and the norm of an operator in $\mathcal{B}(L^{2}(\mathbb{R}^{d}))$
are both denoted by $\|\cdot\|$. If $p\neq 2$, we denote by $\|\cdot\|_{p}$
the $L^{p}$ norms of a function, and
$\|A\|_{\mathcal{L}^{p}}=(\operatorname{Tr}[(A^{*}A)^{p/2}])^{1/p}$ is the
$p$-th Schatten norm of the operator $A\in\mathcal{B}(L^{2}(\mathbb{R}^{d}))$,
for $p\in[1,\infty)$. The corresponding Banach spaces of compact operators
with finite $\mathcal{L}^{p}$ norm are denoted by $\mathcal{L}^{p}$. Sobolev
spaces are denoted by $W^{n,p}$ or by $H^{n}$ if $p=2$, with
$\|\cdot\|_{W^{n,p}}$ and $\|\cdot\|_{H^{n}}$ denoting the corresponding
Sobolev norms. For the Fourier transform in $\mathbb{R}^{d}$ we use the
notation, and the convention,

$\mathcal{F}f(k)=\hat{f}(k)=\frac{1}{(2\pi)^{d/2}}\int_{\mathbb{R}^{d}}e^{\mathrm{i}k\cdot
x}f(x)\,\mathrm{d}x\,.$

For quantities $a,b$ depending on $N$ and/or $t$, we use the notation
$a\lesssim b$ to say that there exists a constant $C>0$ with $a\leq Cb$. If
$a\lesssim b$ and $b\lesssim a$ we write $a\sim b$.

###  1.3 The model

In this article we consider the time evolution of an initially trapped Bose
gas, prepared in a positive temperature state (an approximate Gibbs state),
after the trapping potential has been switched off. We start by introducing
our set-up.

We will be interested in quantum states with a fluctuating particle number.
These states are naturally defined on the bosonic Fock space

$\mathscr{F}(L^{2}(\mathbb{R}^{3}))=\bigoplus_{n=0}^{\infty}L^{2}_{\mathrm{sym}}(\mathbb{R}^{3n})\,,$
(1.1)

where $L^{2}_{\mathrm{sym}}(\mathbb{R}^{3n})$ denotes the closed linear
subspace of $L^{2}(\mathbb{R}^{3n})$ consisting of those functions
$\Psi(x_{1},...,x_{n})$ that are invariant under any permutation of the
coordinates $x_{1},...,x_{N}\in\mathbb{R}^{3}$. As usual, we have
$L^{2}_{\mathrm{sym}}(\mathbb{R}^{0})=\mathbb{C}$. By $a_{x}^{*}$ and $a_{x}$
we denote the creation and annihilation operators (actually operator-valued
distributions) on $\mathscr{F}$ that create and annihilate a particle at point
$x\in\mathbb{R}^{3}$, respectively. They satisfy the canonical commutation
relations (CCR)

$[a_{x},a_{y}^{*}]=\delta(x-y),\quad\quad[a_{x},a_{y}]=0=[a_{x}^{*},a_{y}^{*}]\,,$

where $\delta(x)$ denotes the delta distribution with unit mass at $x=0$.
Quantum states on the bosonic Fock space with an expected number of particles
equal to $N$ are elements of the set

$\mathcal{S}_{N}=\\{\Gamma\in\mathcal{B}(\mathscr{F})\ |\ \Gamma\geq
0,\operatorname{Tr}\Gamma=1,\operatorname{Tr}[\mathcal{N}\Gamma]=N\\}\,,$

where $\mathcal{N}=\int_{\mathbb{R}^{3}}a_{x}^{*}a_{x}\,\mathrm{d}x$ denotes
the particle number operator on $\mathscr{F}$. The one-particle reduced
density matrix (1-pdm) $\gamma$ of a state $\Gamma\in\mathcal{S}_{N}$ is
defined via its integral kernel by

$\gamma(x,y)=\operatorname{Tr}[a_{y}^{*}a_{x}\Gamma]$

and defines a nonnegative operator on $L^{2}(\mathbb{R}^{3})$ with
$\operatorname{Tr}\gamma=N$.

The time evolution we are interested in is governed by the Heisenberg equation

$\mathrm{i}\partial_{t}\Gamma_{t}=[\mathcal{H}_{N},\Gamma_{t}]\,,$ (1.2)

where

$\mathcal{H}_{N}=\int_{\mathbb{R}^{3}}\nabla a_{x}^{*}\nabla
a_{x}\,\mathrm{d}x+\frac{1}{2N}\int_{\mathbb{R}^{6}}v(x-y)a^{*}_{x}a^{*}_{y}a_{y}a_{x}\,\mathrm{d}x\,\mathrm{d}y$
(1.3)

and $[A,B]=AB-BA$ denotes the commutator of two linear operators $A$ and $B$.
The prefactor $N^{-1}$ in front of $v$ guarantees that, for particles in the
condensate, the kinetic and the potential energy are comparable. Any solution
to (1.2) with initial condition $\Gamma_{0}\in\mathcal{S}_{N}$ can be written
in terms of the strongly continuous unitary group
$e^{-\mathrm{i}\mathcal{H}_{N}t}$ as

$\Gamma_{t}=e^{-\mathrm{i}\mathcal{H}_{N}t}\Gamma_{0}e^{\mathrm{i}\mathcal{H}_{N}t}\,.$
(1.4)

The fact that $\mathcal{H}_{N}$ commutes with $\mathcal{N}$ guarantees
$\operatorname{Tr}[\mathcal{N}\Gamma_{t}]=\operatorname{Tr}[\mathcal{N}\Gamma_{0}]$.
As explained in Section 1.1, we are interested in the solution (1.4) for
initial data describing equilibrium states of trapped gases at positive
temperature. Let us now recall some known properties of such states.

###  1.4 Equilibrium states of trapped Bose gases

Equilibrium states of trapped Bose gases can be defined as minimizers, in
$\mathcal{S}_{N}$, of the Gibbs free energy functional

$\mathcal{F}(\Gamma)=\operatorname{Tr}[\mathcal{H}_{N}^{\mathrm{trap}}\Gamma]-TS(\Gamma)\,,\quad\text{
with the von Neumann entropy }\quad
S(\Gamma)=-\operatorname{Tr}[\Gamma\ln(\Gamma)]\,.$ (1.5)

Here

$\mathcal{H}_{N}^{\mathrm{trap}}=\mathcal{H}_{N}+\int_{\mathbb{R}^{3}}w(x)a_{x}^{*}a_{x}\,\mathrm{d}x=\int_{\mathbb{R}^{3}}a_{x}^{*}\left(-\Delta+w(x)\right)a_{x}\,\mathrm{d}x+\frac{1}{2N}\int_{\mathbb{R}^{6}}v(x-y)a^{*}_{x}a^{*}_{y}a_{y}a_{x}\,\mathrm{d}x\,\mathrm{d}y$
(1.6)

denotes the Hamiltonian of the trapped system and $T>0$ is the temperature.
The trapping potential $w(x)$ satisfies $w(x)\to\infty$ for $|x|\to\infty$.
The minimum of (1.5) is attained for the Gibbs state

$\Gamma_{\mathrm{G}}=\frac{\exp(-(\mathcal{H}_{N}^{\mathrm{trap}}-\mu\mathcal{N})/T)}{\operatorname{Tr}[\exp(-(\mathcal{H}_{N}^{\mathrm{trap}}-\mu\mathcal{N})/T)]}\,,$

where the chemical potential $\mu\in\mathbb{R}$ is chosen (depending on $N$)
so that $\operatorname{Tr}[\mathcal{N}\Gamma]=N$. This leads to the (grand
canonical) free energy

$F(T,N)=\min_{\Gamma\in\mathcal{S}_{N}}\mathcal{F}(\Gamma)=-T\ln\operatorname{Tr}[\exp(-(\mathcal{H}_{N}^{\mathrm{trap}}-\mu\mathcal{N})/T)]+\mu
N\,.$ (1.7)

Before we discuss what is known about approximate minimizers of the Gibbs free
energy functional $\mathcal{F}$, let us briefly recall some well-known facts
about the ideal (non-interacting) trapped Bose gas, that is, the system
described by the Hamiltonian in (1.6) with $v=0$. In the following we assume
that the trapping potential $w$ satisfies

$\lim_{|x|\to\infty}\frac{w(x)}{L|x|^{s}}=1$ (1.8)

for two constants $s>0$ and $L>0$. The one-particle density matrix
$\gamma^{\mathrm{id}}$ of the Gibbs state $\Gamma_{\mathrm{G}}$ of the ideal
gas equals

$\gamma^{\mathrm{id}}=\frac{1}{\exp((-\Delta+w-\mu_{N})/T)-1}\,.$ (1.9)

The chemical potential $\mu_{N}$ satisfies $\mu_{N}<e$, where $e$ denotes the
lowest eigenvalue of $-\Delta+w$, and is chosen s.t.
$\operatorname{Tr}\gamma^{\mathrm{id}}=N$ holds. By definition, the system
displays BEC if and only if the largest eigenvalue $N_{0}(T,N)$ of
$\gamma^{\mathrm{id}}$ is of order $N$ in the limit $N\to\ \infty$, that is,
if and only if

$\lim_{N\to\infty}\frac{N_{0}(T,N)}{N}=g>0\,,\quad\text{ where }\quad
N_{0}(T,N)=\frac{1}{\exp((e-\mu_{N})/T)-1}\,.$ (1.10)

The trapped ideal Bose gas is well-known to exhibit a BEC phase transition
with critical temperature given by

$T_{\mathrm{c}}(s)=\frac{N^{1/\alpha}}{(\kappa\alpha\Gamma(\alpha)\zeta(\alpha))^{1/{\alpha}}}\,,\quad\text{
where
}\quad\kappa=\frac{2L^{-3/s}}{3\pi}\int_{0}^{1}(1-x)^{3/2}x^{2}\,\mathrm{d}x\quad\text{
and }\quad\alpha=\frac{6+3s}{2s}\,.$ (1.11)

Here $\Gamma(\alpha)$ denotes the Gamma function and $\zeta(\alpha)$ is the
Riemann zeta function. More precisely, let us define
$t_{\mathrm{c}}(s)=[(\kappa\alpha\Gamma(\alpha)\zeta(\alpha))^{1/{\alpha}}]^{-1}$,
that is, $t_{\mathrm{c}}(s)=N^{-1/\alpha}T_{\mathrm{c}}(s)$. If $T=\lambda
N^{1/\alpha}$ with $\lambda\geq 0$, then we have
$g=[1-(\lambda/t_{\mathrm{c}}(s))^{\alpha}]_{+}$ with $[x]_{+}=\max\\{x,0\\}$.
This statement follows from the definition of $N_{0}(T,N)$ and the eigenvalue
asymptotics of $-\Delta+w$, see e.g. [64, Theorem XIII.81]111The theorem in
the reference is stated for $s>1$ but the proof applies for all $s>0$..

An important property of the ideal gas at temperature $T=\lambda N^{1/\alpha}$
with fixed $\lambda\in(0,t_{\mathrm{c}}(s))$ is that its density
$\gamma^{\mathrm{id}}(x,x)$ displays a two-scale structure. The condensed
particles are described by the ground state of $-\Delta+w$. Accordingly, their
density varies on a length scale of order one, independent of $N$. Excited
particles (the thermal cloud) are characterized by the much larger length
scale $R\gg 1$ because their energy per particle is of the order $T$. The
length $R$ can be obtain by equating $T=R^{s}$, that is, $R\sim N^{2/(6+3s)}$.
The density $\varrho_{\mathrm{th}}$ of the thermal cloud is proportional to
$\varrho_{\mathrm{th}}\sim N/R^{3}\sim T^{3/2}\sim N^{3s/(6+3s)}\ll N$, that
is, the thermal cloud is much more dilute than the condensate.

As recently shown in [24], this two scale structure is preserved if we add a
weak interaction to the Hamiltonian, as in (1.6). To make this statement
precise, we introduce for a function $\phi$ in the form domain of $-\Delta+w$
and $g\in[0,1]$ the Hartree energy functional

$\mathcal{E}^{\mathrm{H}}(\phi)=\langle\phi,(-\Delta+w)\phi\rangle+\frac{g}{2}\int_{\mathbb{R}^{6}}|\phi(x)|^{2}v(x-y)|\phi(y)|^{2}\,\mathrm{d}x\,\mathrm{d}y\,.$
(1.12)

In the following we assume that $v$ is such that the functional
$\mathcal{E}^{\mathrm{H}}$ is bounded from below and admits a unique minimizer
in the set of functions with $L^{2}(\mathbb{R}^{3})$-norm equal to one. We
denote the minimum and the minimizer of $\mathcal{E}^{\mathrm{H}}$ by
$E^{\mathrm{H}}(g)$ and $\phi^{\mathrm{H}}$, respectively. Using Theorem 1.3
and Lemma 7.1 in [24], it can be shown222The result in the reference is stated
for the special choice $w(x)=x^{2}$. However, the techniques that have been
used to prove this result also apply to the more general trapping potentials
we introduced in the discussion of the ideal gas. that for $0\leq T\lesssim
T_{\mathrm{c}}(s)$, with $T_{\mathrm{c}}(s)$ defined in (1.11) as the critical
temperature of the ideal trapped gas, the free energy in (1.7) satisfies

$\lim_{N\to\infty}N^{-1}\left|F(T,N)-F_{0}(T,N)-N_{0}(T,N)E^{\mathrm{H}}(g)\right|=0\,.$
(1.13)

Here $N_{0}(N,T)$ and $g$ are chosen as in (1.10), and $F_{0}(T,N)$ denotes
the free energy in (1.7) with $v=0$. The effect of the interaction can be
observed, to leading order in $N$, only in the condensate because the thermal
cloud is much more dilute. Another consequence of the results in [24] is that
for any approximately minimising sequence $\Gamma_{N}\in\mathcal{S}_{N}$
satisfying

$\lim_{N\to\infty}N^{-1}\left|\mathcal{F}(\Gamma_{N})-F_{0}(T,N)-N_{0}(T,N)E^{\mathrm{H}}(g)\right|=0\,,$
(1.14)

the corresponding 1-pdm $\gamma_{N}$ is such that

$\lim_{N\to\infty}N^{-1}\left\|\gamma_{N}-\gamma^{\mathrm{id}}Q-N_{0}(T,N)|\phi^{\mathrm{H}}\rangle\langle\phi^{\mathrm{H}}|\
\right\|_{\mathcal{L}^{1}}=0$ (1.15)

holds, where $Q=1-|\psi_{0}\rangle\langle\psi_{0}|$ with the ground state
$\psi_{0}$ of $-\Delta+w$. That is, to leading order in $N$, $\gamma_{N}$ is
the 1-pdm of the ideal gas, where the eigenfunction of its largest eigenvalue
(the condensate wave function) has been replaced by the minimizer of the
Hartree energy functional. This, in particular, implies that the interacting
systems displays a BEC phase transition with critical temperature given by the
one of the ideal gas to leading order. It also indicates that the density of
the Gibbs state $\Gamma_{\mathrm{G}}$ has the same two-scale structure as that
of the ideal trapped gas.

###  1.5 Construction of the initial data

In this section we define our class of initial states. Our definition is
motivated by the discussion in Section 1.4.

Let $\phi\in L^{2}(\mathbb{R}^{3})$ and let
$\gamma\in\mathcal{B}(L^{2}(\mathbb{R}^{3}))$ be a positive trace class
operator. We think of $\phi$ as a condensate wave function and of $\gamma$ as
a 1-pdm describing excitations of the condensate. The particle number
$n(\phi,\gamma)$ of the pair $(\phi,\gamma)$ is defined by

$n(\phi,\gamma)=\int_{\mathbb{R}^{3}}|\phi(x)|^{2}\,\mathrm{d}x+\operatorname{Tr}\gamma\,.$
(1.16)

Let $G$ be the unique quasi-free state on $\mathscr{F}$ that satisfies
$[G,\mathcal{N}]=0$ and has $\gamma$ as 1-pdm. We also define the Weyl
transformation

$W(\phi)=\exp(a^{*}(\phi)-a(\phi))\,.$ (1.17)

The above definitions allow us to associate to each pair $(\phi,\gamma)$ the
state

$\Gamma(\phi,\gamma)=W(\phi)GW(\phi)^{*}$ (1.18)

on the bosonic Fock space. The kernel of its 1-pdm is given by

$\operatorname{Tr}[a_{y}^{*}a_{x}\Gamma(\phi,\gamma)]=\phi(x)\overline{\phi(y)}+\gamma(x,y)\,,$
(1.19)

which, in particular, implies
$\operatorname{Tr}[\mathcal{N}\Gamma(\phi,\gamma)]=n(\phi,\gamma)$. The state
$\Gamma(\phi,\gamma)$ is still a quasi-free state, in the sense that higher
order correlations can be computed through the Wick theorem, summing over all
possible partitions in groups of one or two creation and annihilation
operators. Choosing $\phi=\sqrt{N_{0}(T,N)}\phi^{\mathrm{H}}$ and
$\gamma=\gamma^{\mathrm{id}}Q$, it is simple to check that
$\Gamma(\phi,\gamma)$ provides a good approximation for the free energy of the
system in the sense of (1.14). In our main theorem, we will describe the
evolution in (1.4) for initial data given by suitable perturbations of a
quasi-free state of the form $\Gamma(\phi,\gamma)$, under appropriate
assumptions on $(\phi,\gamma)$.

To this end, we first use the Araki–Woods representation, see [3, 23], to
recast a mixed state $\Gamma\in\mathcal{S}_{N}$ as a vector $\Psi$ on the
doubled Fock space $\mathscr{F}\otimes\mathscr{F}$. This will allow us to
apply the formalism of Bogoliubov transformations and fluctuation dynamics.
Note that a similar strategy was used in [7] to study the mean-field evolution
of fermionic mixed states. The vector $\Psi$ is conveniently defined in terms
of the spectral decomposition
$\Gamma=\sum_{\alpha=1}^{\infty}\lambda_{\alpha}|\psi_{\alpha}\rangle\langle\psi_{\alpha}|$
as

$\Psi=\sum_{\alpha=1}^{\infty}\lambda_{\alpha}^{1/2}\psi_{\alpha}\otimes\overline{\psi_{\alpha}},$
(1.20)

where we define the complex conjugate of $\psi_{\alpha}$ in the position-space
representation. It satisfies

$\operatorname{Tr}[A\Gamma]=\langle\Psi,\left(A\otimes\mathds{1}\right)\Psi\rangle$
(1.21)

for $A\in\mathcal{B}(\mathscr{F})$. It should be noted that this
representation is far from being unique. When we replace the vectors
$\psi_{\alpha}$ in the second tensor factor in (1.20) by elements of another
orthonormal basis of $\mathscr{F}$, we obtain an equivalent representation of
$\Gamma(\phi,\gamma)$ on $\mathscr{F}\otimes\mathscr{F}$ in the sense that
(1.21) continues to hold.

For the next step in our construction we need the unitary equivalence

$\mathcal{U}\mathscr{F}(L^{2}(\mathbb{R}^{3}))\otimes\mathscr{F}(L^{2}(\mathbb{R}^{3}))=\mathscr{F}(L^{2}(\mathbb{R}^{3})\oplus
L^{2}(\mathbb{R}^{3}))$

defined by the relations $\mathcal{U}\Omega\otimes\Omega=\widetilde{\Omega}$,
where $\widetilde{\Omega}$ denotes the vacuum vector of
$\mathscr{F}(L^{2}(\mathbb{R}^{3})\oplus L^{2}(\mathbb{R}^{3}))$,

$\mathcal{U}\big{(}a(f)\otimes\mathds{1}\big{)}\mathcal{U}^{*}=a(f\oplus
0)\eqqcolon
a_{\mathrm{\ell}}(f),\quad\text{and}\quad\mathcal{U}\big{(}\mathds{1}\otimes
a(f)\big{)}\mathcal{U}^{*}=a(0\oplus f)\eqqcolon a_{r}(f).$ (1.22)

Thus, every state $\Gamma\in\mathcal{S}_{N}$ can be represented by a
normalized vector
$\widetilde{\Psi}=\mathcal{U}\Psi\in\mathscr{F}(L^{2}(\mathbb{R}^{3})\oplus
L^{2}(\mathbb{R}^{3}))$ satisfying

$\operatorname{Tr}[A\Gamma]=\langle\widetilde{\Psi}\mathcal{U}(A\otimes\mathds{1})\mathcal{U}^{*}\widetilde{\Psi}\rangle$

If $A$ is expressed in terms of creation and annihilation operators $a^{*},a$,
then $\mathcal{U}(A\otimes 1)\mathcal{U}^{*}$ is obtained from $A$ by
replacing all operators $a^{*}$ with $a^{*}_{\ell}$ and all operators $a$ with
$a_{\ell}$. Let us now consider, in particular, the state (1.18). An important
observation is that the quasi-free state $G$ can be described, on
$\mathcal{F}(L^{2}(\mathbb{R}^{3})\oplus L^{2}(\mathbb{R}^{3}))$, by the
vector $\widetilde{\Psi}_{G}=\mathcal{T}(\gamma)\widetilde{\Omega}$, where

$\mathcal{T}(\gamma)=\exp\left(\int_{\mathbb{R}^{6}}k_{\gamma}(x,y)a^{*}_{\mathrm{\ell},x}a^{*}_{r,y}\,\mathrm{d}x\,\mathrm{d}y-\text{h.c.}\right)$
(1.23)

and $k_{\gamma}(x,y)$ denotes the integral kernel of the operator
$k_{\gamma}=\mathrm{arcsinh}(\sqrt{\gamma})$. We observe that
$\mathcal{T}(\gamma)$ is a Bogoliubov transformation and satisfies the
relations

$\mathcal{T}^{*}(\gamma)a_{\ell}(f)\mathcal{T}(\gamma)=a_{\ell}(uf)+a_{r}^{*}(\overline{vf})\,,\hskip
42.67912pt\mathcal{T}^{*}(\gamma)a_{r}(f)\mathcal{T}(\gamma)=a_{r}(uf)+a_{\ell}^{*}(\overline{vf})\,,$
(1.24)

where $u=\sqrt{1+\gamma}$ and $v=\sqrt{\gamma}$. With (1.24), it is clear that
$\widetilde{\Psi}_{G}$ is quasi-free (it satisfies the Wick theorem) and that
it has the same 1-pdm as $G$, in the sense that

$\langle\widetilde{\Psi}_{G},a_{\ell,y}^{*}a_{\ell,x}\widetilde{\Psi}_{G}\rangle=\gamma(x,y)\,.$

Recalling (1.18), we find that the quasi-free state $\Gamma(\phi,\gamma)$, now
with a condensate described by the wave function $\phi$, can be described, on
$\mathscr{F}(L^{2}(\mathbb{R}^{3})\oplus L^{2}(\mathbb{R}^{3}))$ by the vector
$\widetilde{\Psi}(\phi,\gamma)=\mathcal{W}(\phi)\mathcal{T}(\gamma)\widetilde{\Omega}$
(sometimes called the purification of $\Gamma(\phi,\gamma)$), in the sense
that

$\operatorname{Tr}[A\Gamma(\phi,\gamma)]=\langle\widetilde{\Omega},\mathcal{T}^{*}(\gamma)\mathcal{W}^{*}(\phi)\mathcal{U}\left(A\otimes\mathds{1}\right)\mathcal{U}^{*}\mathcal{W}(\phi)\mathcal{T}(\gamma)\widetilde{\Omega}\rangle\,.$
(1.25)

Here, we introduced, on $\mathscr{F}(L^{2}(\mathbb{R}^{3})\oplus
L^{2}(\mathbb{R}^{3}))$, the Weyl operator

$\mathcal{W}(\phi)=\exp\left(a^{*}_{\mathrm{\ell}}(\phi)+a^{*}_{r}(\overline{\phi})-\mathrm{h.c.}\right)$
(1.26)

satisfying the relations

$\mathcal{W}^{*}(\phi)a_{\ell}(f)\mathcal{W}(\phi)=a_{\ell}(f)+\langle
f,\phi\rangle\,,\hskip
42.67912pt\mathcal{W}^{*}(\phi)a_{r}(f)\mathcal{W}(\phi)=a_{r}(f)+\langle
f,\overline{\phi}\rangle\,.$ (1.27)

The validity of (1.25) can be checked similarly as the comparable equality for
$G$, using (1.24) and (1.27).

In our main theorem, we will study the dynamics of initial states
$\Gamma_{\xi}(\phi,\gamma)$, represented, on
$\mathscr{F}(L^{2}(\mathbb{R}^{3})\oplus L^{2}(\mathbb{R}^{3}))$, by the
vector $\mathcal{W}(\phi)\mathcal{T}(\gamma)\xi$, for a normalized
$\xi\in\mathscr{F}(L^{2}(\mathbb{R}^{3})\oplus L^{2}(\mathbb{R}^{3}))$. That
is, the state $\Gamma_{\xi}(\phi,\gamma)$ is defined so that, for every
observable $A\in\mathcal{B}(\mathscr{F})$,

$\operatorname{Tr}[A\Gamma_{\xi}(\phi,\gamma)]=\langle\xi,\mathcal{T}^{*}(\gamma)\mathcal{W}^{*}(\phi)\mathcal{U}\left(A\otimes\mathds{1}\right)\mathcal{U}^{*}\mathcal{W}(\phi)\mathcal{T}(\gamma)\xi\rangle\,.$
(1.28)

Our assumptions on $\phi,\gamma,\xi$ will make sure that
$\Gamma_{\xi}(\phi,\gamma)$ is a perturbation of the quasi-free state
$\Gamma(\phi,\gamma)$ (in particular, we will require the expectation of the
number of particles operator and sufficiently many of its moments to stay
bounded, in the state $\xi$, uniformly in $N$). From (1.4), we conclude that
the evolution
$\Gamma_{\xi,t}(\phi,\gamma)=e^{-\mathrm{i}\mathcal{H}_{N}t}\Gamma_{\xi}(\phi,\gamma)e^{\mathrm{i}\mathcal{H}_{N}t}$
of the initial data $\Gamma_{\xi}(\phi,\gamma)$ is such that

$\operatorname{Tr}[A\Gamma_{\xi,t}(\phi,\gamma)]=\langle\xi,\mathcal{T}^{*}(\gamma)\mathcal{W}^{*}(\phi)\exp(\mathrm{i}\mathcal{L}_{N}t)\mathcal{U}\left(A\otimes\mathds{1}\right)\mathcal{U}^{*}\exp(-\mathrm{i}\mathcal{L}_{N}t)\mathcal{W}(\phi)\mathcal{T}(\gamma)\xi\rangle\,,$
(1.29)

with the Liouvillian

$\mathcal{L}_{N}=\mathcal{H}_{N,\mathrm{\ell}}-\mathcal{H}_{N,r}\,.$ (1.30)

The operators $\mathcal{H}_{N,\mathrm{\ell}}$ and $\mathcal{H}_{N,r}$ are
defined as the original Hamiltonian $\mathcal{H}_{N}$ in (1.3) with
$a_{x}^{*},a_{x}$ replaced by $a_{\mathrm{\ell},x}^{*},a_{\mathrm{\ell},x}$
and $a_{r,x}^{*},a_{r,x}$, respectively. In other words, on
$\mathscr{F}(L^{2}(\mathbb{R}^{3})\oplus L^{2}(\mathbb{R}^{3}))$, we are
interested in the evolution $e^{-\mathrm{i}\
\mathcal{L}_{N}t}\mathcal{W}(\phi)\mathcal{T}(\gamma)\xi$ generated by
$\mathcal{L}_{N}$ on the initial data
$\mathcal{W}(\phi)\mathcal{T}(\gamma)\xi$.

Motivated by the properties of approximate equilibrium states of trapped Bose
gases with a trapping potential satisfying (1.8) for some $s>0$, we will make
the following assumptions on the pair $(\phi,\gamma)$.

######  Assumptions 1.1.

We assume that the pair $(\phi,\gamma)$ is s.t. $\phi\in
H^{3}(\mathbb{R}^{3})$ and
$\operatorname{Tr}[(1-\Delta)^{3/2}\gamma(1-\Delta)^{3/2}]<+\infty$. Moreover,
$n(\phi,\gamma)=N$ and the following holds for $s>0$:

(A) The condensate wave function $\phi$ can be written as
$\phi=c(N)\widetilde{\phi}$, with a fixed normalized $\widetilde{\phi}\in
L^{2}(\mathbb{R}^{3})$. The $N$-dependent constant $c(N)$ determines the
expected number of particles in the condensate.

(B) The 1-pdm $\gamma$ obeys

$\int_{\mathbb{R}^{6}}|\hat{\gamma}(p,q)|\,\mathrm{d}p\,\mathrm{d}q\lesssim
T_{\mathrm{c}}^{3/2}(s)$

with $T_{\mathrm{c}}(s)$ as in (1.11), and where $\hat{\gamma}(p,q)$ denotes
the integral kernel of $\gamma$ in Fourier space.

(C) The operator norm of $\gamma$ satisfies

$\|\gamma\|\lesssim T_{\mathrm{c}}(s).$

Part (A) in the above assumption allows us to choose
$\widetilde{\phi}=\phi^{H}$ minimising (1.12) with $g$ in (1.10), and
$c(N)=\sqrt{N_{0}(N,T)}$, as defined in (1.10). Point (B) implies that the
density satisfies $|\gamma(x,x)|\lesssim T_{\mathrm{c}}^{3/2}(s)$, uniformly
in $x\in\mathbb{R}^{3}$. In fact, it also implies (and this will be important
for us) a bound (uniform in $t$) for the density of the free evolution
$e^{\mathrm{i}\Delta t}\gamma e^{-\mathrm{i}\Delta t}$ and more generally

$\sup_{x,y\in\mathbb{R}^{3}}|(e^{\mathrm{i}\Delta t}\gamma
e^{-\mathrm{i}\Delta
t})(x,y)|\leq\int_{\mathbb{R}^{6}}|\hat{\gamma}(p,q)|\,\mathrm{d}p\,\mathrm{d}q\lesssim
T_{\mathrm{c}}^{3/2}(s)\,.$

As we will see in Section 2, point (B) is satisfied by $\gamma^{\mathrm{id}}Q$
with the 1-pdm $\gamma^{\mathrm{id}}$ of the ideal gas in (1.9) with the
choice $w(x)=|x|^{s}$ and with $Q=1-|\psi_{0}\rangle\langle\psi_{0}|$, where
$\psi_{0}$ denotes the ground state of $-\Delta+|x|^{s}$. Point (C) in
Assumption 1.1 is satisfied by $\gamma^{\mathrm{id}}Q$, too.

We conclude this section with a brief discussion of the relevant scales
related to our time evolution. In the whole discussion we assume that
$T=\kappa T_{\mathrm{c}}(s)$ with $\kappa\in(0,1)$. This in particular implies
that we have order $N$ particles in the condensate and order $N$ particles in
the thermal cloud. We are interested in following the dynamics of the
condensate, and therefore consider times of order one. This is motivated by
the fact that the energy per particle in the condensate and the related length
scale are both of order one, independent of $N$, see Section 1.4. In contrast,
the energy per particle in the thermal cloud is proportional to $T$ and the
related velocity is of the order $T^{1/2}\sim N^{s/(6+3s)}\gg 1$. Thus,
particles in the thermal cloud move much faster than particles in the
condensate. For particles in the condensate, the $N^{-1}$ factor appearing in
(1.3) in front of the interaction produces a mean-field potential and leads to
the Hartree equation. In contrast, the typical length scale $R$ of the thermal
cloud scales as $N^{2/(6+3s)}\gg 1$ (see Section 1.4); hence the thermal cloud
is much more dilute and we can expect that, to leading order in $N$, it moves
according to the free evolution. Our main result confirms this heuristic
picture.

###  1.6 Main result

Our main result is captured in the following theorem.

######  Theorem 1.2.

We assume that $v\in L^{1}(\mathbb{R}^{3})\cap W^{1,p}(\mathbb{R}^{3})$ with
$p>3$ satisfies $v(x)=v(-x)$ for a.e. $x\in\mathbb{R}^{3}$ and that
$\hat{v}\in L^{1}(\mathbb{R}^{3})$, where $\hat{v}$ denotes the Fourier
transform of v. Moreover, we assume that the pair $(\phi,\gamma)$ satisfies
Assumption 1.1 with $0<s\leq 3/2$ and that the fluctuation vector $\xi$ is
s.t. $\langle\xi,(\mathcal{N}_{\ell}+\mathcal{N}_{r})^{44}\xi\rangle\lesssim
1$ holds with $\mathcal{N}_{\ell/r}=\int
a^{*}_{\ell/r,x}a_{\ell/r,x}\,\mathrm{d}x$. Then there exists a constant $c>0$
s.t. the 1-pdm $\gamma_{\xi,t}$ of the state $\Gamma_{\xi,t}(\phi,\gamma)$
defined in (1.29) satisfies

$\left\|\gamma_{\xi,t}-e^{\mathrm{i}\Delta t}\gamma e^{-\mathrm{i}\Delta
t}-|\phi_{t}\rangle\langle\phi_{t}|\
\right\|_{\mathcal{L}^{1}}\lesssim\sqrt{N}T_{\mathrm{c}}^{3/4}(s)\exp(c\exp(c\exp(ct))),$
(1.31)

where $\|\cdot\|_{\mathcal{L}^{1}}$ denotes the trace norm. The function
$\phi_{t}$ is the solution to the time-dependent Hartree equation

$\mathrm{i}\partial_{t}\phi_{t}(x)=\left(-\Delta+N^{-1}v\ast|\phi_{t}(x)|^{2}\right)\phi_{t}(x)\quad\text{
with initial datum }\quad\phi_{0}(x)=\phi(x).$ (1.32)

We have the following remarks concerning Theorem 1.2.

######  Remark 1.3.

  1. 1.

The r.h.s. of (1.31) scales for $0<s\leq 3/2$ as $N^{(1+s)/(2+s)}\leq
N^{5/7}\ll N$, which should be compared to
$\operatorname{Tr}\gamma_{\xi,t}=N=\operatorname{Tr}[e^{\mathrm{i}\Delta
t}\gamma e^{-\mathrm{i}\Delta t}+|\phi_{t}\rangle\langle\phi_{t}|]$. That is,
the 1-pdm $e^{\mathrm{i}\Delta t}\gamma e^{-\mathrm{i}\Delta
t}+|\phi_{t}\rangle\langle\phi_{t}|$ yields a good approximation for
$\gamma_{\xi,t}$ for fixed $t>0$ and large $N$.

  2. 2.

The time-dependence of the r.h.s. of (1.31) can be replaced by
$C_{\varepsilon}\exp(c\exp(c(t+t^{3+\varepsilon})))$ for $\varepsilon>0$
provided we have either $v\geq 0$ or $\hat{v}\geq 0$. More details can be
found in Remark 3.6 in Section 3.2 below.

  3. 3.

We can obtain a better approximation for $\gamma_{\xi,t}$ by letting
$\omega=|\phi\rangle\langle\phi|+\gamma$ and considering the solution to the
time-dependent Hartree equation

$\mathrm{i}\partial_{t}\omega_{t}=[-\Delta+N^{-1}v*\varrho_{\omega_{t}}(x),\omega_{t}]$
(1.33)

with $\omega_{t=0}=\omega$ and $\varrho_{\omega_{t}}(x)=\omega_{t}(x,x)$.
Then, proceeding as in our proof of Theorem 1.2, we could show that

$\|\gamma_{\xi,t}-\omega_{t}\|_{\mathcal{L}^{1}}\lesssim\sqrt{NT_{\mathrm{c}}(s)}\exp(c\exp(c\exp(ct)))\,.$

In this case the dependence on $N$ on the r.h.s. is optimal, already for
$t=0$. In fact, taking
$\xi=\mathcal{U}(1+a_{r}^{*}(\overline{\phi}/\|\phi\|))\Omega/\sqrt{2}$, with
$\mathcal{U}$ defined in (1.22), a short computation shows that

$\operatorname{Tr}[\mathcal{N}\Gamma_{\xi}(\phi,\gamma)]-\operatorname{Tr}[\mathcal{N}\Gamma_{\Omega}(\phi,\gamma)]\geq\frac{\langle\phi,\sqrt{\gamma}\phi\rangle}{\|\phi\|}\,.$
(1.34)

Assume $\|\gamma\|\sim T_{c}(s)$ and let $\psi$ denote the normalized
eigenfunction of $\gamma$ associated with its largest eigenvalue. Choosing
$\phi=c(N)\psi$ so that $\|\phi\|\sim\sqrt{N}$, the r.h.s. of (1.34) is of
order $\sqrt{NT_{c}(s)}$.

  4. 4.

The assumption for $\xi$ appearing in Theorem is needed for $s=3/2$. For
$s<3/2$, less moments are sufficient.

  5. 5.

The theorem holds, in particular, for the following choice of $(\phi,\gamma)$:
Let $h=-\Delta+|x|^{s}$ with $0<s\leq 3/2$. The 1-pdm $\gamma$ is defined by
$\gamma=\gamma^{\mathrm{id}}Q$ with $\gamma^{\mathrm{id}}$ in (1.9), where the
projection $Q$ removes the largest eigenvalue of $\gamma^{\mathrm{id}}$ (the
condensate). Moreover, the temperature $T$ in the definition of
$\gamma^{\mathrm{id}}$ obeys $T\lesssim T_{\mathrm{c}}(s)$. The condensate
wave function $\phi$ is given by
$\sqrt{N-\operatorname{Tr}[\gamma^{\mathrm{id}}Q]}$ times a minimizer of the
Hartree energy functional in (1.12). In Section 2 we show that the pair
$(\phi,\gamma)$ satisfies all assumptions of Theorem 1.2.

  6. 6.

The condition $0<s\leq 3/2$ is needed to consider initial data like those
described in Remark $5$, for $T=\kappa T_{c}(s)$, $\kappa\in(0,1)$. This
restriction is required, because, if $s>3/2$, the largest eigenvalue of
$\gamma^{\mathrm{id}}Q$ grows so fast with $N$ that our analysis breaks down.
It is possible to consider traps with $s>3/2$, in particular the harmonic trap
$w(x)=x^{2}$, if we restrict our attention to temperatures $T\lesssim
T_{c}(3/2)\ll T_{c}(2)$. Moreover, our analysis also applies above the
critical temperature, where $\phi=0$ and $\|\gamma\|\sim 1$; in this case, we
can take any $s\in(0;\infty)$ and we only need weaker assumptions on $\xi$
(and we get better bounds).

###  1.7 Proof strategy and organization of the article

For the convenience of the reader we give in this section a short summary of
the main steps leading to a proof of Theorem 1.2.

In Section 2 we show that the initial data in Remark 1.3.5, consisting of the
minimizer of the Hartree energy functional in (1.12) and the 1-pdm of the
ideal gas trapped by the potential $w(x)=|x|^{s}$ with $s\in(0,3/2]$,
satisfies Assumption 1.1.

Our many-body analysis is based on the definition of fluctuation dynamics
around the dynamics generated by the Hartree–Fock–Bogoliubov (HFB) equations,
see e.g. [4, 10]. We do not give a derivation of these equations from quantum
mechanics but rather use them as a technical tool. In Section 3 we introduce
the HFB equations, we recall their well-posedness, as established in [4], and
we show how they can be approximated by the much simpler effective dynamics
appearing in Theorem 1.2, if the initial data satisfies Assumption 1.1. We
also prove a bound for solutions of the HFB equations that guarantees the
diluteness of the thermal cloud during the time evolution and is a crucial
ingredient for our many-body analysis.

Section 4 is devoted to the construction of the fluctation dynamics around the
HFB equations and to the proof of Theorem 1.2. We start by recalling some
basic facts about Weyl and Bogoliubov transformations in Section 4.1, and
afterwards define the fluctuation dynamics in Section 4.2. In Section 4.3 we
prove a bound for the trace-norm difference of $\gamma_{\xi,t}$ and the 1-pdm
related to the solution to the HFB equations in terms of the expected number
of excitations in the time-dependent fluctuation vector. A bound for this
quantity is stated without proof in Proposition 4.3 in Section 4.4. We end
Section 4 by showing how Proposition 4.3 implies Theorem 1.2.

In Section 5 we give a proof of Proposition 4.3 that is based on a Grönwall
argument for the expected number of particles in the fluctuation dynamics, and
we start our analysis in Section 5.1 with the computation of its generator. In
Section 5.2 we introduce a fluctuation dynamics with a cut-off in the particle
number and prove the equivalent of Proposition 4.3 for this dynamics. Section
5.3 is devoted to the proof of weak a-priori bounds for the original
fluctation dynamics without a cut-off, that are later used in Section 5.4 to
show that the two fluctuation dynamics are close in a suitable sense. Finally,
in Section 5.5 we use the results of the previous sections to give a proof of
Proposition 4.3, thus concluding the proof of Theorem 1.2.

##  2 Properties of the 1-pdm of the ideal gas and of Hartree minimizers

In this section we provide an example of a physically relevant initial pair
$(\phi,\gamma)$ satisfying Assumption 1.1. To this end, we choose $\gamma$ as
the 1-pdm of the ideal gas described by the Schrödinger operator
$h=-\Delta+|x|^{s}$ for some $s\in(0,2]$ and $\phi$ as the minimizer of the
Hartree functional (1.12), again with $h=-\Delta+|x|^{s}$ and with a
sufficiently regular interaction potential $v$ (Theorem 1.2 describes the
evolution of perturbations of the state $\Gamma(\phi,\gamma)$ in (1.18), if
$s\in(0;3/2]$). The next proposition shows the validity of Remark 1.3.5.

######  Proposition 2.1.

Let $h=-\Delta+|x|^{s}$, for some $s\in[0,2]$ and

$\gamma=\frac{1}{e^{\beta(h-\mu)}-1}Q$

with a chemical potential $\mu\in\mathbb{R}$ satisfying $h>\mu$ and with
$Q=1-|\psi_{0}\rangle\langle\psi_{0}|$ denoting the projection on the
orthogonal complement of the ground state $\psi_{0}$ of $h$. Then, for every
$\beta>0$ there exists a constant $C<\infty$ (depending on $\beta$), with

$\operatorname{Tr}\,(1-\Delta)^{3/2}\,\gamma(1-\Delta)^{3/2}\leq C\,.$ (2.1)

Moreover, denoting by $\hat{\gamma}(p,q)$ the integral kernel of $\gamma$ in
Fourier space, we find $\|\gamma\|\lesssim T_{\mathrm{c}}(s)$ and

$\int_{\mathbb{R}^{6}}|\hat{\gamma}(p,q)|\,\,\mathrm{d}p\,\mathrm{d}q\lesssim
T_{\mathrm{c}}^{3/2}(s)$ (2.2)

for all $\beta>0$ with $\beta^{-1}\lesssim T_{\mathrm{c}}(s)$.

Let $v\in W^{1,p}(\mathbb{R}^{3})$ for some $p>3$ and let $\phi\in
H^{1}(\mathbb{R}^{3})$ with $\|\phi\|=1$ solve the Hartree equation

$\Big{(}-\Delta+|x|^{s}+g(v*|\phi|^{2})(x)\Big{)}\phi(x)=\mu^{\mathrm{H}}\phi(x)$
(2.3)

with $g\in[0,1]$ in the sense of distributions, for some
$\mu^{H}\in\mathbb{R}$. Then we have
$\|\phi\|_{H^{3}}\lesssim(1+\|v\|_{\infty}+\mu^{\mathrm{H}})^{3/2}$.

In order to prove Prop. 2.1, we will make use of the following lemma.

######  Lemma 2.2.

Let $w\in W^{1,1}_{\text{loc}}(\mathbb{R}^{3})$ be real-valued, with
$w(x)\geq-C$, for some $C\geq 0$. Let $h_{w}=-\Delta+w$, denote $\kappa=1+C$,
and assume the operator inequality

$|\nabla w|^{2}\lesssim(\kappa+h_{w})\,.$ (2.4)

Then, we have

$(1-\Delta)^{2},w^{2}\lesssim(\kappa+h_{w})^{2},\qquad(1-\Delta)^{3}\lesssim(\kappa+h_{w})^{3}\,.$
(2.5)

Remark: In applications, we will use Lemma 2.2, with $w(x)=|x|^{s}$ and with
$w(x)=|x|^{s}+g(v*|\phi|^{2})(x)$, for $s\in[0,2]$.

###### Proof.

Pick $u\in\mathcal{C}_{\mathrm{c}}^{\infty}(\mathbb{R}^{3})$. A short
computation shows

$\langle u,(\kappa-\Delta+w)^{2}u\rangle=\langle
u,[(1-\Delta)^{2}+(C+w)^{2}]u\rangle+2\langle\nabla u,(C+w)\nabla
u\rangle+2\text{Re}\langle u,(\nabla w)\nabla u\rangle+2\langle
u,(C+w)u\rangle\,.$ (2.6)

With (2.4), we find

$|\langle u,(\nabla w)\nabla u\rangle|\leq\|(\nabla w)u\|\,\|\nabla
u\|\lesssim\langle u,(\kappa+h_{w})u\rangle\,.$

Estimating $(\kappa+h_{w})\leq\varepsilon(\kappa+h_{w})^{2}+C_{\varepsilon}$
and using $w\geq-C$ we conclude, from (2.6), that

$\langle u,(\kappa+h_{w})^{2}u\rangle\gtrsim\langle
u,\big{[}(1-\Delta)^{2}+(C+w)^{2}\big{]}u\rangle\,.$ (2.7)

Since $h_{w}$ is essentially self-adjoint on $C^{\infty}_{c}(\mathbb{R}^{3})$
(see for instance [63, Theorem X.28]), this proves the first bound in (2.5).
To show the second bound, consider

$(1-\Delta)^{3}=(1-\Delta)^{2}-\nabla\cdot(1-\Delta)^{2}\,\nabla\,.$

From (2.7), we have $(1-\Delta)^{2}\lesssim(\kappa+h_{w})^{2}$ and therefore
that

$\begin{split}-\nabla\cdot(1-\Delta)^{2}\,\nabla&\lesssim-\nabla\cdot(\kappa+h_{w})^{2}\nabla=-((h_{w}+\kappa)\nabla+[\nabla,h_{w}])\cdot(\nabla(h_{w}+\kappa)-[\nabla,h_{w}])\\\
&\lesssim(\kappa+h_{w})(-\Delta)(\kappa+h_{w})+|[\nabla,h_{w}]|^{2}\,.\end{split}$

With $-\Delta\leq\kappa+h_{w}$, $[\nabla,h_{w}]=\nabla w$ and the assumption
(2.4), we conclude that

$(1-\Delta)^{3}\lesssim(\kappa+h_{w})^{3}\,.$

∎

We will also need the following estimate on the Fourier transform of the heat
kernel associated with $h$.

######  Lemma 2.3.

Let $k_{t}(p,q)$ denote the Fourier transform of the kernel of the operator
$\exp(-th)$ for $t>0$. The function $(p,q)\mapsto k_{t}(p,q)$ is nonnegative
a.e. and satisfies

$\int_{\mathbb{R}^{6}}k_{t}(p,q)\,\mathrm{d}p\,\mathrm{d}q\leq(\pi/t)^{3/2}\,.$
(2.8)

###### Proof.

Let us define $P_{t}(v)$ as $(2\pi)^{-3/2}$ times the Fourier transform of the
function $\exp(-t|x|^{s})$. We have $\int P_{t}(v)\,\mathrm{d}v=1$, and from
Theorem XIII.52 and Example 2 on p. 220 in [64] we know that it is
nonnegative. Accordingly, $P_{t}(v-w)$ defines a probability transition
kernel333The function $P_{t}(v-w)$ is the transition kernel of an $s$-stable
Levy process, see e.g. [46, Example 6.5].. Moreover, we have the Markov
property

$\int_{\mathbb{R}^{3}}P_{t_{1}}(v-z)P_{t_{2}}(z-w)\,\mathrm{d}z=P_{t_{1}+t_{2}}(v-w)$
(2.9)

for $t_{1},t_{2}\geq 0$.

From the Trotter product formula, see e.g. [62, Theorem VII.31], we know that

$\exp(-th)=\lim_{n\to\infty}\left(\exp(-t|x|^{s}/n)\exp(t\Delta/n)\right)^{n}$
(2.10)

holds in the strong operator topology. Let us denote by $K_{n}(v,w)$ the
Fourier transform of the integral kernel of the $n$-dependent operator on the
r.h.s. of the above equation. With the notation $q_{0}=v$ and $q_{n}=w$, we
can write

$\displaystyle 0\leq K_{n}(v,w)=$
$\displaystyle\int\left(\prod_{i=1}^{n}P_{t/n}(q_{i-1}-q_{i})\right)\exp\left(-\sum_{i=1}^{n}\frac{tq_{i}^{2}}{n}\right)\,\mathrm{d}q_{1}\dots\,\mathrm{d}q_{n-1}$
$\displaystyle\leq$
$\displaystyle\int\left(\prod_{i=1}^{n}P_{t/n}(q_{i-1}-q_{i})\right)\left(\frac{1}{n}\sum_{i=1}^{n}\exp\left(-tq_{i}^{2}\right)\right)\,\mathrm{d}q_{1}\dots\,\mathrm{d}q_{n-1}$
$\displaystyle=$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\int
P_{it/n}(v-q)\exp\left(-tq^{2}\right)P_{(n-i)t/n}(q-w)\,\mathrm{d}q\,.$ (2.11)

To come from the first to the second line, we applied Jensen’s inequality. Let
$\phi,\psi\in L^{2}(\mathbb{R}^{3})\cap L^{\infty}(\mathbb{R}^{3})$ be two
nonnegative functions. In combination, (2.10) and (2.11) imply the bound

$\displaystyle
0\leq\int_{\mathbb{R}^{6}}\phi(v)k_{t}(v,w)\psi(w)\,\mathrm{d}v\,\mathrm{d}w$
$\displaystyle=\lim_{n\to\infty}\int_{\mathbb{R}^{6}}\phi(v)K_{n}(v,w)\psi(w)\,\mathrm{d}v\,\mathrm{d}w$
$\displaystyle\leq\liminf_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}\int_{\mathbb{R}^{9}}\phi(v)P_{it/n}(v-q)\exp\left(-tq^{2}\right)P_{(n-i)t/n}(q-w)\psi(w)\,\mathrm{d}v\,\mathrm{d}q\,\mathrm{d}w$
$\displaystyle\leq\|\psi\|_{\infty}\
\|\phi\|_{\infty}\int_{\mathbb{R}^{3}}\exp\left(-tq^{2}\right)\,\mathrm{d}q=\|\psi\|_{\infty}\
\|\phi\|_{\infty}\ \left(\pi/t\right)^{3/2}.$ (2.12)

The first bound in (2.12) implies that $k_{t}(p,q)\geq 0$ holds for a.e.
$(p,q)\in\mathbb{R}^{6}$. The bound in (2.8) follows from (2.12), when we take
the supremum over all functions $\phi,\psi$ with $0\leq\phi,\psi\leq 1$ on
both sides. This proves Lemma 2.3. ∎

We are now ready to show Proposition 2.1.

###### Proof of Proposition 2.1.

To show (2.1), we remark that, with the choice $w(x)=|x|^{s}$, for
$s\in(0;2]$, we have $|\nabla w(x)|=s|x|^{s-1}$ and therefore

$|\nabla w(x)|^{2}\lesssim\big{[}w(x)+|x|^{-2}+1\big{]}\lesssim(1+h)$

with $h=-\Delta+w(x)$. In the last step, we applied Hardy’s inequality. Thus,
we can apply Lemma 2.2 to estimate

$\operatorname{Tr}\,(1-\Delta)^{3/2}\gamma(1-\Delta)^{3/2}\leq\Big{\|}(1-\Delta)^{3/2}\frac{1}{(1+h)^{3/2}}\Big{\|}^{2}\,\operatorname{Tr}\,(1+h)^{3}\gamma\,.$

The claim follows by noticing (with $\\{e_{j}\\}_{j=0}^{\infty}$ indicating
the eigenvalues of $h$) that

$\operatorname{Tr}\,(1+h)^{3}\gamma=\sum_{j=1}^{\infty}\frac{1}{e^{\,\beta(e_{j}-\mu)}-1}(1+e_{j})^{3}\lesssim\operatorname{Tr}e^{-ch}\lesssim\int_{\mathbb{R}^{3}}e^{-cw(x)/2}dx<\infty$
(2.13)

for some constant $c>0$, where we used $e_{1}>e_{0}>\mu$ in the first step,
see [64, Theorem XIII.47], and, in the second step, the bound for the trace of
the propagator of the heat equation below Eq. (6.5) in [13].

The bound $\|\gamma\|\lesssim T_{c}(s)$ follows from $(h-\mu)Q\geq
e_{1}-e_{0}=:\Delta e>0$, which implies that

$\|\gamma\|\leq\frac{1}{\exp(\beta\Delta e)-1}\leq\frac{1}{\beta\Delta
e}\lesssim T_{\mathrm{c}}(s)\,.$

Next, we show (2.2). We use here the notation $\psi_{j}$ for the eigenfunction
of $h$ corresponding to the eigenvalue $e_{j}$, for $j\in\mathbb{N}$. The
identity $(\exp(x)-1)^{-1}=\sum_{\alpha=1}^{\infty}\exp(-\alpha x)$ for $x>0$
allows us to write

$\gamma=\sum_{\alpha=1}^{\infty}\left(\exp(-\beta(h-\mu)\alpha)-\exp(-\beta(e_{0}-\mu)\alpha)|\psi_{0}\rangle\langle\psi_{0}|\right)\,.$

Using the above representation of $\gamma$ and the notation $k_{t}$ from Lemma
2.3, we estimate

$\displaystyle|\hat{\gamma}(p,q)|\leq$
$\displaystyle\sum_{\alpha=1}^{M}\exp(\beta\mu\alpha)k_{\beta\alpha}(p,q)+\sum_{\alpha=M+1}^{\infty}\sum_{i=1}^{\infty}\exp(-\beta\alpha(e_{i}-\mu))|\hat{\psi}_{i}(p)|\
|\hat{\psi}_{i}(q)|$
$\displaystyle+\sum_{\alpha=1}^{M}\exp(-\beta(e_{0}-\mu)\alpha)|\hat{\psi}_{0}(p)|\
|\hat{\psi}_{0}(q)|\,,$ (2.14)

where $M\in\mathbb{N}$. An application of Lemma 2.3 shows that the
$L^{1}(\mathbb{R}^{6})$-norm of the first term on the r.h.s. is bounded by

$\sum_{\alpha=1}^{M}\exp(\beta\mu\alpha)\int_{\mathbb{R}^{6}}k_{\beta\alpha}(p,q)\,\mathrm{d}p\,\mathrm{d}q\leq\sum_{\alpha=1}^{M}\exp(\beta
e_{0}\alpha)\left(\frac{\pi}{\beta\alpha}\right)^{3/2}\leq\left(\pi/\beta\right)^{3/2}\exp(\beta
e_{0}M)\sum_{\alpha=1}^{\infty}\alpha^{-3/2}.$ (2.15)

To obtain the first bound we used $e_{0}-\mu>0$, and for the second one
$e_{0}>0$. The $L^{1}(\mathbb{R}^{6})$-norm of the term in the last line is
bounded by

$\sum_{\alpha=1}^{M}\exp(-\beta(e_{0}-\mu)\alpha)\int_{\mathbb{R}^{6}}|\hat{\psi}_{0}(p)|\
|\hat{\psi}_{0}(q)|\leq M\|\hat{\psi}_{0}\|_{1}^{2}\,.$ (2.16)

The norm on the r.h.s. is bounded by the $H^{2}$-norm of $\psi_{0}$. Since
$\psi_{0}$ is an eigenfunction of $h$, we conclude with Lemma 2.2 that it is
finite. It remains to consider the second term on the r.h.s. of (2.14). To
that end, with $\Delta e=e_{1}-e_{0}>0$ we estimate

$\displaystyle\sum_{\alpha=M+1}^{\infty}\sum_{i=1}^{\infty}$
$\displaystyle\exp(-\beta\alpha(e_{i}-\mu))|\hat{\psi}_{i}(p)|\
|\hat{\psi}_{i}(q)|$
$\displaystyle=\sum_{i=1}^{\infty}\frac{e^{-\beta(M+1)(e_{i}-\mu)}}{1-e^{-\beta(e_{i}-\mu)}}|\hat{\psi}_{i}(p)|\
|\hat{\psi}_{i}(q)|\leq\frac{1}{\beta\Delta e}\sum_{i=1}^{\infty}e^{-\beta
M(e_{i}-\mu)}|\hat{\psi}_{i}(p)|\ |\hat{\psi}_{i}(q)|\,.$ (2.17)

Integrating over $p,q$, the r.h.s. is bounded by

$\displaystyle\frac{1}{\beta\Delta e}\sum_{i=1}^{\infty}e^{-\beta
M(e_{i}-\mu)}\left(\int|\hat{\psi}_{i}(p)|\,\mathrm{d}p\right)^{2}\lesssim\frac{1}{\beta\Delta
e}\sum_{i=1}^{\infty}e^{-\beta
M(e_{i}-\mu)}\,\langle\psi_{i},(1-\Delta)^{2}\psi_{i}\rangle\lesssim\frac{1}{\beta\Delta
e}\sum_{i=1}^{\infty}e^{-\beta M(e_{i}-\mu)}(1+e_{i})^{2}\,.$

Choosing $M\in\mathbb{N}$ with $M>\beta^{-1}$, (2.2) follows from the last
equation (proceeding as in (2.13) to bound the sum over $i$), together with
(2.14), (2.15), and (2.16).

Finally, let us show the bound on $\|\phi\|_{H^{3}}$, for a normalized $\phi$
solving the Hartree equation (2.3). Let $w(x)=|x|^{s}+g(v*|\phi|^{2})(x)$ and
$h_{w}=-\Delta+w$. Since $v\in L^{\infty}(\mathbb{R}^{3})$ (by Sobolev
embedding), we have $\|v*|\phi|^{2}\|_{\infty}\leq\|v\|_{\infty}$. We also
have (choosing $p^{\prime}<3/2$ such that $1/p+1/p^{\prime}=1$)

$\|\nabla(v*|\phi|^{2})\|_{\infty}\leq\|\nabla
v\|_{p}\|\phi\|_{2p^{\prime}}^{2}\lesssim\|\nabla
v\|_{p}\|\phi\|^{2}_{H_{1}}\,.$

Therefore, we can apply Lemma 2.2, with $w=|x|^{s}+g(v*|\phi|^{2})(x)$ and
$\kappa=1+\|v\|_{\infty}$ to show that

$\|(1-\Delta)^{3/2}\phi\|\lesssim\|(\kappa+h_{w})^{3/2}\phi\|=(\kappa+\mu^{\mathrm{H}})^{3/2}\,,$

which proves the claim. ∎

##  3 The Hartree–Fock–Bogoliubov equations

As explained in Section 1.7, we use the HFB equations to define our
fluctuation dynamics in Section 4.2. Here we introduce the equations and
collect some of their properties that we need for the many-body analysis.

For the triple $(\phi_{t},\gamma_{t},\alpha_{t})$ consisting of a condensate
wave function $\phi_{t}\in L^{2}(\mathbb{R}^{3})$, a positive trace class
operator $\gamma_{t}$ (a 1-pdm), and a pairing function $\alpha_{t}\in
L^{2}(\mathbb{R}^{6})$, the HFB equations take the form

$\displaystyle\mathrm{i}\partial_{t}\phi_{t}=$ $\displaystyle
h(\gamma_{t})\phi_{t}+k(\alpha_{t}^{\phi_{t}})\overline{\phi_{t}}$
$\displaystyle\mathrm{i}\partial_{t}\gamma_{t}=$
$\displaystyle\;[h(\gamma_{t}^{\phi_{t}}),\gamma_{t}]+k(\alpha_{t}^{\phi_{t}})\alpha_{t}^{*}-\alpha_{t}k(\alpha_{t}^{\phi_{t}})^{*}$
$\displaystyle\mathrm{i}\partial_{t}\alpha_{t}=$
$\displaystyle\;[h(\gamma_{t}^{\phi_{t}}),\alpha_{t}]_{+}+[k(\alpha_{t}^{\phi_{t}}),\gamma_{t}]_{+}+k(\alpha_{t}^{\phi_{t}}),$
(3.1)

where $[A,B]_{+}=A\overline{B}^{*}+B\overline{A}^{*}$,
$\gamma^{\phi}=\gamma+|\phi\rangle\langle\phi|$, and
$\alpha^{\phi}=\alpha+|\phi\rangle\langle\overline{\phi}|$. By $\overline{A}$
we denote the operator with the complex conjugate integral kernel in position
space. Moreover, we use the notations

$h(\gamma)=-\Delta+b(\gamma),\qquad k(\alpha)=\frac{1}{N}v\sharp\alpha\,,$
(3.2)

where

$b(\gamma)=\frac{1}{N}v\ast\varrho_{\gamma}+\frac{1}{N}v\sharp\gamma\,,$

$v\sharp\sigma$ denotes the operator with kernel $v(x-y)\sigma(x,y)$, and the
density associated with the 1-pdm $\gamma$ is given by
$\varrho_{\gamma}(x)=\gamma(x,x)$. Here $v$ is an interaction potential
satisfying suitable assumptions to be specified later. The mean-field scaling
we are interested in is reflected in the factor $N^{-1}$ multiplying $v$ in
the equations.

The HFB equations arise naturally as a quasi-free approximation of the many-
body dynamics (which does not preserve quasi-freeness), for quasi-free initial
data. For more details concerning this argument and more information about the
HFB equations we refer the reader to [4, 10].

###  3.1 Well-posedness

In this section we state conditions that guarantee the well-posedness of the
HFB equations. We define $M=(1-\Delta)^{1/2}$, and for $j\geq 0$ we introduce
the linear spaces

$\mathcal{H}^{j,1}=M^{-j}\mathcal{L}^{1}M^{-j}\quad\text{ and
}\quad\mathcal{H}^{j,2}=\left\\{\alpha\in\mathcal{L}^{2}\,|\,M^{j}\alpha,\alpha
M^{j}\in\mathcal{L}^{2}\right\\}$

with $\mathcal{L}^{j}$ defined in Section 1.2, which are Banach spaces when
endowed with their natural norms

$\|\gamma\|_{\mathcal{H}^{j,1}}=\|M^{j}\gamma
M^{j}\|_{\mathcal{L}^{1}}\quad\text{ and
}\quad\|\alpha\|_{\mathcal{H}^{j,2}}=\|M^{j}\alpha\|_{\mathcal{L}^{2}}+\|\alpha
M^{j}\|_{\mathcal{L}^{2}}\,,$

respectively. We also define the Banach space $(X^{j},\|\cdot\|_{X^{j}})$,
where

$X^{j}=H^{j}\times\mathcal{H}^{j,1}\times\mathcal{H}^{j,2}\,\quad\text{ and
}\quad\|(\phi,\gamma,\alpha)\|_{X^{j}}=\|\phi\|_{H^{j}}+\|\gamma\|_{\mathcal{H}^{j,1}}+\|\alpha\|_{\mathcal{H}^{j,2}}\,.$

The following holds true, see [4, Theorem 5.1].

######  Theorem 3.1.

Let $(\phi_{0},\gamma_{0},\alpha_{0})\in X^{3}$ and suppose that

$\Gamma_{0}^{(1)}=\begin{pmatrix}\gamma_{0}&\alpha_{0}\\\
\overline{\alpha}_{0}&1+\overline{\gamma}_{0}\end{pmatrix}$

is positive, i.e. $\Gamma_{0}^{(1)}\geq 0$ as an operator on
$L^{2}(\mathbb{R}^{3})\otimes L^{2}(\mathbb{R}^{3})$ (this is always the case
if $\gamma_{0}$, $\alpha_{0}$ are the reduced 1-pdm and pairing density of a
state $\Gamma$ on $\mathscr{F}(L^{2}(\mathbb{R}^{3}))$). Assume also that the
interaction potential $v$ satisfies $v(x)=v(-x)$ and $v\in
W^{1,p}(\mathbb{R}^{3})$ for some $p>3$. Then there exists a unique triple
$(\phi_{t},\gamma_{t},\alpha_{t})\in C^{1}([0,\infty),X^{1})\cap
C^{0}([0,\infty),X^{3})$ satisfying (3.1) with initial datum
$(\phi_{0},\gamma_{0},\alpha_{0})$ in the classical sense in $X^{1}$.
Moreover, the number of particles

$n(\phi_{t},\gamma_{t})=\int_{\mathbb{R}^{3}}|\phi_{t}(x)|^{2}\,\mathrm{d}x+\operatorname{Tr}\gamma_{t}$
(3.3)

and the energy

$\displaystyle\mathcal{E}(\phi_{t},\gamma_{t},\alpha_{t})=$
$\displaystyle\;\operatorname{Tr}[h(\gamma_{t}^{\phi_{t}})\gamma_{t}]+\frac{1}{N}\operatorname{Tr}[(v\ast|\phi_{t}|^{2}+v\sharp|\phi_{t}\rangle\langle\phi_{t}|)\gamma_{t}]$
$\displaystyle\;+\frac{1}{2N}\operatorname{Tr}[(v\ast\varrho_{\gamma_{t}}+v\sharp\gamma_{t})\gamma_{t}]+\frac{1}{N}\int_{\mathbb{R}^{6}}v(x-y)|\alpha_{t}(x,y)+\phi_{t}(x)\phi_{t}(y)|^{2}\,\mathrm{d}x\,\mathrm{d}y$
(3.4)

are conserved along the evolution, i.e.
$n(\phi_{t},\gamma_{t})=n(\phi_{0},\gamma_{0})$,
$\mathcal{E}(\phi_{t},\gamma_{t},\alpha_{t})=\mathcal{E}(\phi_{0},\gamma_{0},\alpha_{0})$,
for every $t>0$. Finally, the dynamics is positivity-preserving, that is,

$\Gamma_{t}^{(1)}=\begin{pmatrix}\gamma_{t}&\alpha_{t}\\\
\overline{\alpha}_{t}&1+\overline{\gamma}_{t}\end{pmatrix}\geq 0$ (3.5)

for every $t>0$.

###  3.2 Approximation of the HFB dynamics

Although we use the HFB equations to define our fluctuation dynamics in
Section 4.2, we ultimately want to show convergence of the full many-body
dynamics to a simpler effective dynamics, with $\phi_{t}$ solving the time-
dependent Hartree equation (1.32), $\gamma_{t}=e^{\mathrm{i}\Delta
t}\gamma_{0}e^{-\mathrm{i}\Delta t}$ evolving freely, and $\alpha_{t}=0$. To
reach this goal, we will use the next proposition.

######  Proposition 3.2.

Let the pair $(\phi_{0},\gamma_{0})$ satisfy Assumption 1.1 with
$s\in(0,\infty)$, and let the interaction potential $v$ satisfy the
assumptions of Theorem 1.2. Moreover, let $(\phi_{t},\gamma_{t},\alpha_{t})$
denote the solution to the HFB equations (3.1) with initial datum
$(\phi_{0},\gamma_{0},0)$. Then there exists a constant $c>0$ independent of
$N$ and $t$ such that

$\displaystyle\left\|\gamma_{t}-\gamma_{t}^{\mathrm{F}}\right\|_{\mathcal{L}^{1}}\lesssim
N^{1/2}T_{\mathrm{c}}^{3/4}(s)t\exp(ct)\,,$
$\displaystyle\left\|\phi_{t}-\phi_{t}^{\mathrm{H}}\right\|\lesssim
T_{\mathrm{c}}^{3/4}(s)t\exp(ct)\,,$ (3.6)

where $\phi_{t}^{\mathrm{H}}$ is the solution to the time-dependent Hartree
equation in (1.32) with initial condition $\phi_{0}$ and
$\gamma_{t}^{\mathrm{F}}=e^{\mathrm{i}\Delta t}\gamma_{0}e^{-\mathrm{i}\Delta
t}$.

######  Remark 3.3.

  1. 1.

Our assumptions on $(\phi_{0},\gamma_{0})$ and $v$ guarantee that the HFB
equations are well-posed, see Theorem 3.1. It is well-known that they also
guarantee the global well-posedness of the time-dependent Hartree equation in
$H^{1}(\mathbb{R}^{3})$, and that the $L^{2}$-mass of the solution is
preserved along the Hartree flow.

  2. 2.

As explained in Section 1.4, the scaling of the critical temperature with the
number of particles is given by $T_{c}(s)\sim
N^{1/\alpha}=N^{\frac{2s}{6+3s}}$. In particular $T_{\mathrm{c}}^{3/4}(s)\ll
N^{1/2}$ for every $s>0$, showing that the bounds given in (3.6) are non-
trivial.

  3. 3.

It is possible to get better rates than those in (3.6) if one compares to a
slightly more complicated evolution for the 1-pdm. Indeed, let
$\omega_{0}=|\phi_{0}\rangle\langle\phi_{0}|+\gamma_{0}$, and define
$\omega_{t}$ to be the solution of the Hartree equation (1.33) with initial
datum $\omega_{t=0}=\omega_{0}$. Then we have the improved bound

$\left\|\
|\phi_{t}\rangle\langle\phi_{t}|+\gamma_{t}-\omega_{t}\right\|_{\mathcal{L}^{1}}\lesssim
N^{1/2}T_{c}^{1/2}(s)t\exp(ct)\,,$ (3.7)

for some constant $c>0$ independent of $t$ and $N$. As a consequence, one gets
a better (indeed, optimal) rate of convergence in (1.31) by replacing
$|\phi_{t}^{\mathrm{H}}\rangle\langle\phi_{t}^{\mathrm{H}}|+e^{i\Delta
t}\gamma_{0}e^{-i\Delta t}$ with $\omega_{t}$ in the statement, see Remark
1.3.3. The estimate (3.7) is proven by a suitable adaptation of the Grönwall-
type argument given below, but since we do not need it for our main result, we
omit the details. We only highlight that the gain from using (1.33) is that it
allows one to replace the direct term $[v\ast|\phi_{t}|^{2},\gamma_{t}]$
appearing on the r.h.s. of (3.11), which we can only bound by
$N^{1/2}T_{\mathrm{c}}^{3/4}(s)$ (see (3.21)), with some simpler terms that
can be shown to be of order $N^{1/2}T_{\mathrm{c}}^{1/2}(s)$.

###### Proof.

We introduce the notation $\theta_{t}=(\phi_{t},\gamma_{t},\alpha_{t})$, and
we write the HFB equations as

$\mathrm{i}\partial_{t}\theta_{t}=A\theta_{t}+f(\theta_{t})\,,$ (3.8)

where

$A\theta=\left(-\Delta\phi,\leavevmode\nobreak\
[-\Delta,\gamma],\leavevmode\nobreak\ [-\Delta,\alpha]_{+}\right)$ (3.9)

is a linear operator on $X^{0}$ with domain $\mathcal{D}(A)=X^{2}$, and $f\in
C(X^{0})$ is defined by

$\begin{split}f(\theta)=\Big{(}b(\gamma)\phi+k(\alpha^{\phi})\overline{\phi},\leavevmode\nobreak\
[b(\gamma^{\phi}),\gamma]+(k(\alpha^{\phi})\alpha^{*}-\alpha
k(\alpha^{\phi})^{*})\gamma,\leavevmode\nobreak\
[b(\gamma^{\phi}),\alpha]_{+}+[k(\alpha^{\phi}),\alpha]_{+}+k(\alpha^{\phi})\Big{)}\,.\end{split}$

Combining (3.8) and (1.32) we get

$\displaystyle\phi_{t}-\phi_{t}^{\mathrm{H}}=$
$\displaystyle-\mathrm{i}\int_{0}^{t}e^{\mathrm{i}(t-s)\Delta}\left(\left(\frac{1}{N}v\ast|\phi_{s}|^{2}\phi_{s}-\frac{1}{N}v\ast|\phi_{s}^{\mathrm{H}}|^{2}\phi_{s}^{\mathrm{H}}\right)+b(\gamma_{s})\phi_{s}+k(\alpha_{s})\overline{\phi_{s}}\right)\,\mathrm{d}s\,,$
(3.10) $\displaystyle\gamma_{t}-\gamma_{t}^{\mathrm{F}}=$
$\displaystyle-\mathrm{i}\int_{0}^{t}e^{\mathrm{i}(t-s)\Delta}\Big{(}[b(\gamma_{s}^{\phi_{s}}),\gamma_{s}]+k(\alpha_{s}^{\phi_{s}})\alpha_{s}^{*}-\alpha_{s}k(\alpha_{s}^{\phi_{s}})^{*}\Big{)}e^{-\mathrm{i}(t-s)\Delta}\,\mathrm{d}s\,,$
(3.11) $\displaystyle\alpha_{t}=$
$\displaystyle-\mathrm{i}\int_{0}^{t}e^{\mathrm{i}(t-s)(\Delta_{x}+\Delta_{y})}\Big{(}[b(\gamma_{s}^{\phi_{s}}),\alpha_{s}]_{+}+[k(\alpha_{s}^{\phi_{s}}),\gamma_{s}]_{+}+k(\alpha_{s}^{\phi_{s}})\Big{)}\,\mathrm{d}s\,.$
(3.12)

Here, we identify $\alpha_{t}$ with its kernel $\alpha_{t}(x,y)$ in
$L^{2}(\mathbb{R}^{6})$ and we denote by $\Delta_{x}+\Delta_{y}$ the Laplacian
acting on this space.

Let us start by deriving a bound for the r.h.s. of (3.10). We use the identity

$\frac{1}{N}v\ast|\phi_{s}|^{2}\phi_{s}-\frac{1}{N}v\ast|\phi_{s}^{\mathrm{H}}|^{2}\phi_{s}^{\mathrm{H}}=\left(\frac{1}{N}v\ast(|\phi_{s}|^{2}-|\phi_{s}^{\mathrm{H}}|^{2})\right)\phi_{s}^{\mathrm{H}}+\left(\frac{1}{N}v\ast|\phi_{s}|^{2}\right)(\phi_{s}-\phi_{s}^{\mathrm{H}})\,,$
(3.13)

Young’s inequality, and the conservation of the particle number for the
Hartree and HFB equations, to see that

$\begin{split}\left\|\frac{1}{N}v\ast|\phi_{s}|^{2}\phi_{s}-\frac{1}{N}v\ast|\phi_{s}^{\mathrm{H}}|^{2}\phi_{s}^{\mathrm{H}}\right\|\lesssim&N^{-1}\|v\|_{\infty}\left(\|\phi_{s}\|^{2}+\|\phi_{s}^{\mathrm{H}}\|^{2}\right)\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|\\\
\lesssim&N^{-1}\|v\|_{\infty}\left(\|\phi_{0}\|^{2}+\|\gamma_{0}\|_{\mathcal{L}^{1}}\right)\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|\lesssim\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|\,.\end{split}$

In the second term on the r.h.s. of (3.10) we write

$b(\gamma_{s})\phi_{s}=b(\gamma_{s}^{\mathrm{F}})\phi_{s}+b(\gamma_{s}-\gamma_{s}^{\mathrm{F}})\phi_{s}\,,$

and we estimate the two terms separately. First, we bound the operator norm of
$b(\gamma_{s}^{\mathrm{F}})=N^{-1}v\ast\varrho_{\gamma_{s}^{\mathrm{F}}}+N^{-1}v\sharp\gamma_{s}^{\mathrm{F}}$.
Using Assumption 1.1 (B), we see that

$\|N^{-1}v\ast\varrho_{\gamma_{s}^{\mathrm{F}}}\|=N^{-1}\|v\ast\varrho_{\gamma_{s}^{\mathrm{F}}}\|_{\infty}\leq
N^{-1}\|v\|_{1}\sup_{x\in\mathbb{R}^{3}}|\gamma_{s}^{\mathrm{F}}(x,x)|\lesssim
N^{-1}T_{\mathrm{c}}^{3/2}(s)\,.$ (3.14)

Moreover, for any $\psi\in L^{2}(\mathbb{R}^{3})$ we can estimate

$\begin{split}N^{-1}|(v\sharp\gamma_{s}^{\mathrm{F}})\psi(x)|\leq&N^{-1}\int_{\mathbb{R}^{3}}|v(x-y)|\
|\gamma_{s}^{\mathrm{F}}(x,y)|\ |\psi(y)|\,\mathrm{d}y\\\
\leq&N^{-1}\sup_{x,y\in\mathbb{R}^{3}}|\gamma_{s}^{\mathrm{F}}(x,y)|\int_{\mathbb{R}^{3}}|v(x-y)|\
|\psi(y)|\,\mathrm{d}y\lesssim
N^{-1}T_{\mathrm{c}}^{3/2}(s)(|v|\ast|\psi|)(x)\,,\end{split}$

which implies

$N^{-1}\|(v\sharp\gamma_{s}^{\mathrm{F}})\psi\|\lesssim
N^{-1}T_{\mathrm{c}}^{3/2}(s)\|\ |v|\ast|\psi|\ \|\lesssim
N^{-1}T_{\mathrm{c}}^{3/2}(s)\|\psi\|\,.$ (3.15)

For the operator norm of $b(\gamma_{s}-\gamma_{s}^{\mathrm{F}})$ we have the
bound

$\displaystyle
N^{-1}\|v\ast(\varrho_{\gamma_{s}}-\varrho_{\gamma_{s}^{\mathrm{F}}})\|=$
$\displaystyle
N^{-1}\|v\ast(\varrho_{\gamma_{s}}-\varrho_{\gamma_{s}^{\mathrm{F}}})\|_{\infty}$
$\displaystyle\leq$ $\displaystyle
N^{-1}\|v\|_{\infty}\int_{\mathbb{R}^{3}}|\gamma_{s}(x,x)-\gamma_{s}^{\mathrm{F}}(x,x)|\,\mathrm{d}x\lesssim
N^{-1}\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\,,$ (3.16)

and

$\displaystyle N^{-1}\|v\sharp(\gamma_{s}-\gamma_{s}^{\mathrm{F}})\|\leq$
$\displaystyle
N^{-1}\left(\int_{\mathbb{R}^{6}}v(x-y)^{2}|\gamma_{s}(x,y)-\gamma_{s}^{\mathrm{F}}(x,y)|^{2}\,\mathrm{d}x\,\mathrm{d}y\right)^{1/2}$
$\displaystyle\leq$ $\displaystyle
N^{-1}\|v\|_{\infty}\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{2}}\lesssim
N^{-1}\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\,,$ (3.17)

where in the last step we used the fact that the Hilbert-Schmidt norm of a
trace-class operator is bounded by its trace norm. Collecting (3.14)–(3.17) we
deduce

$\displaystyle\|b(\gamma_{s}^{\mathrm{F}})\|\lesssim$ $\displaystyle
N^{-1}T_{\mathrm{c}}^{3/2}(s)$
$\displaystyle\|b(\gamma_{s}-\gamma_{s}^{\mathrm{F}})\|\lesssim$
$\displaystyle
N^{-1}\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\,.$ (3.18)

The above bounds and the conservation of the particle number imply

$\begin{split}\|b(\gamma_{s})\phi_{s}\|\lesssim&N^{-1/2}\left(\|\phi_{0}\|^{2}+\|\gamma_{0}\|_{\mathcal{L}^{1}}\right)^{1/2}\left(N^{-1/2}T_{\mathrm{c}}^{3/2}(s)\|v\|_{1}+N^{-1/2}\|v\|_{\infty}\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\right)\\\
\lesssim&N^{-1/2}T_{\mathrm{c}}^{3/2}(s)\|v\|_{1}+N^{-1/2}\|v\|_{\infty}\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\,.\end{split}$

Similarly, we estimate

$\|k(\alpha_{s})\phi_{s}\|\leq
N^{-1}\|v\|_{\infty}\|\alpha_{s}\|_{\mathcal{L}^{2}}\|\phi_{s}\|\lesssim
N^{-1/2}\|v\|_{\infty}\|\alpha_{s}\|_{\mathcal{L}^{2}}\,.$

We have thus shown that

$\|\phi_{t}-\phi_{t}^{\mathrm{H}}\|\lesssim
N^{-1/2}T_{\mathrm{c}}^{3/2}(s)t+\int_{0}^{t}\left(\|\phi_{s}-\phi_{s}^{\mathrm{F}}\|+N^{-1/2}\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}+N^{-1/2}\|\alpha_{s}\|_{\mathcal{L}^{2}}\right)\,\mathrm{d}s\,.$
(3.19)

We now turn to the r.h.s. of (3.11). Using again the bounds in (3.18), we find

$\displaystyle\|[b(\gamma_{s}^{\phi_{s}}),\gamma_{s}]\|_{\mathcal{L}^{1}}\leq$
$\displaystyle\|b(\gamma_{s}^{\phi_{s}})\gamma_{s}\|_{\mathcal{L}^{1}}+\|\gamma_{s}b(\gamma_{s}^{\phi_{s}})\|_{\mathcal{L}^{1}}$
$\displaystyle\leq$ $\displaystyle
2\Big{(}\|\gamma_{s}\|_{\mathcal{L}^{1}}\left(\|b(\gamma_{s}-\gamma_{s}^{\mathrm{F}})\|+\|b(\gamma_{s}^{\mathrm{F}})\|\right)+\|b(|\phi_{s}\rangle\langle\phi_{s}|)\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}+\|b(|\phi_{s}\rangle\langle\phi_{s}|)\|\
\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\Big{)}$
$\displaystyle\lesssim$
$\displaystyle\|v\|_{\infty}\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}+T_{\mathrm{c}}^{3/2}(s)\|v\|_{1}$
$\displaystyle+N^{-1}\|(v\ast|\phi_{s}|^{2})\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}+N^{-1}\|(v\sharp|\phi_{s}\rangle\langle\phi_{s}|)\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}+\|b(|\phi_{s}\rangle\langle\phi_{s}|)\|\
\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\,.$ (3.20)

The first term in the last line is bounded by

$\displaystyle
N^{-1}\|(v\ast|\phi_{s}|^{2})\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\leq$
$\displaystyle
N^{-1}\|\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}^{1/2}\|(v\ast|\phi_{s}|^{2})(\gamma_{s}^{\mathrm{F}})^{1/2}\|_{\mathcal{L}^{2}}$
$\displaystyle\lesssim$ $\displaystyle
N^{-1/2}\left(\int_{\mathbb{R}^{3}}|(v\ast|\phi_{s}|^{2})(x)|^{2}|\gamma_{s}^{\mathrm{F}}(x,x)|\,\mathrm{d}x\right)^{1/2}$
$\displaystyle\lesssim$ $\displaystyle
N^{-1/2}T_{\mathrm{c}}^{3/4}(s)\|v\ast|\phi_{s}|^{2}\|\lesssim
N^{1/2}T_{\mathrm{c}}^{3/4}(s)\|v\|\lesssim N^{1/2}T_{\mathrm{c}}^{3/4}(s)\,.$
(3.21)

The exchange term can be estimated by

$\|v\sharp|\phi_{s}\rangle\langle\phi_{s}|\|_{\mathcal{L}_{1}}\leq\int\hat{v}(p)\|\
|e^{ip\cdot.}\phi_{s}\rangle\langle e^{ip\cdot.}\phi_{s}|\
\|_{\mathcal{L}_{1}}\,\mathrm{d}p\leq\|\hat{v}\|_{1}\|\phi_{s}\|^{2},$

thus,

$N^{-1}\|v\sharp|\phi_{s}\rangle\langle\phi_{s}|\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}_{1}}\lesssim\|\gamma_{s}^{\mathrm{F}}\|\lesssim
T_{c}(s)\,.$

As for the last term on the last line of (3.20), we have

$\begin{split}\|b(|\phi_{s}\rangle\langle\phi_{s}|)\|\
\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\leq
N^{-1}\left(\|v\ast|\phi_{s}|^{2}\|_{\infty}+\|v\sharp|\phi_{s}\rangle\langle\phi_{s}|\
\|\right)\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\lesssim\|v\|_{\infty}\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\,.\end{split}$

Going back to (3.20) we find

$\|[b(\gamma_{s}^{\phi_{s}}),\gamma_{s}]\|_{\mathcal{L}^{1}}\lesssim\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}+N^{1/2}T_{\mathrm{c}}^{3/4}(s)+T_{\mathrm{c}}^{3/2}(s)+T_{\mathrm{c}}(s)\,.$

To bound the last two terms in the integral on the r.h.s. of (3.11), we use
the elementary inequality
$\|\alpha_{s}\|_{\mathcal{L}^{2}}^{2}\leq(1+\operatorname{Tr}\gamma_{s})\operatorname{Tr}\gamma_{s}$,
which follows from (3.5). We find

$\|k(\alpha_{s}^{\phi_{s}})\alpha_{s}^{*}\|_{\mathcal{L}^{1}}\leq\|k(\alpha_{s}^{\phi_{s}})\|_{\mathcal{L}^{2}}\|\alpha_{s}\|_{\mathcal{L}^{2}}\leq
N^{-1}\|v\|_{\infty}\left(\|\phi_{s}\|^{2}+(1+\|\gamma_{s}\|_{\mathcal{L}^{1}}^{2})^{1/2}\right)\|\alpha_{s}\|_{\mathcal{L}^{2}}\lesssim\|\alpha_{s}\|_{\mathcal{L}^{2}}\,,$
(3.22)

and similarly

$\|\alpha_{s}k(\alpha_{s}^{\phi_{s}})^{*}\|_{\mathcal{L}^{1}}\lesssim\|\alpha_{s}\|_{\mathcal{L}^{2}}\,.$

We conclude that

$\|\gamma_{t}-\gamma_{t}^{\mathrm{F}}\|_{\mathcal{L}^{1}}\lesssim\left(N^{1/2}T_{\mathrm{c}}^{3/4}(s)+T_{\mathrm{c}}^{3/2}(s)+T_{\mathrm{c}}(s)\right)t+\int_{0}^{t}\left(\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|_{\mathcal{L}^{1}}+\|\alpha\|_{\mathcal{L}^{2}}\right)\,\mathrm{d}s\,.$
(3.23)

We also need to control the operator norm of (3.11). Thus, we observe that

$\displaystyle\|b(\gamma_{s}^{\phi_{s}})\|\leq$ $\displaystyle
N^{-1}\|v\ast(\varrho_{\gamma_{s}}+|\phi_{s}|^{2})\|_{\infty}+N^{-1}\|v\sharp(\gamma_{s}+|\phi_{s}\rangle\langle\phi_{s}|)\|_{\mathcal{L}^{2}}$
$\displaystyle\lesssim$ $\displaystyle
N^{-1}\|v\|_{\infty}\left(\|\gamma_{s}\|_{\mathcal{L}^{1}}+\|\phi_{s}\|^{2}\right)\lesssim\|v\|_{\infty}\,,$
(3.24)

which implies

$\|[b(\gamma_{s}^{\phi_{s}}),\gamma_{s}]\|\leq 2\|b(\gamma_{s}^{\phi_{s}})\|\
\|\gamma_{s}\|\lesssim\left(\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|+T_{\mathrm{c}}(s)\right)\,.$

Similarly, we find

$\|k(\alpha_{s}^{\phi_{s}})\alpha_{s}^{*}\|\lesssim\|k(\alpha_{s}^{\phi_{s}})\alpha_{s}^{*}\|_{\mathcal{L}^{1}}\lesssim\|\alpha_{s}\|_{\mathcal{L}^{2}}\,,$

where we used (3.22) in the second step. We conclude that

$\|\gamma_{t}-\gamma_{t}^{\mathrm{F}}\|\lesssim
T_{\mathrm{c}}(s)t+\int_{0}^{t}\left(\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|+\|\alpha_{s}\|_{\mathcal{L}^{2}}\right)\,\mathrm{d}s\,.$
(3.25)

Finally, we combine (3.12), (3.24),
$\|k(\alpha_{s}^{\phi_{s}})\|\leq\|k(\alpha_{s}^{\phi_{s}})\|_{\mathcal{L}^{2}}\lesssim\|v\|_{\infty}$,
and $\|AB\|_{\mathcal{L}^{2}}\leq\|A\|\ \|B\|_{\mathcal{L}^{2}}$ to see that

$\displaystyle\|\alpha_{t}\|_{\mathcal{L}^{2}}\leq$
$\displaystyle\int_{0}^{t}\left\|[b(\gamma_{s}^{\phi_{s}}),\alpha_{s}]_{+}+[k(\alpha_{s}^{\phi_{s}}),\gamma_{s}]_{+}+k(\alpha_{s}^{\phi_{s}})\right\|_{\mathcal{L}^{2}}\,\mathrm{d}s$
$\displaystyle\lesssim$ $\displaystyle
T_{\mathrm{c}}(s)t+\int_{0}^{t}\left(\|\gamma_{s}-\gamma_{s}^{\mathrm{F}}\|+\|\alpha_{s}\|_{\mathcal{L}^{2}}\right)\,\mathrm{d}s$
(3.26)

holds.

Let us define the norm

$\|\theta\|_{*}=T_{\mathrm{c}}^{-3/4}(s)\|\phi\|+N^{-1/2}T_{\mathrm{c}}^{-3/4}(s)\|\gamma\|_{\mathcal{L}^{1}}+T_{\mathrm{c}}^{-1}(s)(\|\gamma\|+\|\alpha\|_{\mathcal{L}^{2}})\,.$

In combination, (3.19),(3.23), (3.25), and (3.26) imply

$\|\theta_{t}-\theta_{t,0}\|_{*}\leq
Ct+K\int_{0}^{t}\|\theta_{s}-\theta_{s,0}\|_{*}\,,$

for constants $C,K>0$ independent of $t$ and $N$. Here,
$\theta_{t}=(\phi_{t},\gamma_{t},\alpha_{t})$, while
$\theta_{t,0}=(\phi_{t}^{\mathrm{H}},\gamma_{t}^{\mathrm{F}},0)$. With an
application of Grönwall’s Lemma we conclude that

$\|\theta_{t}-\theta_{t,0}\|_{*}\leq Cte^{Kt}\,,$

which proves the claim. ∎

Assumption 1.1 guarantees the diluteness of the thermal cloud when propagated
with the free time evolution at all times. We now prove that this property
continues to hold for the 1-pdm evolved with the HFB equations up to a double
exponential growth in time. The precise statement is captured in the following
proposition.

######  Proposition 3.4.

Let the pair $(\phi_{0},\gamma_{0})$ satisfy Assumption 1.1 with $s\in(0,2]$,
and let $v$ be as in Theorem 1.2. Moreover, let
$(\phi_{t},\gamma_{t},\alpha_{t})$ denote the solution to the HFB equations in
(3.1), with initial datum $(\phi_{0},\gamma_{0},0)$. Then there exist
constants $c,C>0$ s.t.

$\sup_{x,y\in\mathbb{R}^{3}}|\gamma_{t}(x,y)|\lesssim
T_{\mathrm{c}}^{{3/2}}(s)\exp(c\exp(ct))$ (3.27)

holds for $t>0$ and $N\geq C\exp(c\exp(ct))$.

In order to prove Proposition 3.4 we need two lemmas. The first one concerns
the growth of the $L^{1}$-norm of the Fourier transform of the solution to the
time-dependent Hartree equation in time.

######  Lemma 3.5.

Assume that $\phi_{t}$ is a solution to the time-dependent Hartree equation in
(1.32) with initial condition $\phi_{0}\in H^{2}(\mathbb{R}^{3})$ and $v$
satisfies the same assumptions as in Theorem 1.2. There exists a constant
$c>0$ s.t.

$\|\widehat{\phi}_{t}\|_{1}\lesssim\|\widehat{\phi}_{0}\|_{1}\exp(ct)\,,$

where $\widehat{\phi}_{t}$ denotes the Fourier transform of $\phi_{t}$.

###### Proof.

We use Duhamel’s formula to write $\phi_{t}$ as

$\phi_{t}=e^{\mathrm{i}\Delta
t}\phi_{0}-\mathrm{i}\int_{0}^{t}e^{\mathrm{i}(t-s)\Delta}\left[\left(N^{-1}v\ast|\phi_{s}|^{2}\right)\phi_{s}\right]\,\mathrm{d}s\,.$

When take the $L^{1}$-norm in Fourier space on both sides, we find

$\begin{split}\|\widehat{\phi}_{t}\|_{1}\leq&\|\widehat{\phi}_{0}\|_{1}+N^{-1}\int_{0}^{t}\|\hat{v}\
\widehat{|\phi_{s}|^{2}}\|_{1}\|\widehat{\phi}_{s}\|_{1}\,\mathrm{d}s\\\
\leq&\|\widehat{\phi}_{0}\|_{1}+N^{-1}\|\hat{v}\|_{1}\int_{0}^{t}\|\
\widehat{|\phi_{s}|^{2}}\|_{\infty}\|\widehat{\phi}_{s}\|_{1}\,\mathrm{d}s\\\
\leq&\|\widehat{\phi}_{0}\|_{1}+N^{-1}\|\hat{v}\|_{1}\|\phi_{0}\|^{2}\int_{0}^{t}\|\widehat{\phi}_{s}\|_{1}\,\mathrm{d}s\lesssim\|\widehat{\phi}_{0}\|_{1}+\int_{0}^{t}\|\widehat{\phi}_{s}\|_{1}\,\mathrm{d}s\,,\end{split}$

where we used the conservation of the $L^{2}$-mass in the last step. An
application of Grönwall’s Lemma yields

$\|\widehat{\phi}_{s}\|_{1}\lesssim\|\widehat{\phi}_{0}\|_{1}\exp(ct)\,,$

which proves the claim. ∎

######  Remark 3.6.

If $v\geq 0$ or $\hat{v}\geq 0$ it can be shown that there exists a constant
$C>0$ depending only on $\|\phi_{0}\|_{2}$, $\mathcal{E}^{\mathrm{H}}(\phi)$,
with $\mathcal{E}^{\mathrm{H}}$ defined in (1.12), and $\varepsilon>0$ s.t.

$\|\phi_{t}\|_{H^{2}(\mathbb{R}^{3})}\leq
C(1+|t|^{2+\varepsilon})\|\phi_{0}\|_{H^{2}(\mathbb{R}^{3})}\,.$ (3.28)

This result follows from the analysis in [67, Section 5]. It allows us to
replace the exponential time dependence in Lemma 3.5 by the polynomial time-
dependence on the r.h.s. of (3.28) (because the $H^{2}(\mathbb{R}^{3})$-norm
of a function dominates the $L^{1}(\mathbb{R}^{3})$-norm of its Fourier
transform), the double exponential time-dependence in Proposition 3.4 by
$\exp(c(t+t^{3+\varepsilon}))$, and the triple exponential time-dependence in
Theorem 1.2 by $\exp(c\exp(c(t+t^{3+\varepsilon})))$.

The second lemma is a generalization of Grönwall’s inequality which follows
from [70, Corollary 3.2].

######  Lemma 3.7.

Let $p\in\mathbb{N}\setminus\\{0,1\\}$, $T>0$ and let $a,b,c$ be continuous,
nonnegative and nondecreasing functions on $[0,T]$ such that
$e^{a(t)(p-1)T}<1+b(t)/(a(t)^{p-1}c(t))$, for every $t\in[0,T]$. If $u\in
L^{\infty}[0,T]$ satisfies $u(t)\geq 0$ and

$u(t)\leq
a(t)+b(t)\int_{0}^{t}u(s)\,\mathrm{d}s+c(t)\int_{0}^{t}u^{p}(s)\,\mathrm{d}s$

for a.e. $t\in[0,T]$, then, for a.e. $t\in[0,T]$, we have

$u(t)\leq
a(t)\left[1-\frac{a(t)^{p-1}c(t)}{b(t)}\left(e^{(p-1)b(t)t}-1\right)\right]^{-\frac{1}{p-1}}e^{b(t)t}\,.$

######  Remark 3.8.

As a consequence of the power-law nonlinearity in the integral inequality in
Lemma 3.7, the bound for $u(t)$ blows up as $t$ approaches the smallest time
$T^{*}$ such that

$T^{*}=b(T^{*})^{-1}(p-1)^{-1}\ln(1+b(T^{*})/(a(T^{*})^{p-1}c(T^{*})))\,.$

For times satisfying $0<t\leq T^{*}/2$, the nonlinearity is irrelevant and
functions obeying the integral inequality satisfy the standard Grönwall
estimate

$u(t)\lesssim a(t)\exp(b(t)t)\,.$

We are now prepared to give the proof of Proposition 3.4.

###### Proof of Proposition 3.4.

For a Hilbert-Schmidt operator $\gamma$ we use the notation
$\|\gamma\|_{d}=\int|\hat{\gamma}(p,q)|\,\mathrm{d}p\,\mathrm{d}q$, where
$\hat{\gamma}(p,q)$ denotes the Fourier transform of the integral kernel of
$\gamma$. For $\phi\in L^{2}(\mathbb{R}^{3})$ we also define
$\|\phi\|_{\mathrm{d}}=\|\widehat{\phi}\|_{1}$. We will prove that there
exists a constant $c>0$, independent of $t$ and $N$, such that

$\displaystyle\|\gamma_{t}\|_{\mathrm{d}}\lesssim$ $\displaystyle
T_{\mathrm{c}}^{3/2}(s)\exp(c\exp(ct))\,,$
$\displaystyle\|\alpha_{t}\|_{\mathrm{d}}\lesssim$ $\displaystyle
T_{\mathrm{c}}^{3/2}(s)\exp(c\exp(ct))\,,$
$\displaystyle\|\phi_{t}-\phi_{t}^{\mathrm{H}}\|_{\mathrm{d}}\lesssim$
$\displaystyle\exp(c\exp(ct))\,,$

holds for $t>0$, and $N$ large enough depending on $t$, which implies the
claim. To this end, we start by noticing that the operator
$e^{-\mathrm{i}A(t-s)}$, with $A$ as in (3.9), acts in Fourier space by
multiplication with a complex phase, and therefore preserves the norm
$\|\cdot\|_{\mathrm{d}}$ for $\phi$, $\gamma$, and $\alpha$, for every
$t>s>0$. Thus (3.11) implies, for $t>0$,

$\|\gamma_{t}\|_{\mathrm{d}}\leq\|\gamma_{t}^{\mathrm{F}}\|_{\mathrm{d}}+\int_{0}^{t}\left\|b(\gamma_{s}^{\phi_{s}})\gamma_{s}+k(\alpha_{s}^{\phi_{s}})\alpha_{s}^{*}-\alpha_{s}k(\alpha_{s}^{\phi_{s}})^{*}\right\|_{\mathrm{d}}\,\mathrm{d}s\,,$
(3.29)

with $\gamma^{\mathrm{F}}_{t}$ defined below (3.6). Young’s convolution
inequality implies

$\left\|(N^{-1}v\ast\varrho_{\gamma_{s}^{\phi_{s}}})\gamma_{s}\right\|_{\mathrm{d}}=N^{-1}\left\|\left(\hat{v}\widehat{\varrho_{\gamma_{s}^{\phi_{s}}}}\right)\ast_{1}\hat{\gamma}_{s}(\cdot,\cdot)\right\|_{1}\leq
N^{-1}\|\hat{v}\|_{1}\left(\|\phi_{0}\|^{2}+\|\gamma_{0}\|_{1}\right)\|\gamma_{s}\|_{\mathrm{d}}\lesssim\|\gamma_{s}\|_{\mathrm{d}}\,,$

where $f\ast_{i}\hat{\gamma}_{s}(\cdot,\cdot)$ denotes the convolution of $f$
with the $i$-th variable of the kernel $\hat{\gamma}_{s}(p,q)$. A
straightforward computation shows that, for any Hilbert-Schmidt operator
$\gamma$, we have

$\widehat{(v\sharp\gamma)}(p,q)=(2\pi)^{-3}\int_{\mathbb{R}^{3}}\hat{v}(\xi)\hat{\gamma}(q-\xi,p-\xi)\,\mathrm{d}\xi\,.$
(3.30)

It follows that

$\begin{split}\|N^{-1}\left(v\sharp\gamma_{s}^{\phi_{s}}\right)\gamma_{s}\|_{\mathrm{d}}=&N^{-1}\int_{\mathbb{R}^{6}}\left|\int_{\mathbb{R}^{3}}\hat{v}(\xi)\left(\hat{\gamma}_{s}(z-\xi,p-\xi)+\overline{\widehat{\phi}_{s}(z-\xi)}\widehat{\phi}_{s}(p-\xi)\right)\hat{\gamma}_{s}(z,q)\,\mathrm{d}\xi\,\mathrm{d}z\right|\,\mathrm{d}p\,\mathrm{d}q\\\
\leq&N^{-1}\|\gamma_{s}\|_{\mathrm{d}}\cdot\left(\sup_{z\in\mathbb{R}^{3}}\int_{\mathbb{R}^{6}}|\hat{v}(\xi)|\left(|\hat{\gamma}_{s}(z-\xi,p-\xi)|+|\widehat{\phi}_{s}(z-\xi)\widehat{\phi}_{s}(p-\xi)|\right)\,\mathrm{d}\xi\,\mathrm{d}p\right)\,.\end{split}$

Bounding $|\hat{v}(\xi)|\leq\|v\|_{1}$ and integrating first in $p$ and then
in $\xi$, we find

$\begin{split}\|N^{-1}\left(v\sharp\gamma_{s}^{\phi_{s}}\right)\gamma_{s}\|_{\mathrm{d}}\leq&N^{-1}\|v\|_{1}\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\phi_{s}\|_{\mathrm{d}}^{2}\right)\|\gamma_{s}\|_{\mathrm{d}}\\\
\lesssim&N^{-1}\|v\|_{1}\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}+\|\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}\right)\|\gamma_{s}\|_{\mathrm{d}}\\\
\lesssim&N^{-1}\|\gamma_{s}\|_{\mathrm{d}}^{2}+N^{-1}\|\gamma_{s}\|_{\mathrm{d}}\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}+\exp(ct)\|\gamma_{s}\|_{\mathrm{d}}\,.\end{split}$

In the last step we used Lemma 3.5. Similarly, we see that

$\begin{split}\left\|k(\alpha^{\phi_{s}})\alpha_{s}^{*}-\alpha_{s}k(\alpha^{\phi_{s}})^{*}\right\|_{\mathrm{d}}\leq&N^{-1}\|v\|_{1}\left(\|\alpha_{s}\|_{\mathrm{d}}+\|\phi_{s}\|_{\mathrm{d}}^{2}\right)\|\alpha_{s}\|_{\mathrm{d}}\\\
\leq&N^{-1}\|v\|_{1}\left(\|\alpha_{s}\|_{\mathrm{d}}+\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}+\|\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}\right)\|\alpha_{s}\|_{\mathrm{d}}\\\
\lesssim&N^{-1}\|\alpha_{s}\|_{\mathrm{d}}^{2}+N^{-1}\|\alpha_{s}\|_{\mathrm{d}}\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}+\exp(ct)\|\alpha_{s}\|_{\mathrm{d}}\,.\end{split}$

In combination, (3.29) and the above bounds imply

$\displaystyle\|\gamma_{t}\|_{\mathrm{d}}\leq$
$\displaystyle\|\gamma_{t}^{\mathrm{F}}\|_{\mathrm{d}}+C\exp(ct)\int_{0}^{t}\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)\,\mathrm{d}s$
$\displaystyle+CN^{-1}\int_{0}^{t}\left[\|\gamma_{s}\|_{\mathrm{d}}^{2}+\|\alpha_{s}\|_{\mathrm{d}}^{2}+\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)\right]\,\mathrm{d}s\,.$
(3.31)

Next, we consider the r.h.s. of (3.12). Estimates similars to the ones above
show

$\displaystyle\left\|[b(\gamma_{s}^{\phi_{s}}),\alpha_{s}]_{+}+[k(\alpha_{s}^{\phi_{s}}),\gamma_{s}]_{+}\right\|_{\mathrm{d}}\lesssim$
$\displaystyle
N^{-1}\|v\|_{1}\Big{(}\left(\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}+\|\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}\right)\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)+\|\gamma_{s}\|_{\mathrm{d}}\|\alpha_{s}\|_{\mathrm{d}}\Big{)}$
$\displaystyle\lesssim$ $\displaystyle
N^{-1}\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}+\exp(ct)\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)+N^{-1}\|\gamma_{s}\|_{\mathrm{d}}\|\alpha_{s}\|_{\mathrm{d}}$

and

$\|k(\alpha_{s}^{\phi_{s}})\|_{\mathrm{d}}\leq
N^{-1}\|\hat{v}\|_{1}\left(\|\alpha_{s}\|_{\mathrm{d}}+\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}+\|\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}\right)\lesssim
N^{-1}\|\alpha_{s}\|_{\mathrm{d}}+N^{-1}\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}+\exp(ct)\,,$

which implies

$\displaystyle\|\alpha_{t}\|_{\mathrm{d}}\lesssim$
$\displaystyle\exp(ct)+\exp(ct)\int_{0}^{t}\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)\,\mathrm{d}s$
$\displaystyle+N^{-1}\int_{0}^{t}\left[\|\gamma_{s}\|_{\mathrm{d}}^{2}+\|\alpha_{s}\|_{\mathrm{d}}^{2}+\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}^{2}\left(1+\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)\right]\,\mathrm{d}s\,.$
(3.32)

Finally, we consider (3.10). We use the identity

$|\phi_{s}|^{2}-|\phi_{s}^{\mathrm{H}}|^{2}=\phi_{s}(\overline{\phi_{s}-\phi_{s}^{\mathrm{H}}})+(\phi_{s}-\phi_{s}^{\mathrm{H}})\overline{\phi_{s}^{\mathrm{H}}}\,,$

Lemma 3.5, and Young’s inequality to see that

$\begin{split}\left\|\left(N^{-1}v\ast(|\phi_{s}|^{2}-|\phi_{s}^{\mathrm{H}}|^{2})\right)\phi_{s}^{\mathrm{H}}\right\|_{\mathrm{d}}\leq&N^{-1}\left\|\hat{v}\left(\left(\widehat{\overline{\phi_{s}}}-\widehat{\overline{\phi_{s}^{\mathrm{H}}}}\right)\ast\widehat{\phi_{s}}+\left(\phi_{s}-\phi_{s}^{\mathrm{H}}\right)\widehat{\overline{\phi_{s}^{\mathrm{H}}}}\right)\right\|_{1}\|\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}\\\
\lesssim&N^{-1/2}\exp(ct)\|v\|_{2}\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}\left(\|\phi_{s}\|+\|\phi_{s}^{\mathrm{H}}\|\right)\lesssim\exp(ct)\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}\,,\end{split}$

as well as

$\begin{split}\left\|\left(N^{-1}v\ast|\phi_{s}|^{2}\right)(\phi_{s}-\phi_{s}^{\mathrm{H}})\right\|_{\mathrm{d}}\leq&N^{-1}\|\hat{v}\widehat{|\phi_{s}|^{2}}\|_{1}\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}\lesssim\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}\,.\end{split}$

Out of the remaining terms on the r.h.s. of (3.10), the direct term is bounded
by

$\begin{split}\|N^{-1}v\ast\varrho_{\gamma_{s}}\phi_{s}\|_{\mathrm{d}}=&N^{-1}\|(\hat{v}\widehat{\varrho_{\gamma_{s}}})\ast\widehat{\phi}_{s}\|_{1}\leq
N^{-1}\|\hat{v}\widehat{\varrho_{\gamma_{s}}}\|_{1}\|\widehat{\phi}_{s}\|_{1}\,.\end{split}$

We use the identity
$\widehat{\varrho_{\gamma_{s}}}(\xi)=(2\pi)^{-3/2}\int\hat{\gamma}_{s}(p,p-\xi)\,\mathrm{d}p$
to see that

$\|\widehat{\varrho_{\gamma_{s}}}\|_{1}\leq(2\pi)^{-3/2}\int_{\mathbb{R}^{6}}|\hat{\gamma}_{s}(p,p-\xi)|\,\mathrm{d}p\,\mathrm{d}\xi\lesssim\|\gamma_{s}\|_{\mathrm{d}}\,,$

which implies

$\|N^{-1}v\ast\varrho_{\gamma_{s}}\phi_{s}\|_{\mathrm{d}}\lesssim
N^{-1}\|v\|_{1}\|\gamma_{s}\|_{\mathrm{d}}\left(\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}+\|\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}\right)\lesssim
N^{-1}\|\gamma_{s}\|_{\mathrm{d}}\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}+N^{-1/2}\exp(ct)\|\gamma_{s}\|_{\mathrm{d}}\,.$

Using (3.30) another time, we can bound the remaining term on the r.h.s. of
(3.10) by

$\begin{split}\|N^{-1}v\sharp(\gamma_{s}+\alpha_{s})\phi_{s}\|_{\mathrm{d}}\leq&N^{-1}\int_{\mathbb{R}^{3}}\left|\int_{\mathbb{R}^{6}}\hat{v}(\xi)\Big{(}\hat{\gamma}(q-\xi,p-\xi)+\hat{\alpha}(q-\xi,p-\xi)\Big{)}\widehat{\phi}_{s}(q)\,\mathrm{d}\xi\,\mathrm{d}q\right|\,\mathrm{d}p\\\
\leq&N^{-1}\|v\|_{1}\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)\left(\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}+\|\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}\right)\\\
\lesssim&N^{-1}\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}+N^{-1/2}\exp(ct)\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)\,.\end{split}$

In combination, (3.10) and the above estimates imply

$\|\phi_{t}-\phi_{t}^{\mathrm{H}}\|_{\mathrm{d}}\lesssim\exp(ct)\int_{0}^{t}\left[\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}+N^{-1/2}\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)\right]\,\mathrm{d}s+N^{-1}\int_{0}^{t}\|\phi_{s}-\phi_{s}^{\mathrm{H}}\|_{\mathrm{d}}\left(\|\gamma_{s}\|_{\mathrm{d}}+\|\alpha_{s}\|_{\mathrm{d}}\right)\,\mathrm{d}s\,.$
(3.33)

Let us define
$u(t)=T_{\mathrm{c}}^{-3/2}(s)\|\gamma_{t}\|_{\mathrm{d}}+T_{\mathrm{c}}^{-3/2}(s)\|\alpha_{t}\|_{\mathrm{d}}+\|\phi_{t}-\phi_{t}^{\mathrm{H}}\|_{\mathrm{d}}$.
We collect the bounds in (3.31), (3.32), and (3.33), and use the assumption
$s\leq 2$ to find

$\begin{split}u(t)\lesssim&1+\exp(ct)T_{\mathrm{c}}^{-3/2}(s)+\exp(ct)\int_{0}^{t}u(z)\,\mathrm{d}z+N^{-1}T_{\mathrm{c}}^{3/2}(s)\int_{0}^{t}u(z)^{2}\,\mathrm{d}z+N^{-1}\int_{0}^{t}u(z)^{3}\,\mathrm{d}z\\\
\lesssim&\exp(ct)+\exp(ct)\int_{0}^{t}u(z)\,\mathrm{d}z+N^{-1}\int_{0}^{t}u(z)^{3}\,\mathrm{d}z\,,\end{split}$

for every $t>0$ and some $c>0$ independent of $t$ and $N$. For any $t>0$, the
assumptions of Lemma 3.7 are satisfied with $T=t$, as long as $N\geq
C\exp(c\exp(ct))$, and we conclude that, under this condition,

$u(t)\lesssim\exp(c\exp(ct))\,,$

which proves the claim. ∎

##  4 Fluctuation dynamics

As explained in Section 1.7, the proof of Theorem 1.2 is based on the
construction of a suitable fluctuation dynamics on the double Fock space
$\mathscr{F}({\mathfrak{h}}\oplus{\mathfrak{h}})$ with
${\mathfrak{h}}=L^{2}(\mathbb{R}^{3})$. In this section we construct the
fluctuation dynamics and discuss some of its properties. We start by recalling
a few well known facts about Weyl operators and Bogoliubov transformations.

###  4.1 Weyl operators and Bogoliubov transformations

We recall the definition of the Weyl operator $W(\phi)$ acting on
$\mathscr{F}({\mathfrak{h}})$ in (1.17) as well as the Weyl operator
$\mathcal{W}(\phi)$ acting on
$\mathscr{F}({\mathfrak{h}}\oplus{\mathfrak{h}})$ in (1.26). Both operators
act as shifts on creation and annihilation operators, i.e.

$W(\phi)^{*}a_{x}W(\phi)=a_{x}+\phi(x)\,,\qquad
W(\phi)^{*}a_{x}^{*}W(\phi)=a_{x}^{*}+\overline{\phi(x)}\,,$ (4.1)

and

$\displaystyle\mathcal{W}(\phi)^{*}a_{\ell,x}\mathcal{W}(\phi)$
$\displaystyle=a_{\ell,x}+\phi(x)\,,\qquad\mathcal{W}(\phi)^{*}a_{\ell,x}^{*}\mathcal{W}(\phi)=a_{\ell,x}^{*}+\overline{\phi(x)}\,,$
$\displaystyle\mathcal{W}(\phi)^{*}a_{r,x}\mathcal{W}(\phi)$
$\displaystyle=a_{r,x}+\overline{\phi(x)}\,,\qquad\mathcal{W}(\phi)^{*}a_{r,x}^{*}\mathcal{W}(\phi)=a_{r,x}^{*}+\phi(x)\,.$
(4.2)

To express (4.2) in a more compact form, we denote
$\mathbf{x}=(\sigma,x)\in\\{\ell,r\\}\times\mathbb{R}^{3}$ and for a given
$\phi\in{\mathfrak{h}}$ we define
$\boldsymbol{\phi}:\\{\ell,r\\}\times\mathbb{R}^{3}\to\mathbb{C}$ by

$\boldsymbol{\phi}(\mathbf{x})=\boldsymbol{\phi}(\sigma,x)=\begin{cases}\phi(x),&\sigma=\ell\,,\\\
\overline{\phi(x)},&\sigma=r\,.\end{cases}$

This allows us to write (4.2) as

$\mathcal{W}(\phi)^{*}a_{\mathbf{x}}\mathcal{W}(\phi)=a_{\mathbf{x}}+\boldsymbol{\phi}(\mathbf{x})\,,\qquad\mathcal{W}(\phi)^{*}a^{*}_{\mathbf{x}}\mathcal{W}(\phi)=a^{*}_{\mathbf{x}}+\overline{\boldsymbol{\phi}(\mathbf{x})}\,.$
(4.3)

Next we recall some well known facts about Bogoliubov transformations. We are
interested in the case where the one-particle Hilbert space is given by
${\mathfrak{h}}$ or by
${\mathfrak{h}}\oplus{\mathfrak{h}}={\mathfrak{h}}^{2}$. The complex conjugate
$\overline{A}$ of an operator $A\in\mathcal{B}({\mathfrak{h}}^{n})$ with
$n\in\\{1,2\\}$ is defined as the operator whose integral kernel, in position
space, is the complex conjugate of the kernel of $A$. For $n=1,2$ a bounded
linear map
$\nu:({\mathfrak{h}}^{n}\oplus{\mathfrak{h}}^{n})\to({\mathfrak{h}}^{n}\oplus{\mathfrak{h}}^{n})$
of the form

$\nu=\begin{pmatrix}U&\overline{V}\\\ V&\overline{U}\end{pmatrix}$ (4.4)

is called a symplectomorphism if

$\nu^{*}\mathcal{S}\nu=\mathcal{S}\quad\text{ and
}\quad\nu\mathcal{S}\nu^{*}=\mathcal{S}\quad\text{ with
}\quad\mathcal{S}=\begin{pmatrix}\mathds{1}&0\\\
0&-\mathds{1}\end{pmatrix}\,.$ (4.5)

Here $\mathds{1}$ denotes the identity operator on ${\mathfrak{h}}^{n}$. If
$V$ is a Hilbert–Schmidt operator, $\nu$ is implementable on
$\mathscr{F}({\mathfrak{h}}^{n})$, see [4, 68], i.e. there exists a unitary
$T_{\nu}$ on $\mathscr{F}({\mathfrak{h}}^{n})$ such that

$\displaystyle T_{\nu}^{*}a(f)T_{\nu}$
$\displaystyle=a(Uf)+a^{*}(\overline{V}\ \overline{f})\,,$ $\displaystyle
T_{\nu}^{*}a^{*}(f)T_{\nu}$ $\displaystyle=a^{*}(Uf)+a(\overline{V}\
\overline{f})\,,$ (4.6)

hold for $f\in{\mathfrak{h}}^{n}$. The operator $T_{\nu}$ is called the
Bogoliubov transformation corresponding to $\nu$.

For Bogoliubov transformations on $\mathscr{F}({\mathfrak{h}}^{2})$ it is
convenient to introduce a compact notation, similar to (4.3). For a
symplectomorphism
$\nu:({\mathfrak{h}}^{2}\oplus{\mathfrak{h}}^{2})\to({\mathfrak{h}}^{2}\oplus{\mathfrak{h}}^{2})$
of the form (4.4), with $U,V\in\mathcal{B}({\mathfrak{h}}^{2})$, we denote by
$U(\mathbf{x},\mathbf{y})$, $V(\mathbf{x},\mathbf{y})$ the kernels444Since
$V\in\mathcal{L}^{2}({\mathfrak{h}}^{2})$ its kernel is a function in
$L^{2}(\mathbb{R}^{3})\oplus L^{2}(\mathbb{R}^{3})$. The kernel of $U$ is a
distribution. of $U,V$, that is,

$U(\mathbf{x},\mathbf{y})=U_{\sigma,\sigma^{\prime}}(x,y)\,,$

where $\mathbf{x}=(\sigma,x)$ and $\mathbf{y}=(\sigma^{\prime},y)$, and the
same for $V$. We also define
$U_{\mathbf{x}}(\mathbf{y})=U(\mathbf{y},\mathbf{x})$,
$V_{\mathbf{x}}(\mathbf{y})=V(\mathbf{y},\mathbf{x})$. With this notation, we
can write

$\displaystyle T_{\nu}^{*}a_{\mathbf{x}}T_{\nu}$
$\displaystyle=a(U_{\mathbf{x}})+a^{*}(\overline{V_{\mathbf{x}}})\,,$
$\displaystyle T_{\nu}^{*}a_{\mathbf{x}}^{*}T_{\nu}$
$\displaystyle=a^{*}(U_{\mathbf{x}})+a(\overline{V_{\mathbf{x}}})\,.$ (4.7)

###  4.2 Construction of the fluctuation dynamics

In this section we use the Hartree–Fock–Bogoliubov (HFB) equations in (3.1) to
define our fluctuation dynamics. As explained in Section 1.7, we do not intend
to derive the HFB equations, but rather use them as an intermediate step
between the many-body evolution and the simple effective dynamics appearing in
Theorem 1.2 (or in Remark 1.3.3).

First of all, we discuss how our reference state, associated with the solution
of the HFB equations, can be conveniently expressed as a vector in
$\mathscr{F}({\mathfrak{h}}^{2})$. Let $(\phi,\gamma)$ satisfy Assumption 1.1
and let $(\phi_{t},\gamma_{t},\alpha_{t})$ be the solution of the HFB
equations (3.1), with initial datum $(\phi,\gamma,0)$. We denote by $G_{t}$
the density matrix of the quasi-free state on $\mathscr{F}(\mathfrak{h})$,
with 1-pdm $\gamma_{t}$ and pairing density $\alpha_{t}$, and we define its
generalized 1-pdm by

$\Gamma_{t}^{(1)}=\begin{pmatrix}\gamma_{t}&\alpha_{t}\\\
\overline{\alpha}_{t}&1+\overline{\gamma}_{t}\end{pmatrix}$ (4.8)

(the condensate will be added later on, through conjugation with a Weyl
operator). As proven in [4, Proposition 3.9], we can write

$\Gamma_{t}^{(1)}=\mathcal{U}_{t}^{*}\Gamma_{0}^{(1)}\mathcal{U}_{t}$ (4.9)

with

$\Gamma_{0}^{(1)}=\begin{pmatrix}\gamma&0\\\
0&1+\overline{\gamma}\end{pmatrix}$

and an implementable symplectomorphism

$\mathcal{U}_{t}=\begin{pmatrix}\pi_{t}&\overline{\vartheta}_{t}\\\
\vartheta_{t}&\overline{\pi}_{t}\end{pmatrix}$ (4.10)

satisfying the equation

$\mathrm{i}\partial_{t}\mathcal{U}^{*}_{t}=\mathcal{S}\Lambda(\Gamma_{t}^{(1)})\mathcal{U}_{t}^{*}\quad\text{
with
}\quad\Lambda(\Gamma_{t}^{(1)})=\begin{pmatrix}h(\gamma_{t}^{\phi_{t}})&k(\alpha_{t}^{\phi_{t}})\\\
\overline{k}(\alpha_{t}^{\phi_{t}})&\overline{h}(\gamma_{t}^{\phi_{t}})\end{pmatrix}\,,$
(4.11)

with $h(\gamma),k(\alpha)$ defined in (3.2) and the initial condition
$\mathcal{U}_{0}=\mathds{1}$ (with $\mathcal{S}$ defined in (4.5)). Raising
$\mathcal{U}_{t}$ to ${\mathfrak{h}}^{2}$, we obtain the map
$R_{t}:{\mathfrak{h}}^{2}\to{\mathfrak{h}}^{2}$, defined by

$R_{t}=\begin{pmatrix}\Pi_{t}&\overline{\Theta}_{t}\\\
\Theta_{t}&\overline{\Pi}_{t}\end{pmatrix}\,,$

where we defined

$\Pi_{t}=\begin{pmatrix}\pi_{t}&0\\\
0&\overline{\pi}_{t}\end{pmatrix}\quad\text{ and
}\quad\Theta_{t}=\begin{pmatrix}\vartheta_{t}&0\\\
0&\overline{\vartheta}_{t}\end{pmatrix}\,.$

As $\mathcal{U}_{t}$, the map $R_{t}$ is an implementable symplectomorphism,
for all $t\in\mathbb{R}$. Hence, we find a family of unitary operators
$\mathcal{R}_{t}$ on $\mathscr{F}({\mathfrak{h}}^{2})$, with

$\mathcal{R}_{t}^{*}a_{\mathbf{x}}\mathcal{R}_{t}=a(\Pi_{t,{\mathbf{x}}})+a^{*}(\overline{\Theta_{t,\mathbf{x}}})\quad\text{
and
}\quad\mathcal{R}_{t}^{*}a^{*}_{\mathbf{x}}\mathcal{R}_{t}=a^{*}(\Pi_{t,{\mathbf{x}}})+a(\overline{\Theta_{t,\mathbf{x}}})\,.$

Let us also introduce the unitary family
$\mathcal{T}_{t}=\mathcal{R}_{t}\mathcal{T}(\gamma)$, with
$\mathcal{T}(\gamma)$ as defined in (1.23). We observe that $\mathcal{T}_{t}$
is again a family of Bogoliubov transformations, satisfying

$\mathcal{T}^{*}_{t}a_{\mathbf{x}}\mathcal{T}_{t}=a(U_{t,{\mathbf{x}}})+a^{*}(\overline{V_{t,\mathbf{x}}})\quad\text{
and
}\quad\mathcal{T}^{*}_{t}a^{*}_{\mathbf{x}}\mathcal{T}_{t}=a^{*}(U_{t,{\mathbf{x}}})+a(\overline{V_{t,\mathbf{x}}})$
(4.12)

with

$U_{t}=U_{0}\Pi_{t}+\overline{V}_{0}\Theta_{t}\,,\quad\quad
V_{t}=V_{0}\Pi_{t}+\overline{U}_{0}\Theta_{t}\,,$ (4.13)

and

$U_{0}=\begin{pmatrix}\sqrt{1+\gamma}&0\\\
0&\overline{\sqrt{1+\gamma}}\end{pmatrix},\quad\quad
V_{0}=\begin{pmatrix}0&\overline{\sqrt{\gamma}}\\\
\sqrt{\gamma}&0\end{pmatrix}\,.$ (4.14)

In other words, $\mathcal{T}_{t}$ implements the symplectomorphism

$T_{t}=T_{0}R_{t}=\begin{pmatrix}U_{t}&\overline{V}_{t}\\\
V_{t}&\overline{U}_{t}\end{pmatrix}\,.$ (4.15)

The following lemma provides us with bounds for the operator norms of $U_{t}$
and $V_{t}$.

######  Lemma 4.1.

Under the assumptions of Proposition 3.2, there exists a constant $c>0$ such
that

$\|U_{t}\|\lesssim\sqrt{T_{\mathrm{c}}(s)}\exp(ct)\quad\text{ and
}\quad\|V_{t}\|\lesssim\sqrt{T_{\mathrm{c}}(s)}\exp(ct)$ (4.16)

hold. Here $\|\cdot\|$ denotes the norm of a linear operator on
${\mathfrak{h}}\oplus{\mathfrak{h}}$.

###### Proof.

We first prove a bound on the operators $\pi_{t}$, $\vartheta_{t}$ in (4.10),
from which the claim will follow. To do this we observe that, as a consequence
of (4.11), $\pi_{t}$, $\vartheta_{t}$ satisfy the system of ODEs

$\displaystyle\mathrm{i}\partial_{t}\pi_{t}^{*}$
$\displaystyle=h(\gamma_{t}^{\phi_{t}})\pi_{t}^{*}+k(\alpha_{t}^{\phi_{t}})\overline{\vartheta_{t}^{*}}\,,$
$\displaystyle\mathrm{i}\partial_{t}\vartheta_{t}^{*}$
$\displaystyle=h(\gamma_{t}^{\phi_{t}})\vartheta_{t}^{*}+k(\alpha_{t}^{\phi_{t}})\overline{\pi_{t}^{*}}\,.$
(4.17)

Let $\psi$ with $\|\psi\|\leq 1$ belong to an appropriate dense subset of
$L^{2}(\mathbb{R}^{3})$. Using (4.17) and the fact that
$h(\gamma_{t}^{\phi_{t}})$ is self-adjoint we get

$\mathrm{i}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|\pi_{t}\psi\|^{2}=\mathrm{i}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\langle\psi,\pi_{t}\pi_{t}^{*}\psi\rangle=\langle\psi,\left(\pi_{t}k(\alpha_{t}^{\phi_{t}})\overline{\vartheta_{t}^{*}}-\overline{\vartheta_{t}}k(\alpha_{t}^{\phi_{t}})^{*}\pi_{t}^{*}\right)\psi\rangle=2\mathrm{i}\
\textrm{Im}\langle\psi,\pi_{t}k(\alpha_{t}^{\phi_{t}})\overline{\vartheta_{t}^{*}}\psi\rangle\,,$

thus we can estimate

$\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|\pi_{t}\psi\|^{2}\leq
2\|k(\alpha_{t}^{\phi_{t}})\|\left(\|\pi_{t}\psi\|^{2}+\|\vartheta_{t}\psi\|^{2}\right)\,.$
(4.18)

Proceeding in the same way, we also find

$\frac{\,\mathrm{d}}{\,\mathrm{d}t}\|\vartheta_{t}\psi\|^{2}\leq
2\|k(\alpha_{t}^{\phi_{t}})\|\left(\|\pi_{t}\psi\|^{2}+\|\vartheta_{t}\psi\|^{2}\right)\,.$
(4.19)

The operator norm of $k(\alpha_{t}^{\phi_{t}})$ is easily bounded by

$\|k(\alpha_{t}^{\phi_{t}})\|\leq\|k(\alpha_{t}^{\phi_{t}})\|_{\mathcal{L}^{2}}\leq
N^{-1}\|v\|_{\infty}\left(\int\left(|\alpha_{t}(x,y)|^{2}+|\phi_{t}(x)|^{2}|\phi_{t}(y)|^{2}\right)\,\mathrm{d}x\,\mathrm{d}y\right)^{1/2}\leq
c/2\,,$

for some constant $c>0$ independent of $t$, $N$ and $\psi$. To obtain the
bound, we used
$\|\alpha_{t}\|_{\mathcal{L}^{2}}^{2}\leq\|\gamma_{t}\|_{\mathcal{L}^{1}}^{2}+\|\gamma_{t}\|_{\mathcal{L}^{1}}\leq
N(N+1)$ and $\|\phi_{t}\|^{2}\leq N$. Inserting this into (4.18) and (4.19) we
find

$\frac{\,\mathrm{d}}{\,\mathrm{d}t}\left(\|\pi_{t}\psi\|^{2}+\|\vartheta_{t}\psi\|^{2}\right)\leq
c\left(\|\pi_{t}\psi\|^{2}+\|\vartheta_{t}\psi\|^{2}\right)\,,$

which together with the initial condition $\pi_{0}=\mathds{1}$,
$\vartheta_{0}=0$ and an application of Grönwall’s Lemma lets us conclude

$\|\pi_{t}\psi\|^{2}\leq\exp(ct)\,,\qquad\|\vartheta_{t}\psi\|^{2}\leq\exp(ct)\,.$
(4.20)

Therefore, taking the $\sup$ over $\psi$ we find

$\|\pi_{t}\|^{2}\leq\exp(ct)\,,\qquad\|\vartheta_{t}\|^{2}\leq\exp(ct)\,.$
(4.21)

Using (4.21), (4.13), (4.14) and point (C) of Assumption 1.1, we deduce
(4.16). ∎

We claim that $\mathcal{T}_{t}\Omega\in\mathscr{F}({\mathfrak{h}}^{2})$
describes exactly the quasi-free state on $\mathscr{F}({\mathfrak{h}})$ with
generalized 1-pdm (4.8). In fact,

$\displaystyle\langle\mathcal{U}^{*}\mathcal{T}_{t}\Omega,\left(a^{*}(g)a(f)\otimes\mathds{1}\right)\mathcal{U}^{*}\mathcal{T}_{t}\Omega\rangle$
$\displaystyle=\langle\mathcal{T}(\gamma)\Omega,\left(a_{\ell}^{*}(\pi_{t}g)+a_{\ell}(\overline{\vartheta_{t}}\overline{g})\right)\left(a_{\ell}(\pi_{t}f)+a^{*}_{\ell}(\overline{\vartheta_{t}}\overline{f})\right)\mathcal{T}(\gamma)\Omega\rangle$
$\displaystyle=\langle
f,\left[\pi_{t}^{*}\gamma\pi_{t}+\vartheta_{t}^{*}(\mathds{1}+\overline{\gamma})\vartheta_{t}\right]g\rangle$
$\displaystyle=\langle f,\gamma_{t}g\rangle\,,$ (4.22)

with $\mathcal{U}$ in (1.22). To come to the last line we used (4.9). A
similar computation shows

$\langle\mathcal{U}^{*}\mathcal{T}_{t}\Omega,\left(a(g)a(f)\otimes\mathds{1}\right)\mathcal{U}^{*}\mathcal{T}_{t}\Omega\rangle=\langle
f,\alpha_{t}\overline{g}\rangle\,.$ (4.23)

Next, we add the condensate. To this end, we apply the Weyl operator
$\mathcal{W}(\phi_{t})$. Proceeding as in (4.22) and (4.23), and using (4.2),
we find that the vector
$\mathcal{W}(\phi_{t})\mathcal{T}_{t}\Omega\in\mathscr{F}({\mathfrak{h}}^{2})$
describes the quasi-free state associated with the solution of the HFB
equations. With this reference state, we are now ready to define the
fluctuation dynamics as

$\mathcal{U}^{\mathrm{fluct}}(t,s)=\mathcal{T}_{t}^{*}\mathcal{W}(\phi_{t})^{*}\exp(-\mathrm{i}\mathcal{L}_{N}(t-s))\mathcal{W}(\phi_{s})\mathcal{T}_{s}$
(4.24)

with $\mathcal{L}_{N}$ in (1.30).

###  4.3 Bound by the expectation of the number of particles

For $\xi\in\mathscr{F}({\mathfrak{h}}^{2})$, we define the fluctuation vector
$\xi_{t}=\mathcal{U}^{\mathrm{fluct}}(t,0)\xi$. By definition, we have

$\mathcal{W}(\phi_{t})\mathcal{T}_{t}\xi_{t}=e^{-\mathrm{i}\mathcal{L}_{N}t}\mathcal{W}(\phi)\mathcal{T}_{0}\xi=e^{-\mathrm{i}\mathcal{L}_{N}t}\mathcal{W}(\phi)\mathcal{T}(\gamma)\xi\,,$
(4.25)

which is exactly the vector in $\mathscr{F}({\mathfrak{h}}^{2})$ associated
with the state $\Gamma_{\xi,t}(\phi,\gamma)$ defined in (1.29) and considered
in Theorem 1.2. Hence, the 1-pdm of the state $\Gamma_{\xi,t}(\phi,\gamma)$
can be written as

$\gamma_{\xi,t}(x,y)=\operatorname{Tr}[a^{*}_{y}a_{x}\Gamma_{\xi,t}(\phi,\gamma)]=\langle\xi_{t},\mathcal{T}^{*}_{t}\mathcal{W}^{*}(\phi_{t})a^{*}_{\ell,y}a_{\ell,x}\mathcal{W}(\phi_{t})\mathcal{T}_{t}\xi_{t}\rangle\,.$

A short computation that uses (4.3), (4.7) and

$\gamma_{t}(x,y)=\sum_{\sigma\in\\{\ell,r\\}}\int_{\mathbb{R}^{3}}\overline{V_{t,\ell}((\sigma,z),x)}V_{t,\ell}((\sigma,z),y)\,\mathrm{d}z$
(4.26)

allows us to write the kernel of $\gamma_{\xi,t}$ as

$\displaystyle\gamma_{\xi,t}(x,y)$
$\displaystyle=\phi_{t}(x)\overline{\phi_{t}(y)}+\gamma_{t}(x,y)+\langle\xi_{t},\left(a^{*}(U_{t,\ell,y})a(U_{t,\ell,x})+a^{*}(\overline{V_{t,\ell,x}})a(\overline{V_{t,\ell,y}})\right)\xi_{t}\rangle$
$\displaystyle+\langle\xi_{t},\left(a(\overline{V_{t,\ell,y}})a(U_{t,\ell,x})+a^{*}(U_{t,\ell,y})a^{*}(\overline{V_{t,\ell,x}})\right)\xi_{t}\rangle$
$\displaystyle+\langle\xi_{t},\left\\{\phi_{t}(x)\left(a^{*}(U_{t,\ell,y})+a(\overline{V_{t,\ell,y}})\right)+\overline{\phi_{t}(y)}\left(a(U_{t,\ell,x})+a^{*}(\overline{V_{t,\ell,x}})\right)\right\\}\xi_{t}\rangle\,.$
(4.27)

With this representation for $\gamma_{\xi,t}$, we obtain the following lemma,
which bounds the trace norm of
$\gamma_{\xi,t}-|\phi_{t}\rangle\langle\phi_{t}|-\gamma_{t}$ in terms of the
expectation of the number of particles operator

$\mathcal{N}=\mathcal{N}_{\ell}+\mathcal{N}_{r}+5=\,\mathrm{d}\Gamma(\mathds{1})+5\,.$
(4.28)

######  Lemma 4.2.

Under the assumptions of Theorem 1.2 and with $U_{t},V_{t}$ as in (4.13), we
have the bound

$\displaystyle\|\gamma_{\xi,t}-|\phi_{t}\rangle\langle\phi_{t}|-\gamma_{t}\|_{\mathcal{L}^{1}}\leq$
$\displaystyle\;(\|U_{t}\|^{2}+\|V_{t}\|^{2})\langle\xi_{t},\mathcal{N}\xi_{t}\rangle+2N^{1/2}(\|U_{t}\|^{2}+\|V_{t}\|^{2})^{1/2}\langle\xi_{t},\mathcal{N}\xi_{t}\rangle^{1/2}$
$\displaystyle+2\|\phi_{t}\|\cdot(\|U_{t}\|^{2}+\|V_{t}\|^{2})^{1/2}\langle\xi_{t},\mathcal{N}\xi_{t}\rangle^{1/2}\,.$
(4.29)

###### Proof.

Let $J\in\mathcal{B}(h)$. We take the trace of (4.27) against $J$: this will
allow us to prove (4.29) by duality. Denoting by $J(x,y)$ the integral kernel
of $J$, the contribution of the third term on the r.h.s. of (4.27) is given by

$\displaystyle\int_{\mathbb{R}^{6}}J(x,y)\left\langle\xi_{t},\left[a^{*}(U_{t,\ell,y})a(U_{t,\ell,x})+a^{*}(\overline{V_{t,\ell,x}})a(\overline{V_{t,\ell,y}})\right]\xi_{t}\right\rangle\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle\hskip
28.45274pt=\int\left[\overline{U_{t,\ell}(\mathbf{z}_{2},x)}J(x,y)U_{t,\ell}(\mathbf{z}_{1},y)+\overline{V_{t,\ell}(\mathbf{z}_{2},x)}J(x,y)V_{t,\ell}(\mathbf{z}_{1},y)\right]\left\langle\xi_{t},a^{*}_{\mathbf{z}_{1}}a_{\mathbf{z}_{2}}\xi_{t}\right\rangle\,\mathrm{d}\mathbf{z}_{1}\,\mathrm{d}\mathbf{z}_{2}\text{d}x\,\mathrm{d}y$
$\displaystyle\hskip
28.45274pt=\left\langle\xi_{t},\,\mathrm{d}\Gamma(J_{1}+J_{2})\xi_{t}\right\rangle\,,$
(4.30)

where we used the notation
$\int\,\mathrm{d}\mathbf{z}=\sum_{\sigma\in\\{\ell,r\\}}\int_{\mathbb{R}^{3}}\,\mathrm{d}z$.
The operators $J_{1}$ and $J_{2}$ are defined by

$\displaystyle J_{1}(\mathbf{z}_{1},\mathbf{z}_{2})=$
$\displaystyle\int_{\mathbb{R}^{6}}\overline{U_{t,\ell}(\mathbf{z}_{2},x)}J(x,y)\overline{U_{t,\ell}^{*}(y,\mathbf{z}_{1})}\text{d}x\,\mathrm{d}y,\quad
J_{2}(\mathbf{z}_{1},\mathbf{z}_{2})=\int_{\mathbb{R}^{6}}\overline{V_{t,\ell}(\mathbf{z}_{2},x)}J(x,y)\overline{V_{t,\ell}^{*}(y,\mathbf{z}_{1})}\text{d}x\,\mathrm{d}y\,.$
(4.31)

The absolute value of the r.h.s. of (4.30) is bounded by

$\displaystyle\left|\left\langle\xi_{t},\text{d}\Gamma(J_{1}+J_{2})\xi_{t}\right\rangle\right|\leq(\|J_{1}\|+\|J_{2}\|)\left\langle\xi_{t},\mathcal{N}\xi_{t}\right\rangle\,.$
(4.32)

With

$\|J_{1}\|\leq\|U_{t}\|^{2}\|J\|\quad\text{ and
}\quad\|J_{1}\|\leq\|V_{t}\|^{2}\|J\|$ (4.33)

we conclude that

$\left|\int_{\mathbb{R}^{6}}J(x,y)\left\langle\xi_{t},\left[a^{*}(U_{t,\ell,y})a(U_{t,\ell,x})+a^{*}(\overline{V_{t,\ell,x}})a(\overline{V_{t,\ell,y}})\right]\xi_{t}\right\rangle\text{d}x\,\mathrm{d}y\right|\leq\|J\|\left(\|U_{t}\|^{2}+\|V_{t}\|^{2}\right)\left\langle\xi_{t},\mathcal{N}\xi_{t}\right\rangle\,.$
(4.34)

Next we consider the contribution of the fourth term on the r.h.s. of (4.27).

Defining $J_{3},J_{4}:{\mathfrak{h}}\to{\mathfrak{h}}\oplus{\mathfrak{h}}$ by

$J_{3}(\mathbf{z},y)=\int_{\mathbb{R}^{3}}U_{t,l}(\mathbf{z},x)\overline{J(x,y)}\,\mathrm{d}x\,,\quad\quad
J_{4}(\mathbf{z},x)=\int_{\mathbb{R}^{3}}J(x,y)U_{t,l}(\mathbf{z},y)\,\mathrm{d}y\,,$
(4.35)

we write

$\displaystyle\int_{\mathbb{R}^{6}}J(x,y)\langle\xi_{t},\left[a(\overline{V_{t,\ell,y}})a(U_{t,\ell,x})+a^{*}(U_{t,\ell,y})a^{*}(\overline{V_{t,\ell,x}})\right]\xi_{t}\rangle\,\mathrm{d}x\,\mathrm{d}y$
$\displaystyle\hskip
85.35826pt=\int_{\mathbb{R}^{3}}\langle\xi_{t},a(\overline{V_{t,\ell,y}})a(J_{3,y})\xi_{t}\rangle\,\mathrm{d}y+\int_{\mathbb{R}^{3}}\langle\xi_{t},a^{*}(J_{4,x})a^{*}(\overline{V_{t,\ell,x}})\xi_{t}\rangle\,\mathrm{d}x\,.$
(4.36)

The absolute value of the first term on the r.h.s. of (4.36) is bounded by

$\displaystyle\left(\int_{\mathbb{R}^{3}}\langle\xi_{t},a(\overline{V_{t,\ell,y}})a^{*}(\overline{V_{t,\ell,y}})\xi_{t}\rangle\,\mathrm{d}y\right)^{1/2}\left(\int_{\mathbb{R}^{3}}\langle\xi_{t},a^{*}(J_{3,y})a(J_{3,y})\xi_{t}\rangle\,\mathrm{d}y\right)^{1/2}\,.$
(4.37)

Using the canonical commutation relations, (4.26), and
$\operatorname{Tr}\gamma_{t}\leq N$, we see that the first factor is not
larger than

$\|V_{t}\|\ \langle\xi_{t},\mathcal{N}_{r}\xi_{t}\rangle^{1/2}+N^{1/2}\,.$
(4.38)

For the second factor we have the bound

$\|J_{3}\|\ \langle\xi_{t},\mathcal{N}\xi_{t}\rangle^{1/2}\leq\|J\|\
\|U_{t}\|\ \langle\xi_{t},\mathcal{N}\xi_{t}\rangle^{1/2}\,.$ (4.39)

The second term on the r.h.s. of (4.36) can be bounded similarly by the same
expression, and we thus find

$\displaystyle\left|\int_{\mathbb{R}^{6}}J(x,y)\langle\xi_{t},\left[a(\overline{V_{t,\ell,y}})a(U_{t,\ell,x})+a^{*}(U_{t,\ell,y})a^{*}(\overline{V_{t,\ell,x}})\right]\xi_{t}\rangle\,\mathrm{d}x\,\mathrm{d}y\right|$
$\displaystyle\hskip 142.26378pt\leq 2\left(\|J\|\ \|U_{t}\|\ \|V_{t}\|\
\langle\xi_{t},\mathcal{N}\xi_{t}\rangle+N^{1/2}\|J\|\ \|U_{t}\|\
\langle\xi_{t},\mathcal{N}\xi_{t}\rangle^{1/2}\right)\,.$ (4.40)

It remains to provide bounds for the terms in the third line of (4.27).

We define the operator $J_{5}$ via its integral kernel

$J_{5}(\mathbf{z},x)=\int_{\mathbb{R}^{3}}\overline{V_{t,l}(\mathbf{z},y)}\overline{J(x,y)}\,\mathrm{d}y$
(4.41)

and write

$\displaystyle\left|\int_{\mathbb{R}^{6}}\phi_{t}(x)J(x,y)\langle\xi_{t},\left(a^{*}(U_{t,\ell,y})+a(\overline{V_{t,\ell,y}})\right)\xi_{t}\rangle\,\mathrm{d}(x,y)\right|$
$\displaystyle=\left|\int_{\mathbb{R}^{3}}\phi_{t}(x)\langle\xi_{t},\left(a^{*}(J_{4,x})+a(J_{5,x})\right)\xi_{t}\rangle\,\mathrm{d}x\right|$
$\displaystyle\leq\|\phi_{t}\|_{2}\
\|J\|\left(\|U_{t}\|^{2}+\|V_{t}\|^{2}\right)^{1/2}\langle\xi_{t},\mathcal{N}\xi_{t}\rangle^{1/2}.$
(4.42)

The other term in the second line of (4.27) can be estimated in the same way.
In combination, (4.27), (4.34), (4.40), and (4.42) prove the claim. ∎

###  4.4 Proof of Theorem 1.2

In the next proposition, whose proof is deferred to the next section, we
control the growth of the expectation of the number of particles operator.

######  Proposition 4.3.

Let $v$ satisfy the same assumptions as in Theorem 1.2, and let the pair
$(\phi,\gamma)$ satisfy Assumption 1.1 with $0<s\leq 3/2$. For every
$k\in\mathbb{N}$, there exist constants $c,C>0$ independent of $N,t$ such that

$\displaystyle\langle\xi_{t},\mathcal{N}\xi_{t}\rangle\leq$
$\displaystyle\exp(c\exp(c\exp(ct)))\langle\xi,\mathcal{N}\xi\rangle$
$\displaystyle+CN^{(22-k)/14}\exp(c\exp(c\exp(ct)))\left(\langle\xi,\mathcal{N}^{4}\xi\rangle^{1/2}+N^{10/7}\right)\big{\langle}\xi,\mathcal{N}^{k+2}\xi\big{\rangle}^{1/2}$
(4.43)

for every $N\in\mathbb{N}$, $t>0$.

######  Remark 4.4.

Similar bounds can be established for higher moments of the number of
particles operator on the fluctuation vector
$\langle\xi_{t},\mathcal{N}^{j}\xi_{t}\rangle$, for any $j\geq 2$. The proof
is essentially the same as for (4.43). We do not pursue this generalization
here because we do not need it for our main result, and we want to keep the
notation as simple as possible.

Next, we show how Proposition 4.3 implies Theorem 1.2.

###### Proof of Theorem 1.2.

We combine (4.29) and Lemma 4.1 to see that

$\|\gamma_{\xi,t}-|\phi_{t}\rangle\langle\phi_{t}|-\gamma_{t}\|_{\mathcal{L}^{1}}\lesssim\exp(ct)\left(T_{\mathrm{c}}(s)\langle\xi_{t},\mathcal{N}\xi_{t}\rangle+\sqrt{NT_{\mathrm{c}}(s)}\langle\xi_{t},\mathcal{N}\xi_{t}\rangle^{1/2}\right)\,.$
(4.44)

An application of Proposition 4.3 with $k=42$ proves

$\langle\xi_{t},\mathcal{N}\xi_{t}\rangle\lesssim\exp(c\exp(c\exp(ct)))\,,$
(4.45)

under the assumptions of Theorem 1.2. Plugging (4.45) into (4.44) and using
the bound

$\left\|\
|\phi_{t}\rangle\langle\phi_{t}|+\gamma_{t}-\left(|\phi_{t}^{\mathrm{H}}\rangle\langle\phi_{t}^{\mathrm{H}}|+e^{\mathrm{i}\Delta
t}\gamma e^{-\mathrm{i}\Delta t}\right)\right\|_{\mathcal{L}^{1}}\lesssim
N^{1/2}T_{c}^{3/4}(s)\exp(ct)\,,$ (4.46)

which follows from (3.6), we get (1.31), concluding the proof of the theorem.
∎

######  Remark 4.5.

Remark 1.3.3 is proven by using the bound (3.7) instead of (4.46) at the end
of the proof of Theorem 1.2.

##  5 Study of the fluctuation dynamics

The aim of this section is to prove Proposition 4.3. We will apply a Grönwall
argument, and therefore start our analysis with the computation of the
generator of the fluctuation dynamics $\mathcal{U}^{\mathrm{fluct}}(t,s)$ in
(4.24).

###  5.1 Generator of the fluctuation dynamics

The propagator of the fluctuations satisfies the equation

$\mathrm{i}\partial_{t}\mathcal{U}^{\mathrm{fluct}}(t,s)=\mathcal{G}_{N,t}\mathcal{U}^{\mathrm{fluct}}(t,s)$
(5.1)

with the time-dependent generator

$\mathcal{G}_{N,t}=\left(\mathrm{i}\partial_{t}\mathcal{T}_{t}^{*}\right)\mathcal{T}_{t}+\mathcal{T}_{t}^{*}\left(\mathrm{i}\partial_{t}\mathcal{W}_{t}^{*}\right)\mathcal{W}_{t}\mathcal{T}_{t}+\mathcal{T}_{t}^{*}\mathcal{W}_{t}^{*}\mathcal{L}_{N}\mathcal{W}_{t}\mathcal{T}_{t}\,,$
(5.2)

where we used the shorthand $\mathcal{W}_{t}=\mathcal{W}(\phi_{t})$. In the
following proposition we compute $\mathcal{G}_{N,t}$. To keep formulas at a
reasonable length, we introduce the notation

$v(\mathbf{x},\mathbf{y})=\begin{cases}v(x-y)&\text{ if
}\sigma=\sigma^{\prime}=\ell\,,\\\ -v(x-y)&\text{ if
}\sigma=\sigma^{\prime}=r\,,\\\ 0&\text{ instead}\,,\end{cases}$ (5.3)

where $\mathbf{x}=(\sigma,x)$ and $\mathbf{y}=(\sigma^{\prime},y)$. We also
recall the notation $\int
f(\mathbf{x})\,\mathrm{d}\mathbf{x}=\sum_{\sigma=l,r}\int
f(\sigma,x)\,\mathrm{d}x$ for
$f:\\{\ell,r\\}\times\mathbb{R}^{3}\to\mathbb{C}$, which was introduced in
Section 4 and will be used throught the present Section.

######  Proposition 5.1.

Let $v$ be an interaction potential satisfying the same assumptions as in
Theorem 1.2, and let $\phi\in H^{3}(\mathbb{R}^{3})$,
$\gamma\in\mathcal{H}^{3,1}$. The generator of the fluctuation dynamics
$\mathcal{U}^{\mathrm{fluct}}(t,s)$ can be written as

$\mathcal{G}_{N,t}=\mathcal{I}_{N,t}^{(1)}+\mathcal{I}_{N,t}^{(2)}+\mathcal{I}_{N,t}^{(3)}+\mathcal{I}_{N,t}^{(4)}\,,$
(5.4)

where

$\displaystyle\mathcal{I}_{N,t}^{(1)}=$ $\displaystyle\;\frac{1}{2N}\int
v(\mathbf{x},\mathbf{y})\Big{(}a^{*}(U_{t,\mathbf{x}})a^{*}(U_{t,\mathbf{y}})a(U_{t,\mathbf{y}})a(U_{t,\mathbf{x}})$
$\displaystyle\;+a^{*}(U_{t,\mathbf{x}})a^{*}(\overline{V_{t,\mathbf{y}}})a(\overline{V_{t,\mathbf{y}}})a(U_{t,\mathbf{x}})+a^{*}(U_{t,\mathbf{y}})a^{*}(\overline{V_{t,\mathbf{x}}})a(\overline{V_{t,\mathbf{x}}})a(U_{t,\mathbf{y}})$
$\displaystyle\;+a^{*}(\overline{V_{t,\mathbf{x}}})a^{*}(\overline{V_{t,\mathbf{y}}})a(\overline{V_{t,\mathbf{y}}})a(\overline{V_{t,\mathbf{x}}})+2a^{*}(U_{t,\mathbf{x}})a^{*}(\overline{V_{t,\mathbf{y}}})a(\overline{V_{t,\mathbf{y}}})a(U_{t,\mathbf{x}})\Big{)}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}$
(5.5)

collects the terms that are quartic in $a$ and $a^{*}$ and commute with
$\mathcal{N}$,

$\displaystyle\mathcal{I}_{N,t}^{(2)}=\;$ $\displaystyle\frac{1}{2N}\int
v(\mathbf{x},\mathbf{y})\Big{(}a^{*}(U_{t,\mathbf{x}})a^{*}(U_{t,\mathbf{y}})a^{*}(\overline{V_{t,\mathbf{y}}})a^{*}(\overline{V_{t,\mathbf{x}}})+a^{*}(U_{t,\mathbf{x}})a^{*}(U_{t,\mathbf{y}})a^{*}(\overline{V_{t,\mathbf{y}}})a(U_{t,\mathbf{x}})$
$\displaystyle\;+a^{*}(U_{t,\mathbf{x}})a^{*}(\overline{V_{t,\mathbf{y}}})a^{*}(\overline{V_{t,\mathbf{x}}})a(\overline{V_{t,\mathbf{y}}})+a^{*}(U_{t,\mathbf{x}})a^{*}(U_{t,\mathbf{y}})a^{*}(\overline{V_{t,\mathbf{x}}})a(U_{t,\mathbf{y}})$
$\displaystyle\;+a^{*}(U_{t,\mathbf{y}})a^{*}(\overline{V_{t,\mathbf{y}}})a^{*}(\overline{V_{t,\mathbf{x}}})a(\overline{V_{t,\mathbf{x}}})+\operatorname{h.c.}\Big{)}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}$
(5.6)

contains the quartic terms that do not commute with $\mathcal{N}$,

$\displaystyle\mathcal{I}_{N,t}^{(3)}=\;$ $\displaystyle\frac{1}{N}\int
v(\mathbf{x},\mathbf{y})\Big{\\{}\boldsymbol{\phi}(\mathbf{x})\Big{(}a^{*}(U_{t,\mathbf{x}})a^{*}(U_{t,\mathbf{y}})a^{*}(\overline{V_{t,\mathbf{y}}})+a^{*}(U_{t,\mathbf{x}})a^{*}(U_{t,\mathbf{y}})a(U_{t,\mathbf{y}})$
$\displaystyle\;+a^{*}(U_{t,\mathbf{x}})a^{*}(\overline{V_{t,\mathbf{y}}})a(\overline{V_{t,\mathbf{y}}})+a^{*}(U_{t,\mathbf{y}})a^{*}(\overline{V_{t,\mathbf{y}}})a(\overline{V_{t,\mathbf{x}}})+a^{*}(\overline{V_{t,\mathbf{y}}})a(\overline{V_{t,\mathbf{x}}})a(\overline{V_{t,\mathbf{y}}})$
$\displaystyle\;+a^{*}(U_{t,\mathbf{x}})a(\overline{V_{t,\mathbf{y}}})a(U_{t,\mathbf{y}})+a^{*}(U_{t,\mathbf{y}})a(\overline{V_{t,\mathbf{x}}})a(U_{t,\mathbf{y}})+a(\overline{V_{t,\mathbf{x}}})a(\overline{V_{t,\mathbf{y}}})a(U_{t,\mathbf{y}})\Big{)}+\operatorname{h.c.}\Big{\\}}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}$
(5.7)

the cubic terms, and

$\mathcal{I}_{N,t}^{(4)}=\operatorname{Tr}(V_{t}LV_{t}^{*})+\frac{1}{2N}\int
v(\mathbf{x},\mathbf{y})|\boldsymbol{\phi}_{t}(\mathbf{x})|^{2}|\boldsymbol{\phi}_{t}(\mathbf{y})|^{2}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}\,,$

with

$L=\begin{pmatrix}-\Delta&0\\\ 0&\Delta\end{pmatrix}\,,$ (5.8)

is a time-dependent constant.

######  Remark 5.2.

Comparing with the literature (especially [7, Proposition 3.1]), it may seem
surprising that the kinetic energy $\,\mathrm{d}\Gamma(\Delta)$ of the
excitations does not appear in (5.4). This is a consequence of a different
definition of the symplectomorphism $T_{t}$ in (4.15), which of course leads
to a different Bogoliubov transformation $\mathcal{T}_{t}$. In particular,
from (4.13), we observe that the HFB propagator $\mathcal{U}_{t}$ only acts
from the right on the initial data (4.14). Effectively, this amounts to
considering the fluctuation dynamics in the interaction picture.

###### Proof.

We start our proof with the computation of
$(\mathrm{i}\partial_{t}\mathcal{T}_{t}^{*})\mathcal{T}_{t}$. Our analysis for
this term follows that in the proof of [7, Proposition 3.1]. We observe that

$\mathcal{T}_{t}^{*}A(f,g)\mathcal{T}_{t}=A(T_{t}(f,g))\qquad$

holds for $f,g\in{\mathfrak{h}}\oplus{\mathfrak{h}}$, with the notation
$A(f,g)=a(f)+a^{*}(\overline{g})$. If we differentiate both sides w.r.t. time
we find

$(\mathrm{i}\partial_{t}\mathcal{T}_{t}^{*})\mathcal{T}_{t}A(T_{t}(f,g))+A(T_{t}(f,g))\mathcal{T}_{t}^{*}(\mathrm{i}\partial_{t}\mathcal{T}_{t})=-A((\mathrm{i}\partial_{t}T_{t})(f,g))\,.$

We use
$0=\mathrm{i}\partial_{t}(\mathcal{T}_{t}^{*}\mathcal{T}_{t})=(\mathrm{i}\partial_{t}\mathcal{T}_{t}^{*})\mathcal{T}_{t}+\mathcal{T}_{t}^{*}(\mathrm{i}\partial_{t}\mathcal{T}_{t})$
to arrive at

$\left[(\mathrm{i}\partial_{t}\mathcal{T}_{t}^{*})\mathcal{T}_{t},A(T_{t}(f,g))\right]=-A((\mathrm{i}\partial_{t}T_{t})(f,g))\,.$
(5.9)

Since $T_{t}$ is a symplectomorphism and (5.9) holds for
$f,g\in{\mathfrak{h}}\oplus{\mathfrak{h}}$, we conclude that
$(\mathrm{i}\partial_{t}\mathcal{T}_{t}^{*})\mathcal{T}_{t}$ is quadratic in
creation and annihilation operators. That is, we can write it as

$(\mathrm{i}\partial_{t}\mathcal{T}_{t}^{*})\mathcal{T}_{t}=\int
C_{t}(\mathbf{x},\mathbf{y})a^{*}_{\mathbf{x}}a_{\mathbf{y}}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}+\frac{1}{2}\left(\int
D_{t}(\mathbf{x},\mathbf{y})a_{\mathbf{x}}^{*}a_{\mathbf{y}}^{*}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}+\operatorname{h.c.}\right)\,,$
(5.10)

where $C_{t}(\mathbf{x},\mathbf{y})$ and $D_{t}(\mathbf{x},\mathbf{y})$ denote
the kernels of two operators $C_{t}$ and $D_{t}$ on
${\mathfrak{h}}\oplus{\mathfrak{h}}$. The operator $C_{t}$ is self-adjoint,
i.e.

$C_{t}(\mathbf{y},\mathbf{x})=\overline{C_{t}(\mathbf{x},\mathbf{y})}\,,$

and we assume

$D(\mathbf{x},\mathbf{y})=D(\mathbf{y},\mathbf{x})\,,$

which will be justified a-posteriori. We use (5.9) to determine the kernels. A
straightforward computation shows that

$\begin{split}\left[\int
C_{t}(\mathbf{x},\mathbf{y})a^{*}_{\mathbf{x}}a_{\mathbf{y}}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y},A(T_{t}(f,0))\right]=&\;-a(C_{t}U_{t}f)+a^{*}(C_{t}\overline{V_{t}}\
\overline{f})\,,\\\ \left[\frac{1}{2}\int
D_{t}(\mathbf{x},\mathbf{y})a^{*}_{\mathbf{x}}a^{*}_{\mathbf{y}}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y},A(T_{t}(f,0))\right]=&\;-a^{*}(D_{t}\overline{U_{t}}\
\overline{f})\,,\\\
\left[\frac{1}{2}\int\overline{D_{t}(\mathbf{x},\mathbf{y})}a_{\mathbf{x}}a_{\mathbf{y}}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y},A(T_{t}(f,0))\right]=&\;a(D_{t}V_{t}f)\,.\end{split}$

We insert this into (5.9) and obtain the following system of equations for
$C_{t},D_{t}$:

$\begin{pmatrix}C_{t}&D_{t}\end{pmatrix}\begin{pmatrix}-U_{t}&\overline{V_{t}}\\\
V_{t}&-\overline{U_{t}}\end{pmatrix}=\begin{pmatrix}-\mathrm{i}\partial_{t}U_{t}&\mathrm{i}\partial_{t}\overline{V_{t}}\end{pmatrix}.$

Using (4.5) and (4.11), we see that

$\begin{split}\begin{pmatrix}C_{t}&D_{t}\end{pmatrix}=&\begin{pmatrix}i\partial_{t}U_{t}&-i\partial_{t}\overline{V_{t}}\end{pmatrix}\begin{pmatrix}U_{t}^{*}&V_{t}^{*}\\\
\overline{V^{*}_{t}}&\overline{U^{*}_{t}}\end{pmatrix}\\\
=&-\begin{pmatrix}U_{t}&\overline{V_{t}}\end{pmatrix}\begin{pmatrix}H_{t}&K_{t}\\\
\overline{K_{t}}&\overline{H_{t}}\end{pmatrix}\begin{pmatrix}U_{t}^{*}&V_{t}^{*}\\\
\overline{V^{*}_{t}}&\overline{U^{*}_{t}}\end{pmatrix}\end{split}$

holds, and hence

$\displaystyle C_{t}=$
$\displaystyle\;-\left(U_{t}H_{t}U_{t}^{*}+U_{t}K_{t}\overline{V^{*}_{t}}+\overline{V_{t}}\
\overline{K_{t}}\ U_{t}^{*}+\overline{V_{t}}\ \overline{H_{t}}\
\overline{V_{t}^{*}}\right)\,,$ $\displaystyle D_{t}=$
$\displaystyle\;-\left(U_{t}H_{t}V^{*}_{t}+U_{t}K_{t}\overline{U_{t}^{*}}+\overline{V_{t}}\
\overline{K_{t}}V_{t}^{*}+\overline{V_{t}}\ \overline{H_{t}}\
\overline{U_{t}^{*}}\right)\,.$ (5.11)

Observe that $D_{t}$ is indeed symmetric, that is,
$\overline{D_{t}^{*}}=D_{t}$.

Next, we compute
$\mathcal{T}_{t}^{*}(\mathrm{i}\partial_{t}\mathcal{W}_{t}^{*})\mathcal{W}_{t}\mathcal{T}_{t}$.
A short computation shows

$(\mathrm{i}\partial_{t}\mathcal{W}_{t}^{*})\mathcal{W}_{t}=-\left(a^{*}(\mathrm{i}\partial_{t}\boldsymbol{\phi}_{t})+a(\mathrm{i}\partial_{t}\boldsymbol{\phi}_{t})\right)\,,$

which implies

$\mathcal{T}_{t}^{*}(\mathrm{i}\partial_{t}\mathcal{W}_{t}^{*})\mathcal{W}_{t}\mathcal{T}_{t}=-\left(a^{*}(U_{t}(\mathrm{i}\partial_{t}\boldsymbol{\phi}_{t}))+a^{*}(\overline{V_{t}}(-\mathrm{i}\partial_{t}\overline{\boldsymbol{\phi}_{t}}))+\operatorname{h.c.}\right)\,.$
(5.12)

We write $\mathcal{L}_{N}=\mathcal{L}_{0}+\mathcal{V}_{N}$ with

$\begin{split}\mathcal{L}_{0}=&\;\,\mathrm{d}\Gamma_{\ell}(-\Delta)-\,\mathrm{d}\Gamma_{r}(-\Delta)=\,\mathrm{d}\Gamma(L)\,,\\\
\mathcal{V}_{N}=&\;\frac{1}{2N}\int
v(x-y)a^{*}_{\ell,x}a^{*}_{\ell,y}a_{\ell,y}a_{\ell,x}\,\mathrm{d}x\,\mathrm{d}y-\frac{1}{2N}\int
v(x-y)a^{*}_{r,x}a^{*}_{r,y}a_{r,y}a_{r,x}\,\mathrm{d}x\,\mathrm{d}y\\\
=&\;\frac{1}{2N}\int
v(\mathbf{x},\mathbf{y})a^{*}_{\mathbf{x}}a^{*}_{\mathbf{y}}a_{\mathbf{y}}a_{\mathbf{x}}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}\,,\end{split}$

and $L$ as defined in (5.8). A straightforward calculation shows that

$\mathcal{W}_{t}^{*}\mathcal{L}_{0}\mathcal{W}_{t}=\mathcal{L}_{0}+a(L\boldsymbol{\phi}_{t})+a^{*}(L\boldsymbol{\phi}_{t})\,.$

Using (4.7) we can compute

$\begin{split}\mathcal{T}_{t}^{*}\mathcal{L}_{0}\mathcal{T}_{t}=&\;\,\mathrm{d}\Gamma(U_{t}LU_{t}^{*}+\overline{V_{t}}L\overline{V_{t}^{*}})+\operatorname{Tr}(V_{t}LV_{t}^{*})\\\
&\;+\int
a_{\mathbf{x}}^{*}a_{\mathbf{y}}^{*}(U_{t}LV_{t}^{*})(\mathbf{x},\mathbf{y})\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}+\operatorname{h.c.}\,,\end{split}$

as well as

$\mathcal{T}_{t}^{*}[a(L\boldsymbol{\phi}_{t})+a^{*}(L\boldsymbol{\phi}_{t})]\mathcal{T}_{t}=a(U_{t}L\boldsymbol{\phi}_{t}+\overline{V_{t}}L\overline{\boldsymbol{\phi}_{t}})+a^{*}(U_{t}L\boldsymbol{\phi}_{t}+\overline{V_{t}}L\overline{\boldsymbol{\phi}_{t}})\,.$

We conclude that

$\displaystyle\mathcal{T}_{t}^{*}\mathcal{W}_{t}^{*}\mathcal{L}_{0}\mathcal{W}_{t}\mathcal{T}_{t}=$
$\displaystyle\;a(U_{t}L\boldsymbol{\phi}_{t}+\overline{V_{t}}L\overline{\boldsymbol{\phi}_{t}})+a^{*}(U_{t}L\boldsymbol{\phi}_{t}+\overline{V_{t}}L\overline{\boldsymbol{\phi}_{t}})+\operatorname{Tr}(V_{t}LV_{t}^{*})$
$\displaystyle\;+\,\mathrm{d}\Gamma(U_{t}LU_{t}^{*}+\overline{V_{t}}L\overline{V_{t}^{*}})+\left(\int
a_{\mathbf{x}}^{*}a_{\mathbf{y}}^{*}(U_{t}LV_{t}^{*})(\mathbf{x},\mathbf{y})\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}+\operatorname{h.c.}\right)\,.$
(5.13)

Let us continue with the computation of
$\mathcal{T}_{t}^{*}\mathcal{W}_{t}^{*}\mathcal{V}_{N}\mathcal{W}_{t}\mathcal{T}_{t}$.
Using (4.2), we obtain

$\begin{split}\mathcal{W}_{t}^{*}a_{\mathbf{x}}^{*}a_{\mathbf{y}}^{*}a_{\mathbf{y}}a_{\mathbf{x}}\mathcal{W}_{t}=&\;a_{\mathbf{x}}^{*}a_{\mathbf{y}}^{*}a_{\mathbf{y}}a_{\mathbf{x}}+\left(\boldsymbol{\phi}_{t}(\mathbf{x})a_{\mathbf{x}}^{*}a_{\mathbf{y}}^{*}a_{\mathbf{y}}+\overline{\boldsymbol{\phi}_{t}(\mathbf{x})}a_{\mathbf{y}}^{*}a_{\mathbf{y}}a_{\mathbf{x}}+(\mathbf{x}\longleftrightarrow\mathbf{y})\right)\\\
&\;+\left(|\boldsymbol{\phi}_{t}(\mathbf{y})|^{2}a_{\mathbf{x}}^{*}a_{\mathbf{x}}+\boldsymbol{\phi}_{t}(\mathbf{x})\overline{\boldsymbol{\phi}_{t}(\mathbf{y})}a_{\mathbf{x}}^{*}a_{\mathbf{y}}+(\mathbf{x}\longleftrightarrow\mathbf{y})\right)+\left(\boldsymbol{\phi}_{t}(\mathbf{x})\boldsymbol{\phi}_{t}(\mathbf{y})a_{\mathbf{x}}^{*}a^{*}_{\mathbf{y}}+\operatorname{h.c.}\right)\\\
&\;+\left(a_{\mathbf{x}}^{*}\boldsymbol{\phi}_{t}(\mathbf{x})|\boldsymbol{\phi}_{t}(\mathbf{y})|^{2}+a_{\mathbf{x}}\overline{\boldsymbol{\phi}_{t}(\mathbf{x})}|\boldsymbol{\phi}_{t}(\mathbf{y})|^{2}+(\mathbf{x}\longleftrightarrow\mathbf{y})\right)+|\boldsymbol{\phi}_{t}(\mathbf{x})|^{2}|\boldsymbol{\phi}_{t}(\mathbf{y})|^{2}\,,\end{split}$

where $(\mathbf{x}\longleftrightarrow\mathbf{y})$ is a shorthand notation for
the expression preceding it with $\mathbf{x}$ and $\mathbf{y}$ exchanged. We
conjugate both sides with $\mathcal{T}_{t}$, apply (4.7), bring all terms to
normal order and integrate the result against the potential
$\frac{1}{2N}v(\mathbf{x},\mathbf{y})$. We omit the details of this lenghty
but straightforward computation, which leads us to

$\mathcal{T}_{t}^{*}\mathcal{W}_{t}^{*}\mathcal{V}_{N}\mathcal{W}_{t}\mathcal{T}_{t}=\sum_{i=1}^{3}\mathcal{I}_{N,t}^{(i)}+\mathcal{Q}+\mathcal{L}+C\,,$

where the operators $\mathcal{I}_{N,t}^{(i)}$ are defined below (5.4).
Moreover,

$\begin{split}\mathcal{Q}=&\;\frac{1}{2N}\int
v(\mathbf{x},\mathbf{y})\Bigg{\\{}\boldsymbol{\phi}_{t}(\mathbf{x})\boldsymbol{\phi}_{t}(\mathbf{y})\Big{[}a^{*}(U_{t,\mathbf{x}})a^{*}(U_{t,\mathbf{y}})+a^{*}(U_{t,\mathbf{x}})a(\overline{V_{t,\mathbf{y}}})\\\
&\qquad\qquad\qquad\qquad\qquad+a^{*}(U_{t,\mathbf{y}})a(\overline{V_{t,\mathbf{x}}})+a(\overline{V_{t,\mathbf{x}}})a(\overline{V_{t,\mathbf{y}}})\Big{]}+\operatorname{h.c.}\Bigg{\\}}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}\\\
&\;+\frac{1}{N}\int
v(\mathbf{x},\mathbf{y})|\boldsymbol{\phi}_{t}(\mathbf{y})|^{2}\Big{[}a^{*}(U_{t,\mathbf{x}})a(U_{t,\mathbf{x}})+a^{*}(U_{t,\mathbf{x}})a^{*}(\overline{V_{t,\mathbf{x}}})+a(\overline{V_{t,\mathbf{x}}})a(U_{t,\mathbf{x}})+a^{*}(\overline{V_{t,\mathbf{x}}})a(\overline{V_{t,\mathbf{x}}})\Big{]}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}\\\
&\;+\frac{1}{N}\int
v(\mathbf{x},\mathbf{y})\boldsymbol{\phi}_{t}(\mathbf{x})\overline{\boldsymbol{\phi}_{t}(\mathbf{y})}\Big{[}a^{*}(U_{t,\mathbf{x}})a(U_{t,\mathbf{y}})+a^{*}(U_{t,\mathbf{x}})a^{*}(\overline{V_{t,\mathbf{y}}})+a(\overline{V_{t,\mathbf{x}}})a(U_{t,\mathbf{y}})+a^{*}(\overline{V_{t,\mathbf{y}}})a(\overline{V_{t,\mathbf{x}}})\Big{]}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}\\\
&\;+\frac{1}{2N}\int
v(\mathbf{x},\mathbf{y})\Bigg{\\{}\left(U_{t,\mathbf{y}},\overline{V_{t,\mathbf{x}}}\right)\Big{[}a^{*}(U_{t,\mathbf{y}})a(\overline{V_{t,\mathbf{x}}})+a^{*}(U_{t,\mathbf{x}})a^{*}(U_{t,\mathbf{y}})\\\
&\qquad\qquad\qquad\qquad\qquad+a(\overline{V_{t,\mathbf{x}}})a(\overline{V_{t,\mathbf{y}}})+a^{*}(U_{t,\mathbf{x}})a(\overline{V_{t,\mathbf{y}}})\Big{]}+\operatorname{h.c.}\Bigg{\\}}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}\\\
&\;+\frac{1}{N}\int
v(\mathbf{x},\mathbf{y})\left(\overline{V_{t,\mathbf{y}}},\overline{V_{t,\mathbf{y}}}\right)\Big{[}a^{*}(U_{t,\mathbf{x}})a(U_{t,\mathbf{x}})+a^{*}(\overline{V_{t,\mathbf{x}}})a(\overline{V_{t,\mathbf{x}}})+a^{*}(U_{t,\mathbf{x}})a^{*}(\overline{V_{t,\mathbf{x}}})+a(\overline{V_{t,\mathbf{x}}})a(U_{t,\mathbf{x}})\Big{]}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}\\\
&\;+\frac{1}{N}\int
v(\mathbf{x},\mathbf{y})\left(\overline{V_{t,\mathbf{y}}},\overline{V_{t,\mathbf{x}}}\right)\Big{[}a^{*}(U_{t,\mathbf{x}})a^{*}(\overline{V_{t,\mathbf{y}}})+a(\overline{V_{t,\mathbf{x}}})a(U_{t,\mathbf{y}})+a^{*}(\overline{V_{t,\mathbf{y}}})a(\overline{V_{t,\mathbf{x}}})+a^{*}(U_{t,\mathbf{x}})a(U_{t,\mathbf{y}})\Big{]}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}\end{split}$

contains the quadratic terms,

$\displaystyle\mathcal{L}=\frac{1}{N}\int$ $\displaystyle
v(\mathbf{x},\mathbf{y})\Bigg{\\{}\boldsymbol{\phi}_{t}(\mathbf{x})\Big{[}|\boldsymbol{\phi}_{t}(\mathbf{y})|^{2}+\left(\overline{V_{t,\mathbf{y}}},\overline{V_{t,\mathbf{y}}}\right)\Big{]}\Big{(}a^{*}(U_{t,\mathbf{x}})+a(\overline{V_{t,\mathbf{x}}})\Big{)}$
$\displaystyle\;+\left(\overline{V_{t,\mathbf{x}}},U_{t,\mathbf{y}}\right)\Big{(}a(U_{t,\mathbf{y}})+a^{*}(\overline{V_{t,\mathbf{y}}})\Big{)}+\left(\overline{V_{t,\mathbf{x}}},\overline{V_{t,\mathbf{y}}}\right)\Big{(}a^{*}(U_{t,\mathbf{y}})+a(\overline{V_{t,\mathbf{y}}})\Big{)}\Big{]}+\operatorname{h.c.}\Bigg{\\}}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}$
(5.14)

the linear terms, and

$C=\frac{1}{2N}\int
v(\mathbf{x},\mathbf{y})|\boldsymbol{\phi}_{t}(\mathbf{x})|^{2}|\boldsymbol{\phi}_{t}(\mathbf{y})|^{2}\,\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{y}$

is a time-dependent constant. Here we used the notation $(\cdot,\cdot)$ for
the inner product in $(L^{2}(\mathbb{R}^{3})\oplus
L^{2}(\mathbb{R}^{3}),\,\mathrm{d}\mathbf{z})$, that is,

$(f,g)=\int\overline{f(\mathbf{z})}g(\mathbf{z})\,\mathrm{d}\mathbf{z}\,.$
(5.15)

Eq. (4.9) implies

$(\overline{V_{t,\mathbf{y}}},\overline{V_{t,\mathbf{x}}})=V_{t}^{*}V_{t}(\mathbf{x},\mathbf{y})\eqqcolon\widetilde{\gamma}_{t}(\mathbf{x},\mathbf{y})\quad\text{
and
}\quad(U_{t,\mathbf{y}},\overline{V_{t,\mathbf{x}}})\eqqcolon\widetilde{\alpha}_{t}(\mathbf{x},\mathbf{y})\,,$
(5.16)

where, $\widetilde{\gamma}_{t}(\mathbf{x},\mathbf{y})$ and
$\widetilde{\alpha}_{t}(\mathbf{x},\mathbf{y})$ denote the kernels of the full
1-pdm
$\widetilde{\gamma}_{t}\in\mathcal{B}({\mathfrak{h}}\oplus{\mathfrak{h}})$ and
of the full pairing function
$\widetilde{\alpha}_{t}\in\mathcal{B}({\mathfrak{h}}\oplus{\mathfrak{h}})$ of
the vector state $\Psi_{t}=\mathcal{T}_{t}\Omega$ on
$\mathscr{F}({\mathfrak{h}}\oplus{\mathfrak{h}})$, respectively. It should be
highlighted that $\widetilde{\gamma}_{t}((x,\ell),(y,\ell))=\gamma_{t}(x,y)$
and $\widetilde{\alpha}_{t}((x,\ell),(y,\ell))=\alpha_{t}(x,y)$. We use
$\overline{\widetilde{\gamma}_{t}(\mathbf{x},\mathbf{y})}=\widetilde{\gamma}_{t}(\mathbf{y},\mathbf{x})$
and