Content-Type: multipart/mixed; boundary="-------------0810240337487"

This is a multi-part message in MIME format.

---------------0810240337487
Content-Type: text/plain; name="08-198.keywords"
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment; filename="08-198.keywords"

Hartree-Fock, fermions, mean-field limit, Egorov theorem, Coulomb interaction
---------------0810240337487
Content-Type: application/x-tex; name="notes.tex"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline; filename="notes.tex"

\documentclass[a4paper]{article}


\setlength{\textwidth}{14.5cm}
\setlength{\oddsidemargin}{1cm}
\setlength{\evensidemargin}{0.5cm}
\setlength{\topmargin}{-1cm}
\setlength{\textheight}{23.5cm}
\usepackage{cite}

\usepackage{amsmath} %[intlimits]
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{amsxtra}
\usepackage{bbm}
\usepackage{mathrsfs}
\usepackage{ifthen}

\setlength{\unitlength}{1cm}

\newcommand{\ul}[1]{\underline{#1} \!\,} %underline
\newcommand{\ol}[1]{\overline{#1} \!\,} %overline
\newcommand{\op}[1]{\Hat{#1}} %operator
\newcommand{\vecop}[1]{\Hat{\vec{#1}}} %vector operator
\newcommand{\define}{\textbf}

%special symbols
\newcommand{\mi}{\mathrm{i}}
\newcommand{\md}{\mathrm{d}}
\newcommand{\me}{\mathrm{e}}
\newcommand{\ee}{\mathrm{e}}
\newcommand{\ii}{\mathrm{i}}
\newcommand{\dd}{\mathrm{d}}
\newcommand{\ddd}{\dd \mspace{-4mu} \bar{\,} \mspace{1mu}}

\newcommand{\deq}{\mathrel{\mathop:}=}
\newcommand{\eqd}{=\mathrel{\mathop:}}
\newcommand{\id}{\mspace{2mu}\mathrm{i}\mspace{-0.6mu}\mathrm{d}} %identity map
\newcommand{\grad}{\vec{\nabla}}
\newcommand{\s}{\mspace{-0.4mu}}
\newcommand{\umat}{\mathbbmss{1}} %unit matrix
\renewcommand{\epsilon}{\varepsilon}

\newcommand{\EE} {\mathrm{E}}
\newcommand{\PP} {\mathrm{P}}
\newcommand{\HH} {\mathrm{H}}
\newcommand{\LL}{\mathrm{L}}


\newcommand{\floor}[1] {\lfloor {#1} \rfloor}
\newcommand{\ceil}[1]  {\lceil  {#1} \rceil}
\newcommand{\genarg} {\,\cdot\,}

\newcommand{\R} {\mathbb{R}}
\newcommand{\C} {\mathbb{C}}
\newcommand{\N} {\mathbb{N}}
\newcommand{\Z} {\mathbb{Z}}
\newcommand{\Q} {\mathbb{Q}}

%%%%%%%%%%%%%%%%%%%%%%%   Paranthesis   %%%%%%%%%%%%%%%%%%%%%
\newcommand{\p}[1]{({#1})}
\newcommand{\pb}[1]{\bigl({#1}\bigr)}
\newcommand{\pB}[1]{\Bigl({#1}\Bigr)}
\newcommand{\pbb}[1]{\biggl({#1}\biggr)}
\newcommand{\pBB}[1]{\Biggl({#1}\Biggr)}
\newcommand{\pa}[1]{\left({#1}\right)}

\newcommand{\q}[1]{[{#1}]}
\newcommand{\qb}[1]{\bigl[{#1}\bigr]}
\newcommand{\qB}[1]{\Bigl[{#1}\Bigr]}
\newcommand{\qbb}[1]{\biggl[{#1}\biggr]}
\newcommand{\qBB}[1]{\Biggl[{#1}\Biggr]}
\newcommand{\qa}[1]{\left[{#1}\right]}

\newcommand{\h}[1]{\{{#1}\}}
\newcommand{\hb}[1]{\bigl\{{#1}\bigr\}}
\newcommand{\hB}[1]{\Bigl\{{#1}\Bigr\}}
\newcommand{\hbb}[1]{\biggl\{{#1}\biggr\}}
\newcommand{\hBB}[1]{\Biggl\{{#1}\Biggr\}}
\newcommand{\ha}[1]{\left\{{#1}\right\}}

%absolute value
\newcommand{\abs}[1]{\lvert #1 \rvert}
\newcommand{\absb}[1]{\big\lvert #1 \big\rvert}
\newcommand{\absB}[1]{\Big\lvert #1 \Big\rvert}
\newcommand{\absbb}[1]{\bigg\lvert #1 \bigg\rvert}
\newcommand{\absBB}[1]{\Bigg\lvert #1 \Bigg\rvert}
\newcommand{\absa}[1]{\left\lvert #1 \right\rvert}

%norm
\newcommand{\norm}[1]{\lVert #1 \rVert}
\newcommand{\normb}[1]{\big\lVert #1 \big\rVert}
\newcommand{\normB}[1]{\Big\lVert #1 \Big\rVert}
\newcommand{\normbb}[1]{\bigg\lVert #1 \bigg\rVert}
\newcommand{\normBB}[1]{\Bigg\lVert #1 \Bigg\rVert}
\newcommand{\norma}[1]{\left\lVert #1 \right\rVert}

%average
\newcommand{\avg}[1]{\langle #1 \rangle}
\newcommand{\avgb}[1]{\big\langle #1 \big\rangle}
\newcommand{\avgB}[1]{\Big\langle #1 \Big\rangle}
\newcommand{\avgbb}[1]{\bigg\langle #1 \bigg\rangle}
\newcommand{\avgBB}[1]{\Bigg\langle #1 \Bigg\rangle}
\newcommand{\avga}[1]{\left\langle #1 \right\rangle}

%scalar product with angle brackets
\newcommand{\scalar}[2]{\langle{#1} \mspace{2mu}, {#2}\rangle}
\newcommand{\scalarb}[2]{\big\langle{#1} \mspace{2mu}, {#2}\big\rangle}
\newcommand{\scalarB}[2]{\Big\langle{#1} \,\mspace{2mu},\, {#2}\Big\rangle}
\newcommand{\scalarbb}[2]{\bigg\langle{#1} \,\mspace{2mu},\, {#2}\bigg\rangle}
\newcommand{\scalarBB}[2]{\Bigg\langle{#1} \,\mspace{2mu},\, {#2}\Bigg\rangle}
\newcommand{\scalara}[2]{\left\langle{#1} \,\mspace{2mu},\, {#2}\right\rangle}

%Dirac bra-c-ket notation
\newcommand{\bra}[1]{\langle #1 |}
\newcommand{\brab}[1]{\big\langle #1 \big|}
\newcommand{\braB}[1]{\Big\langle #1 \Big|}
\newcommand{\brabb}[1]{\bigg\langle #1 \bigg|}
\newcommand{\braBB}[1]{\Bigg\langle #1 \Bigg|}

\newcommand{\ket}[1]{| #1 \rangle}
\newcommand{\ketb}[1]{\big| #1 \big\rangle}
\newcommand{\ketB}[1]{\Big| #1 \Big\rangle}
\newcommand{\ketbb}[1]{\bigg| #1 \bigg\rangle}
\newcommand{\ketBB}[1]{\Bigg| #1 \Bigg\rangle}

\newcommand{\dscalar}[2]{\langle #1 \,|\, #2 \rangle}
\newcommand{\dscalarb}[2]{\big\langle #1 \,\big|\, #2 \big\rangle}
\newcommand{\dscalarB}[2]{\Big\langle #1 \,\Big|\, #2 \Big\rangle}
\newcommand{\dscalarbb}[2]{\bigg\langle #1 \,\bigg|\, #2 \bigg\rangle}
\newcommand{\dscalarBB}[2]{\Bigg\langle #1 \,\Bigg|\, #2 \Bigg\rangle}

\newcommand{\dexp}[3]{\langle #1 \,|\, #2 \,|\, #3 \rangle}
\newcommand{\dexpb}[3]{\big\langle #1 \,\big|\, #2 \,\big|\, #3 \big\rangle}
\newcommand{\dexpB}[3]{\Big\langle #1 \,\Big|\, #2 \,\Big|\, #3 \Big\rangle}
\newcommand{\dexpbb}[3]{\bigg\langle #1 \,\bigg|\, #2 \,\bigg|\, #3 \bigg\rangle}
\newcommand{\dexpBB}[3]{\Bigg\langle #1 \,\Bigg|\, #2 \,\Bigg|\, #3 \Bigg\\rangle}

%Operators
\DeclareMathOperator*{\slim}{s-lim}
\DeclareMathOperator*{\wlim}{w-lim}
\DeclareMathOperator*{\wstarlim}{w*-lim}
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator{\tr}{Tr}

%\DeclareMathOperator{\Tr}{Tr}
\DeclareMathOperator{\var}{var}
\DeclareMathOperator{\supp}{supp}
\DeclareMathOperator{\cotanh}{cotanh}
\DeclareMathOperator{\U}{U}
\DeclareMathOperator{\SU}{SU}
\DeclareMathOperator{\SO}{SO}
\DeclareMathOperator{\su}{su}
\DeclareMathOperator{\re}{Re}
\DeclareMathOperator{\im}{Im}
\DeclareMathOperator{\dom}{\mathcal{D}}
\DeclareMathOperator{\ran}{\mathcal{R}}
\DeclareMathOperator{\keroperator}{\mathcal{N}}
\DeclareMathOperator*{\esssup}{ess\,sup}
\DeclareMathOperator{\ad}{ad}
\DeclareMathOperator{\Ad}{Ad}
\DeclareMathOperator{\sgn}{sgn}
\renewcommand{\ker}{\keroperator}

\renewcommand{\leq}{\leqslant}
\renewcommand{\geq}{\geqslant}

\theoremstyle{definition} %plain, definition, remark
\newtheorem{definition}{Definition}[section]
\newtheorem*{definition*}{Definition}
\newtheorem{example}{Example}[section]
\newtheorem*{example*}{Example}
\newtheorem{remark}{Remark}[section]
\newtheorem*{remarks}{Remarks}
\newtheorem*{remark*}{Remark}
\newtheorem*{notation}{Notation}
\theoremstyle{plain} %plain, definition, remark
\newtheorem{theorem}{Theorem}[section]
\newtheorem*{theorem*}{Theorem}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem*{lemma*}{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem*{corollary*}{Corollary}




\newcommand{\quant}{^{\!\widehat{\;\;\;\;}}_N}
\newcommand{\A}{\mathrm{A}}
\newcommand{\B}{\mathrm{B}}
\renewcommand{\p}[1]{\left(#1 \right)}
\renewcommand{\q}[1]{\left[#1 \right]}
\renewcommand{\h}[1]{\left\{#1 \right\}}

\begin{document}
\title{A Microscopic Derivation of the Time-Dependent Hartree-Fock Equation with Coulomb Two-Body Interaction}
\author{
\\
\\
Antti Knowles \qquad J\"urg Fr\"ohlich
\\
\\
Institute of Theoretical Physics
\\
ETH H\"onggerberg,
\\
CH-8093 Z\"urich, 
\\
Switzerland.
}
\maketitle

\begin{abstract}
We study the dynamics of a Fermi gas with a Coulomb interaction potential, and show that, in a mean-field limiting regime, the dynamics is described by the Hartree-Fock equation. This extends previous work of Bardos et al.\ \cite{BardosGolse} to the case of unbounded interaction potentials. We also express the mean-field limit as a ``superhamiltonian'' system, and state our main result in terms of a Heisenberg-picture dynamics of observables. This is a Egorov-type theorem.
\end{abstract}

\section{Introduction}
The Hartree-Fock equation is a fundamental tool, used throughout physics and chemistry, for describing a system consisting of a large number of fermions. Despite its importance for both conceptual and numerical applications, many questions surrounding it remain unsolved.  One area in which significant progress has been made is the microscopic justification of the static Hartree-Fock equation, which is known to yield the correct asymptotic ground state energy of large atoms and molecules; see \cite{LiebSimon1, LiebSimon2, Bach, Fefferman1, Fefferman2, GrafSolovej}. The time-dependent Hartree-Fock equation, which is supposed to describe the dynamics of a large Fermi system, has received less attention. To our knowledge, the only work in which this equation is derived from microscopic Hamiltonian dynamics is \cite{BardosGolse}. The Cauchy problem for the time-dependent Hartree-Fock equation has also been studied in the literature; see \cite{Bove, Chadam} and especially \cite{Zagatti}, where the Cauchy problem is solved for singular interaction potentials.

A key assumption in \cite{BardosGolse} is that the interaction potential be bounded. A goal of this article is to extend the result of \cite{BardosGolse} to a class of singular interaction potentials, which includes the physically relevant Coulomb potential. We also describe how this mean-field result can be formulated as a Egorov-type theorem.

A system of $N$ fermions is described by a wave function $\psi_N(x_1, \dots, x_N) \in \bigwedge^N \LL^2(\R^3, \dd x)$ which is totally antisymmetric in its arguments. The dynamics of $\psi_N$ is governed by the usual Schr\"odinger equation. In order to obtain a mean-field limit, the Schr\"odinger equation is rescaled with $N$. In this article we adopt the scaling of \cite{BardosGolse}. The Schr\"odinger equation reads
\begin{equation} \label{schroedinger equation}
\ii \partial_t \psi_N(t) \;=\; H_N \psi_N(t)\,,
\end{equation}
where the $N$-particle Hamiltonian $H_N$ is defined by
\begin{equation} \label{mean-field Hamiltonian}
H_N \;\deq\; \sum_{i = 1}^N h_i + \frac{1}{N} \sum_{i < j} w(x_i - x_j)\,.
\end{equation}
Here, $h_i$ is a one-particle Hamiltonian acting on the $i$'th particle, typically of the form $h_i = -\Delta_i + v(x_i)$, where $\Delta$ is the three-dimensional Laplacian and $v$ is some external potential; $w$ is the interaction potential. Under the assumptions on $v$ and $w$ we make below, it is easy to see that $H_N$ is a well-defined self-adjoint operator with domain $\bigwedge^N \HH^2(\R^3)$.

We briefly sketch our main result. Consider a sequence of $N$ orthonormal orbitals $\varphi_1, \dots, \varphi_N$, where $\varphi_i$ is a one-particle wave function. This defines an $N$-particle fermionic state through the Slater determinant
\begin{equation*}
\psi_N \;\deq\; \varphi_1 \wedge \cdots \wedge \varphi_N\,.
\end{equation*}
Let $\psi_N(t)$ be the solution of the Schr\"odinger equation \eqref{schroedinger equation} with initial state $\psi_N$.
In general, $\psi_N(t)$ is no longer a Slater determinant for $t \neq 0$.
However, one expects that this holds asymptotically for large $N$:
\begin{equation*}
\psi_N(t) \;\approx\; \varphi_1(t) \wedge \cdots \wedge \varphi_N(t)\,.
\end{equation*}
Here the orbitals $\varphi_1(t), \dots, \varphi_N(t)$ are supposed to solve the Hartree-Fock equation
\begin{equation} \label{rescaled Hartree-Fock}
\ii \partial_t \varphi_i \;=\; h \varphi_i + \frac{1}{N} \sum_{j  = 1}^N (w * \abs{\varphi_j}^2) \varphi_i - \frac{1}{N} \sum_{j = 1}^N \p{w * (\varphi_i \bar{\varphi}_j)} \varphi_j\,.
\end{equation}
Our main result (Theorem \ref{theorem: main result} below) is a precise formulation of this asymptotic behaviour.

This result is of some physical relevance for studying the dynamics of excited states of electrons in large atoms or molecules in the Born-Oppenheimer approximation. Consider a molecule consisting of $K$ nuclei at fixed positions $R_1, \dots, R_K \in \R^3$, as well as $N$ electrons. The Hamiltonian is given by
\begin{equation*}
\sum_{i = 1}^N \pBB{-\Delta_i - \sum_{k = 1}^K \frac{e_N^2 N z_k}{\abs{x_i - R_k}}} + \sum_{1 \leq i < j \leq N} \frac{e_N^2}{\abs{x_i - x_j}}\,.
\end{equation*}
Here, $e_N$ is the elementary electric charge which we rescale with $N$. The electric charge of nucleus $k$ is $e_N N z_k$, where $z_1, \dots, z_K$ are constants chosen so that $\sum_{k = 1}^K z_k = 1$. This means that the molecule is electrically neutral. If we choose $e_N = e_0 /\sqrt{N}$, for some fixed $e_0$, the Hamiltonian becomes
\begin{equation} \label{Hamiltonian of large molecule}
\sum_{i = 1}^N \pBB{-\Delta_i - \sum_{k = 1}^K \frac{e_0^2 z_k}{\abs{x_i - R_k}}} + \frac{1}{N} \sum_{1 \leq i < j \leq N} \frac{e_0^2}{\abs{x_i - x_j}}\,.
\end{equation}
The scaling of the elementary electric charge $e_N$ may be justified by the fact that the fine structure constant $\alpha = e_N^2 = e_0^2 / N$ is small, i.e.\ $N = O(\alpha^{-1})$. In fact, $\alpha \approx 1/137$.

One problem in the above model, as well as in the works \cite{LiebSimon1, LiebSimon2, Bach, Fefferman1, Fefferman2, GrafSolovej}, is that, as $N$ becomes large, relativistic effects should be taken into account. Indeed, a simple argument shows that the average speed of the innermost electron of an atom with atomic number $Z$ behaves like $Z \alpha$ (in units where the speed of light $c = 1$). Another problem in applying the time-dependent Hartree-Fock theory to the dynamics of excited states is that the interaction with the radiation field is neglected. This interaction is responsible for the relaxation of excited states to the ground state of the molecule.


A physical scenario that is quite different from the large atom or molecule described above is an interacting Fermi-gas confined to a box of fixed size. As discussed in \cite{NarnhoferSewell, ElgartErdosSchleinYau}, the natural scaling in this situation may be viewed as a combination of mean-field and semiclassical scalings. This problem was first studied in \cite{NarnhoferSewell, Spohn}. The authors show that the limiting dynamics is governed by the Vlasov equation. These results were somewhat sharpened in \cite{ElgartErdosSchleinYau}, where the authors compare the Hamiltonian dynamics with the dynamics of the Hartree equation, and derive estimates on the rate of convergence.

Finally, we outline the key ideas of our proof. It relies on the diagrammatic Schwinger-Dyson expansion and Kato smoothing estimates developed in \cite{bosons}. The main steps are:
\begin{itemize}
\item[(a)]
Use the Schwinger-Dyson expansion to express the Hamiltonian time evolution of a $p$-particle observable. 
\item[(b)]
Show that, in the limit $N \to \infty$, only the tree terms of the Schwinger-Dyson expansion survive.
\item[(c)]
Show that the time evolution of a $p$-particle observable under the Hartree-Fock equation converges to the tree terms of the Schwinger-Dyson series as $N \to \infty$.
\end{itemize}
Steps (a) and (b) have been addressed in \cite{bosons}. Thus, the proof below consists in doing step (c).

The article is organised as follows. In Section \ref{section: Hartree-Fock equation} we introduce the Hartree-Fock equation, discuss its Hamiltonian structure and prove a Schwinger-Dyson series for its time evolution. In Section \ref{section: density matrix Hartree-Fock} we rewrite the Hartree-Fock equation using density matrices. Section \ref{section: Slater determinants} is devoted to a discussion of the key properties of Slater determinants. After these preparations, we state and prove our main result in Section \ref{section: convergence}. The final Section \ref{section: Egorov} is devoted to a Egorov-type theorem, which describes the microscopic dynamics as a quantisation of a classical ``superhamiltonian'' theory.

\subsubsection*{Conventions}
In the following, the expression ``$A(t)$ holds for small times'' is understood to mean that there is a constant $T$ such that $A(t)$ is true for all $\abs{t} < T$. The precise value of $T$ can always be inferred from the context. To simplify notation, we assume in the following that $t \geq 0$.

The norm of a Hilbert space $\mathcal{H}$ is denoted by $\norm{\cdot}$. We denote by
\begin{equation*}
\mathcal{H}^{(n)}_\pm \;\deq\; P_\pm \mathcal{H}^{\otimes n}
\end{equation*}
the symmetric/antisymmetric subspaces of the tensor product space $\mathcal{H}^{\otimes n}$. Here, $P_\pm$ is the orthogonal projector onto the symmetric/antisymmetric subspace. The Banach space of bounded operators on $\mathcal{H}$ with operator norm is denoted by $(\mathcal{B}(\mathcal{H}), \norm{\cdot})$, and the Banach space of trace-class operators on $\mathcal{H}$ with trace norm is denoted by $(\mathcal{L}^1(\mathcal{H}), \norm{\cdot}_1)$.

We use the notation $a^{(p)}_{i_1 \dots i_p} \in \mathcal{B}(\mathcal{H}^{\otimes n})$ to denote a $p$-particle operator $a^{(p)} \in \mathcal{B}(\mathcal{H}^{\otimes p})$ acting on the particles $i_1, \dots, i_p \in \{1, \dots, n\}$ in $n$-particle space. Similarly $\tr_{i_1 \dots i_p}$ denotes a partial trace over the degrees of freedom of particles $i_1, \dots, i_p$.

A time subscript of the form $(\cdot)_t$ is always understood to mean time evolution up to time $t$ of $(\cdot)$ with respect to the appropriate free dynamics. We shall explain this in greater detail whenever this notation is used.

The symbol $C$ is reserved for a constant whose dependence on some parameters may be indicated. The value of $C$ need not be the same from one equation to the next.

\section{The Hartree-Fock equation} \label{section: Hartree-Fock equation}
For simplicity of notation, we only consider spinless fermions in the following; the one-particle Hilbert space is $\mathcal{H} \deq \LL^2(\R^3, \dd x) \equiv \LL^2(\R^3)$. Merely cosmetic modifications extend our results to the case of spin-$s$ fermions for which the one-particle Hilbert space is $\LL^2(\R^3) \otimes \C^{2s + 1}$. To fix ideas, we consider the free Hamiltonian $h \deq - \Delta$ and a Coulomb two-body interaction potential $w(x) \deq \abs{x}^{-1}$. By a simple extension of the results of \cite{bosons}, Section 8, our results remain valid for a free Hamiltonian of the form $h = -\Delta + v$ and a two-body interaction potential $w$, where $w$ is even and $v,w \in \LL^\infty(\R^3) + \LL^3_w(\R^3)$ are both real. Here, $\LL^p_w$ denotes the weak $\LL^p$-space (see e.g.\ \cite{ReedSimonII}). In particular, we may treat Hamiltonians of the form \eqref{Hamiltonian of large molecule}.

\subsection{Some notation}
It is convenient to state the time-dependent Hartree-Fock equation in terms of an infinite sequence of orbitals $\Psi = (\psi_i)_{i \in \N}$ which is an element of the Hilbert space
\begin{equation*}
\tilde{\mathcal{H}} \;\deq\; l^2(\N; \LL^2(\R^3)) \;=\; l^2(\N) \otimes \LL^2(\R^3)\,.
\end{equation*}
To simplify notation, we set $\alpha = (x,i)$ and write $\Psi(\alpha) = \psi_i(x)$. Furthermore, we abbreviate
\begin{equation*}
\int \dd \alpha \;\deq\; \sum_{i \in \N} \int \dd x\,, \qquad
\delta(\alpha - \alpha') \;\deq\; \delta_{i i'} \delta(x - x')\,.
\end{equation*}
The scalar product on $\tilde{\mathcal{H}}$ is then given by
\begin{equation*}
\scalar{\Psi}{\Psi'} \;=\; \int \dd \alpha \; \bar{\Psi}(\alpha) \Psi'(\alpha)\,.
\end{equation*}
Let $a^{(p)} \in \mathcal{B}(\mathcal{H}^{\otimes p})$ and define $\tilde{a}^{(p)} \in \mathcal{B}(\tilde{\mathcal{H}}^{\otimes p})$ through
\begin{equation*}
\tilde{a}^{(p)} \;\deq\; \umat_{(l^2(\N))^{\otimes p}} \otimes a^{(p)}\,.
\end{equation*}
We have the identity 
\begin{equation} \label{norm of tilde operator}
\norm{\tilde{a}^{(p)}} \;=\; \norm{a^{(p)}}\,.
\end{equation}
Furthermore,
\begin{equation} \label{extended operator calculated}
\scalarb{\Psi^{\otimes p}}{\tilde{a}^{(p)} \Psi^{\otimes p}} \;=\; \sum_{i_1, \dots, i_p \in \N} \scalarb{\psi_{i_1} \otimes \cdots \otimes \psi_{i_p}}{a^{(p)} \, \psi_{i_1} \otimes \cdots \otimes \psi_{i_p}}\,.
\end{equation}

\subsection{Hamiltonian formulation of the Hartree-Fock equation}
The time-dependent Hartree-Fock equation for $\Psi$ reads
\begin{equation} \label{Hartree-Fock}
\ii \partial_t \psi_i \;=\; h \psi_i + \sum_{j \in \N} (w * \abs{\psi_j}^2) \psi_i - \sum_{j \in \N} \p{w * (\psi_i \bar{\psi}_j)} \psi_j\,.
\end{equation}
It is of interest to note that \eqref{Hartree-Fock} is the Hamiltonian equation of motion of a classical Hamiltonian system with phase space $\Gamma \deq l^2(\N) \otimes \HH^1(\R^3)$.

Define the map $\A$ from closed operators $A^{(p)}$ on  $\tilde{\mathcal{H}}_+^{(p)}$ to ``polynomial'' functions on phase space, through
\begin{multline*}
\A(A^{(p)})(\Psi) \;\deq\; \;\scalarb{\Psi^{\otimes p}}{A^{(p)} \Psi^{\otimes p}}\,
\\
=\; \int \dd \alpha_1 \cdots \dd \alpha_p \, \dd \beta_1 \cdots \dd \beta_p \; \bar{\Psi}(\alpha_p) \cdots \bar{\Psi}(\alpha_1) A^{(p)}(\alpha_1, \dots, \alpha_p; \beta_1, \dots, \beta_p) \, \Psi(\beta_1) \cdots \Psi(\beta_p)\,,
\end{multline*}
where $A^{(p)}(\alpha_1, \dots, \alpha_p; \beta_1, \dots, \beta_p)$ is the distribution kernel of $A^{(p)}$ (see \cite{bosons} for details). We denote by $\mathfrak{A}$ the linear hull of functions of the form $\A(A^{(p)})$, with $A^{(p)} \in \mathcal{B}(\tilde{\mathcal{H}}_+^{(p)})$.

The Hamilton function is given by
\begin{equation} \label{orbital Hamiltonian}
H \;\deq\; \A(\tilde{h}) + \frac{1}{2} \A(\tilde{\mathcal{W}})\,,
\end{equation}
where
\begin{equation*}
\mathcal{W} \;\deq\; W(\umat - E)
\end{equation*}
with $(E \Phi)(x_1,x_2) \deq \Phi(x_2,x_1)$ and $W$ is the two-particle operator defined by multiplication by $w(x_1 - x_2)$.
Written out in terms of components, \eqref{orbital Hamiltonian} reads
\begin{equation*}
H(\Psi) \;=\; \sum_{i \in \N} \scalar{\psi_i}{h \psi_i} + \frac{1}{2} \sum_{i,j \in \N}  \pb{\scalar{\psi_i \otimes \psi_j}{W \, \psi_i \otimes \psi_j} - \scalar{\psi_i \otimes \psi_j}{W \, \psi_j \otimes \psi_i}}\,.
\end{equation*}
Using Sobolev-type inequalities, one readily sees that $H$ is well-defined on $\Gamma$.

A short calculation shows that the Hartree-Fock equation is equivalent to
\begin{equation*}
\ii \partial_t \Psi \;=\; \partial_{\bar{\Psi}} H(\Psi)\,.
\end{equation*}
The symplectic form on $\Gamma$ is given by
\begin{equation*}
\omega \;=\; \ii \int \dd \alpha \; \dd \bar{\Psi}(\alpha) \wedge \dd \Psi(\alpha)\,,
\end{equation*}
which induces the Poisson bracket
\begin{equation} \label{Poisson bracket for Hartree-Fock}
\{\Psi(\alpha), \bar{\Psi}(\beta)\} \;=\; \ii \delta(\alpha - \beta) \,, \qquad \{\Psi(\alpha), \Psi(\beta)\} \;=\; \{\bar{\Psi}(\alpha), \bar{\Psi}(\beta)\} \;=\; 0 \,.
\end{equation}
Thus, for two observables $A, B \in \mathfrak{A}$,
\begin{equation*}
\h{A, B}(\Psi) \;=\; \ii \int \dd \alpha \; \p{\frac{\delta A}{\delta \Psi(\alpha)}(\Psi) \frac{\delta B}{\delta \bar{\Psi}(\alpha)}(\Psi) - \frac{\delta B}{\delta \Psi(\alpha)}(\Psi) \frac{\delta A}{\delta \bar{\Psi}(\alpha)}(\Psi)}\,.
\end{equation*}
The Hamiltonian equation of motion on $\Gamma$ is the Hartree-Fock equation \eqref{Hartree-Fock}.

The conservation laws of the Hartree-Fock flow can be understood in terms of symmetries of the Hamiltonian \eqref{orbital Hamiltonian}. One immediately sees that \eqref{orbital Hamiltonian} is invariant under the rotation $\Psi \mapsto (U \otimes \umat_{\LL^2(\R^3)}) \Psi$, where $U \in \mathcal{B}(l^2(\N))$ is unitary. A one-parameter group of such unitary transformations is generated by linear combinations of the functions $\re \scalar{\psi_i}{\psi_j}$ and $\im \scalar{\psi_i}{\psi_j}$, which Poisson-commute with the Hamiltonian $\eqref{orbital Hamiltonian}$. By Noether's principle, it follows that $\scalar{\psi_i}{\psi_j}$ is (at least formally) conserved. The energy $H$ is of course formally conserved as well.

In order to solve the Hartree-Fock equation \eqref{Hartree-Fock} with initial state $\Psi$, we rewrite it as an integral equation
\begin{equation} \label{integral Hartree-Fock}
\psi_i(t) \;=\; \ee^{-\ii t h} \psi_i - \ii \int_0^t \dd s \; \sum_{j \in \N} \pb{(w * \abs{\psi_j(s)}^2) \psi_i(s) - \p{w * (\psi_i(s) \bar{\psi}_j(s))} \psi_j(s)}\,.
\end{equation}
The Cauchy-problem for \eqref{integral Hartree-Fock} was solved in \cite{Zagatti}. We quote the relevant results:
\begin{lemma} \label{lemma: Zagatti}
Let $\Psi \in \tilde{\mathcal{H}}$. Then \eqref{integral Hartree-Fock} has a unique global solution $\Psi(\cdot) \in C(\R; \tilde{\mathcal{H}})$. Furthermore, the quantities $\scalar{\psi_i}{\psi_j}$ are conserved. In particular, $\norm{\Psi(t)} = \norm{\Psi}$.
\end{lemma}

\subsection{A Schwinger-Dyson expansion for the Hartree-Fock equation}
Our main tool is the Schwinger-Dyson expansion for the flow of the Hartree-Fock equation. We use the notation $(\cdot)_t$ to denote free time evolution generated by the free Hamiltonian $\A(\tilde{h})$. Explicitly,
\begin{equation*}
A_t(\psi_1, \psi_2, \dots) \;=\; A \pb{\ee^{-\ii th } \psi_1, \ee^{-\ii th } \psi_2, \dots}\,.
\end{equation*}

\begin{lemma} \label{classical Schwinger-Dyson expansion}
Let $A \in \mathfrak{A}$, $\nu > 0$, and $\Psi(t)$ be the solution of \eqref{integral Hartree-Fock} with initial data $\Psi$. Then, for small times $t$,
\begin{align*}
A(\Psi(t)) &\;=\; A_t(\Psi) + \int_0^t \dd s \; \frac{1}{2} \hb{\A(\tilde{\mathcal{W}}), A_{t-s}}(\Psi(s)) 
\\
&\;=\; \sum_{k = 0}^\infty \frac{1}{2^k} \int_{\Delta^k(t)} \dd \ul{t} \; \hB{\A(\tilde{\mathcal{W}}_{t_k}), \dots \hb{\A(\tilde{\mathcal{W}}_{t_1}), A_t}}(\Psi)\,,
\end{align*}
uniformly for $\Psi \in B_\nu \deq \{\Psi \in \tilde{\mathcal{H}} \,:\, \norm{\Psi}^2 \leq \nu\}$\,.
\end{lemma}

\begin{proof}
The proof of Lemma 7.1 in \cite{bosons} applies with virtually no modifications. One uses \eqref{norm of tilde operator}, the identity
\begin{equation*}
\A(\tilde{\mathcal{W}})_t \;=\; \A(\tilde{\mathcal{W}}_t) \;=\; \A\pb{(W_t (\umat - E))\,\tilde{\;}\,}\,,
\end{equation*}
and $\norm{E} = 1$.
\end{proof}

\section{The density matrix Hartree-Fock equation} \label{section: density matrix Hartree-Fock}
From now on, we only work with orthogonal sequence of orbitals
\begin{equation*}
\mathcal{K} \;\deq\; \{\Psi \in \tilde{\mathcal{H}} \,:\, \scalar{\psi_i}{\psi_j} = 0 ,\, i\neq j\}\,.
\end{equation*}
By Lemma \ref{lemma: Zagatti}, $\Psi \in \mathcal{K}$ implies that $\Psi(t) \in \mathcal{K}$ for all $t$. To each sequence of orbitals $\Psi$ we assign a one-particle density matrix
\begin{equation*}
\gamma_\Psi \;\deq\; \sum_{i \in \N} \ket{\psi_i} \bra{\psi_i}\,.
\end{equation*}
It is easy to see that this defines a mapping from $\mathcal{K}$ onto the set of density matrices
\begin{equation*}
\mathcal{D} \;\deq\; \h{\gamma \in \mathcal{L}^1(\mathcal{H}) \,:\, \gamma \geq 0}\,. 
\end{equation*}
Furthermore,
\begin{equation*}
\norm{\gamma_\Psi}_1 \;=\; \norm{\Psi}^2\,.
\end{equation*}
Conversely, one may recover $\Psi$ from $\gamma_\Psi$, up to ordering of the orbitals, by spectral decomposition. 
Furthermore, \eqref{extended operator calculated} implies that
\begin{equation} \label{observables computed using trace}
\A(\tilde{a}^{(p)})(\Psi) \;=\; \tr (a^{(p)} \gamma_\Psi^{\otimes p})\,.
\end{equation}

Let $\Psi(t)$ be a solution of \eqref{Hartree-Fock} with initial data $\Psi$ and write
\begin{equation*}
\gamma(t) \;=\; \gamma_{\Psi(t)}\,.
\end{equation*}
Then a short calculation shows that
\begin{equation} \label{Hartree-Fock for density matrices}
\ii \partial_t \gamma \;=\; [h, \gamma] + \tr_2 \q{\mathcal{W}, \gamma \otimes \gamma}\,,
\end{equation}
which is the Hartree-Fock equation for density matrices.
As an integral equation in the interaction picture, this reads
\begin{equation} \label{integral Hartree-Fock for density matrices}
\gamma(t) \;=\; \ee^{-\ii th} \, \gamma \, \ee^{\ii th} - \ii \int_0^t \dd s \; \ee^{- \ii (t-s) h} \, \tr_2 \q{\mathcal{W}, \gamma(s) \otimes \gamma(s)} \, \ee^{\ii (t-s) h} \,.
\end{equation}
Sometimes it is convenient to rewrite this using the shorthand
\begin{equation} \label{interaction picture density matrix}
\tilde{\gamma}(t) \;\deq\; \ee^{\ii th} \,\gamma(t) \,\ee^{-\ii th}\,.
\end{equation}
Then \eqref{integral Hartree-Fock for density matrices} is equivalent to
\begin{equation} \label{integral Hartree-Fock for density matrices interaction}
\tilde{\gamma}(t) \;=\; \gamma - \ii \int_0^t \dd s \; \tr_2 \q{\mathcal{W}_s, \tilde{\gamma}(s) \otimes \tilde{\gamma}(s)}\,.
\end{equation}

\begin{lemma} \label{Hartree-Fock solves density matrix Hartree-Fock}
Let $\Psi(t)$ be the solution of \eqref{integral Hartree-Fock}. Then $\gamma_{\Psi(t)}$ solves \eqref{integral Hartree-Fock for density matrices}.
\end{lemma}
\begin{proof}
Let $a^{(1)} \equiv a \in \mathcal{B}(\mathcal{H})$. From 
Lemma \ref{classical Schwinger-Dyson expansion} we get
\begin{equation} \label{integral equation for one-particle observable}
\A(\tilde{a})(\Psi(t)) \;=\; \A(\tilde{a}_t)(\Psi) + \int_0^t \dd s \; \hb{\A(\tilde{\mathcal{W}}), \A(\tilde{a}_{t-s})}(\Psi(s))\,. 
\end{equation}
Now \eqref{Poisson bracket for Hartree-Fock} and \eqref{observables computed using trace} imply
\begin{align*}
\hb{\A(\tilde{\mathcal{W}}), \A(\tilde{a})}(\Psi) &\;=\; \ii \A\pb{\qb{\tilde{\mathcal{W}} , \tilde{a} \otimes \umat}}(\Psi) 
\\
&\;=\; \ii \tr \pb{\qb{\mathcal{W} , a \otimes \umat} \gamma_\Psi \otimes \gamma_\Psi}
\\
&\;=\; - \ii \tr \pb{(a \otimes \umat) \qb{\mathcal{W} , \gamma_\Psi \otimes \gamma_\Psi}}\,.
\end{align*}
Thus \eqref{integral equation for one-particle observable} reads
\begin{align*}
\tr (a \, \gamma_{\Psi(t)}) &\;=\; \tr (a_t \, \gamma_\Psi) - \ii \int_0^t \dd s \; \tr \pb{(a_{t - s} \otimes \umat) \qb{\mathcal{W} , \gamma_{\Psi(s)} \otimes \gamma_{\Psi(s)}}}
\\
&\;=\; \tr \pb{a \, \ee^{\ii t h} \gamma_\Psi \ee^{- \ii th}} - \ii \int_0^t \dd s \; \tr \pB{a \, \ee^{- \ii (t-s) h} \, \tr_2 \q{\mathcal{W}, \gamma_{\Psi(s)} \otimes \gamma_{\Psi(s)}} \, \ee^{\ii (t-s) h} }\,.
\end{align*}
Since $a \in \mathcal{B}(\mathcal{H})$ was arbitrary, this is equivalent to \eqref{integral Hartree-Fock for density matrices}.
\end{proof}

\section{Slater determinants} \label{section: Slater determinants}
The Hartree-Fock equation naturally describes the time evolution of quasi-free states \cite{BratteliRobinson}. Let $\omega_\gamma$ be the quasi-free state corresponding to the one-particle state $\gamma \in \mathcal{D}$. Define 
\begin{equation*}
\gamma^{(p)}(x_1, \dots, x_p; y_1, \dots, y_p) \;\deq\; \omega_\gamma\p{a^*(y_p) \cdots a^*(y_1) a(x_1) \cdots a(x_p)}\,,
\end{equation*}
where $a^*(x), a(x)$ are the usual creation and annihilation operators of the CAR algebra over $\mathcal{H}$ (see e.g.\ \cite{BratteliRobinson}).
The quasifreeness of $\omega_\gamma$ means that
\begin{equation*}
\gamma^{(p)}(x_1, \dots, x_p; y_1, \dots, y_p) \;=\; \det \p{\gamma(x_i; y_j)}_{i,j}\,.
\end{equation*}
In other words, $\gamma^{(p)}$ is the operator kernel of
\begin{equation} \label{p-particle action of quasi-free state}
\gamma^{(p)} \;=\; \gamma^{\otimes p} \, \Sigma_-^{(p)} \,,
\end{equation}
where
\begin{equation*}
\Sigma_-^{(p)} \;\deq\; p! P_-^{(p)}\,.
\end{equation*}
For the following calculations it is convenient to introduce the symbol $\epsilon_{i_1 \dots i_p}^{j_1 \dots j_p}$, which is equal to $\sgn \sigma$ if $i_1, \dots, i_p$ are disjoint and there is a permutation $\sigma \in S_p$ such that $(i_1, \dots, i_p) = (j_{\sigma(1)}, \dots, j_{\sigma(p)})$, and equal to $0$ otherwise. Also, for the remainder of this section, summation over any index appearing twice in an equation is implied.
\begin{lemma} \label{lemma: bound on gamma_p}
Let $\gamma \in \mathcal{D}$ with $\tr \gamma = 1$. Then $\tr \gamma^{(p)} \leq 1$.
\end{lemma}
\begin{proof}
There is an orthonormal basis $(\varphi_i)_{i \in \N}$ and a sequence of nonnegative numbers $(\lambda_i)_{i \in \N}$ such that $\sum_i \lambda_i = 1$ and $\gamma = \sum_i \lambda_i \ket{\varphi_i} \bra{\varphi_i}$. Therefore,
\begin{equation*}
\gamma^{(p)} \;=\; \epsilon_{i_1 \dots i_p}^{j_1 \dots j_p} \lambda_{i_1} \cdots \lambda_{i_p} \ket{\varphi_{i_1} \otimes \cdots \otimes \varphi_{i_p}} \bra{\varphi_{j_1} \otimes \cdots \otimes \varphi_{j_p}}
\end{equation*}
This yields
\begin{align*}
\tr \gamma^{(p)} &\;=\; \epsilon_{i_1 \dots i_p}^{j_1 \dots j_p} \lambda_{i_1} \cdots \lambda_{i_p} \delta_{i_1 k_1} \cdots \delta_{i_p k_p} \delta_{j_1 k_1} \cdots \delta_{j_p k_p}
\\
&\;=\; \sum_{i_1, \dots, i_p \text{ disjoint}} \lambda_{i_1} \cdots \lambda_{i_p}
\\
&\;\leq\; \sum_{i_1, \dots, i_p} \lambda_{i_1} \cdots \lambda_{i_p} \;=\; 1\,.
\end{align*}
\end{proof}

Next, we introduce a special class of quasi-free states, described by \emph{Slater determinants}.
Let $N \in \N$  and take an orthonormal sequence of orbitals $\Phi_N = (\varphi_1, \dots, \varphi_N)$.
We define the Slater determinant
\begin{equation*}
S(\Phi_N) \;\deq\; \varphi_1 \wedge \cdots \wedge \varphi_N \;\deq\; \sqrt{N!} \, P^{(N)}_- \, \varphi_1 \otimes \cdots \otimes \varphi_N \;\in\; \mathcal{H}^{(N)}_-\,.
\end{equation*}
Note that the normalization is chosen so that $\norm{S(\Phi_N)} = 1$. The corresponding $N$-particle density matrix is
\begin{equation*}
\Gamma_N \;\deq\; \ket{S(\Phi_N)} \bra{S(\Phi_N)}\,.
\end{equation*}
One finds for the $p$-particle marginals
\begin{align}
\Gamma_N^{(p)} &\;\deq\; \tr_{p+1 \dots N} \Gamma_N
\notag \\
&\;=\; \tr_{p+1 \dots N} \frac{1}{N!} \, \epsilon_{i_1 \dots i_N} \epsilon_{j_1 \dots j_N} \, \ketb{\varphi_{i_1} \otimes \cdots \otimes \varphi_{i_N}} \brab{\varphi_{j_1} \otimes \cdots \otimes \varphi_{j_N}}
\notag \\ \label{p-particle marginals computed}
&\;=\; \frac{(N-p)!}{N!} \,\epsilon_{i_1 \dots i_p}^{j_1 \dots j_p} \, \ketb{\varphi_{i_1} \otimes \cdots \otimes \varphi_{i_p}} \brab{\varphi_{j_1} \otimes \cdots \otimes \varphi_{j_p}}\,.
\end{align}

In order to relate the sequence $\Phi_N$ to the results of the previous sections, we define the normalised sequence
\begin{equation} \label{definition of Psi_N}
\Psi_N \;\deq\; \frac{1}{\sqrt{N}} (\varphi_1, \dots, \varphi_N, 0, \dots)\,.
\end{equation}
Thus, $\Psi_N \in \mathcal{K}$ and $\norm{\Psi_N} = 1$. It is trivial to check that $\Psi_N(t)$ is a solution of \eqref{Hartree-Fock} if and only if $\Phi_N(t)$ is a solution of \eqref{rescaled Hartree-Fock}. Similarly,
$\Psi_N(t)$ is a solution of \eqref{integral Hartree-Fock} if and only if $\Phi_N(t) = (\varphi_1(t), \dots, \varphi_N(t))$ is a solution of 
\begin{equation} \label{rescaled integral Hartree-Fock}
\varphi_i(t) \;=\; \ee^{-\ii t h} \varphi_i - \frac{\ii}{N} \int_0^t \dd s \; \sum_{j =1}^N \pb{(w * \abs{\varphi_j(s)}^2) \varphi_i(s) - \p{w * (\varphi_i(s) \bar{\varphi}_j(s))} \varphi_j(s)}\,.
\end{equation}

Next, from \eqref{p-particle marginals computed} we find
\begin{equation*}
\gamma_N \;\deq\; \Gamma_N^{(1)} \;=\; \frac{1}{N} \sum_{i = 1}^N \ket{\varphi_i} \bra{\varphi_i} \;=\; \sum_{i \in \N} \ket{\psi_i} \bra{\psi_i} \;=\; \gamma_{\Psi_N}\,.
\end{equation*}
Thus \eqref{p-particle action of quasi-free state} implies that
\begin{equation} \label{relation between marginals}
\gamma_{N}^{(p)} \;=\; \frac{1}{N^p} \, \epsilon_{i_1 \dots i_p}^{j_1 \dots j_p} \, \ketb{\varphi_{i_1} \otimes \cdots \otimes \varphi_{i_p}} \brab{\varphi_{j_1} \otimes \cdots \otimes \varphi_{j_p}} \;=\; \frac{p!}{N^p} \binom{N}{p} \, \Gamma_N^{(p)}\,.
\end{equation}
Thus, Slater determinants determine quasi-free states by their reduced $p$-particle marginals. The normalisation $\frac{p!}{N^p} \binom{N}{p}$ differs slightly from the usual normalisation $1$ of quasi-free states, but in the limit $N \to \infty$ this difference vanishes.
Note also that
\begin{equation} \label{operator norm of gamma}
\norm{\gamma_N} \;=\; \frac{1}{N}\,.
\end{equation}
This is a special case of the well-known statement (see e.g.\ \cite{LiebLoss}) that $\normb{\tr_{2 \dots N} \Gamma} \leq N^{-1}$, for any $N$-particle density matrix $\Gamma$ \,\footnote{This can also be inferred from \eqref{operator norm of gamma} by writing $\Gamma$ as a linear combination of projectors.}. 


\section{Proof of convergence and the mean-field limit} \label{section: convergence}
We now turn to the proof of our main result. We use the graph expansion scheme for the Schwinger-Dyson expansion developed in \cite{bosons}. 

\subsection{The Schwinger-Dyson graph expansion}
For the convenience of the reader we summarise the relevant results of the graph expansion in \cite{bosons}. For details and proofs we refer to \cite{bosons}.
Let $a^{(p)} \in \mathcal{B}(\mathcal{H}^{(p)}_-)$ and define its extension $\widehat{\A}_N(a^{(p)})$ to $N$-particle space $\mathcal{H}^{(N)}_-$ through
\begin{equation} \label{second quantisation of a}
\widehat{\A}_N(a^{(p)}) \;\deq\;
\begin{cases}
\frac{p!}{N^p} \binom{N}{p} P_- (a^{(p)} \otimes \umat^{(n-p)}) P_-\,, &N\geq p\,,
\\
0\,, &N<p\,.
\end{cases}
\end{equation}
An immediate consequence of \eqref{relation between marginals} is
\begin{equation} \label{expectation of second-quantised operator}
\scalarb{S(\Phi_N)}{\widehat{\A}_N(a^{(p)}) \, S(\Phi_N)} \;=\; \tr \p{a^{(p)} \gamma_N^{(p)}}\,.
\end{equation}
The Schwinger-Dyson series for the time evolution of $\widehat{\A}_N(a^{(p)})$ is given by
\begin{equation} \label{loop expansion}
\me^{\mi t H_N} \, \widehat{\A}_N (a^{(p)}) \, \me^{- \mi t H_N} \;=\; \sum_{k = 0}^\infty \sum_{l = 0}^k \frac{1}{N^l} \, \widehat{\A}_N \pb{G^{(k,l)}_t(a^{(p)})}\,.
\end{equation}
The Hamiltonian $H_N$ with mean-field scaling is defined in \eqref{mean-field Hamiltonian}. The series \eqref{loop expansion}
converges in norm for all times $t$; for small times $t$ the convergence is uniform in $N$. The $(p+k-l)$-particle operator $G^{(k,l)}_t(a^{(p)})$ is a shorthand notation for
\begin{equation}
G^{(k,l)}_{t}(a^{(p)}) \;\deq\; \int_{\Delta^k(t)} \md \ul{t} \; G^{(k,l)}_{t,\ul{t}}(a^{(p)})\,,
\end{equation}
where $\ul{t} = (t_1, \dots, t_k)$ and $\Delta^k(t)$ is the $k$-simplex $\{(t_1, \dots, t_k) : 0 < t_k < \dots < t_1 < t\}$. The operators $G^{(k,l)}_{t}(a^{(p)})$ are recursively defined by
\begin{align}
G^{(k,l)}_{t,t_1,\dots,t_k}(a^{(p)}) 
\;=\; \;&\ii P_- \sum_{i = 1}^{p+k-l-1} \qB{W_{i \, p+k-l, t_k}, G^{(k-1,l)}_{t, t_1, \dots, t_{k - 1}}(a^{(p)}) \otimes \umat} P_- 
\notag \\ \label{recursive definition of graphs}
&\ii P_-  \sum_{1 \leq i < j \leq p+k-l} \qB{W_{ij, t_k}, G^{(k-1, l-1)}_{t, t_1, \dots, t_{k - 1}}(a^{(p)})} P_- \,,
\end{align}
as well as $G^{(0,0)}_t(a^{(p)}) \deq a^{(p)}_t$. If $l < 0$ or $l > k$ then $G^{(k,l)}_{t,t_1,\dots,t_k}(a^{(p)}) = 0$. Here, as before, a time subscript denotes free time evolution:
\begin{equation*}
a^{(p)}_t \;\deq\; \ee^{\ii \sum_{i} h_i t} \, a^{(p)} \, \ee^{-\ii \sum_i h_i t}\,.
\end{equation*}

The operator $G^{(k,l)}_{t}(a^{(p)})$ may be written as a sum over graphs\footnote{called graph structures in \cite{bosons}}:
\begin{equation} \label{graph expansion}
G^{(k, l)}_t(a^{(p)}) \;=\; \frac{\mi^k}{2^k} \sum_{\mathcal{Q} \in \mathscr{Q}(p,k,l)} i_{\mathcal{Q}} \int_{\Delta^k_{\mathcal{Q}}(t)} \md \ul{t} \; G^{(k, l)(\mathcal{Q})}_{t,t_1,\dots, t_k}(a^{(p)})\,,
\end{equation}
where $i_{\mathcal{Q}} \in \{0,1\}$, and $\Delta^k_{\mathcal{Q}}(t) \subset [0,t]^k$, and
$\mathscr{Q}(p,k,l)$ is a set of graphs satisfying
\begin{equation} \label{bound on number of graphs}
\abs{\mathscr{Q}(p,k,l)} \;\leq\; 2^k \binom{k}{l} \, \binom{2p + 3k}{k} \, (p+k-l)^l\,.
\end{equation}
The operator $G^{(k, l)(\mathcal{Q})}_{t,t_1,\dots, t_k}(a^{(p)})$ is an \emph{elementary term}, indexed by the graph $\mathcal{Q}$, of the form
\begin{equation} \label{elementary term}
P_- \, W_{i_1 j_1, t_{v_1}} \cdots W_{i_r j_r, t_{v_r}} \,\pb{a_t^{(p)} \otimes \umat^{(k-l)}} \, W_{i_{r+1} j_{r+1}, t_{v_{r+1}}} \cdots W_{i_k j_k, t_{v_k}} P_- \,,
\end{equation}
where $r = 0, \dots, k$.

The operator norm of $G^{(k,l)}_{t}(a^{(p)})$ may now be estimated using the dispersive estimate
\begin{equation} \label{general Kato smoothing}
\int \dd t \; \normb{\abs{x}^{-1} \ee^{\ii t \Delta} \psi}^2 \;\leq\; \pi \norm{\psi}^2\,.
\end{equation}
Going to centre of mass coordinates and using Cauchy-Schwarz, one sees that \eqref{general Kato smoothing} implies
\begin{equation} \label{Kato smoothing}
\int_0^t \dd s \; \normb{W_{i j, s} \, \varphi} \;\leq\; \sqrt{\frac{\pi \kappa^2 t}{2}} \, \norm{\varphi}\,.
\end{equation}
Together with the estimate \eqref{bound on number of graphs}, it is now easy to argue, as in \cite{bosons}, that \eqref{loop expansion} converges uniformly in $N$ for small times $t$. Moreover, the large-$N$ asymptotics of the Schwinger-Dyson series \eqref{loop expansion} is given by the tree terms: for small times $t$ we have
\begin{equation} \label{one loop expansion}
\me^{\mi t H_N} \, \widehat{\A}_N (a^{(p)}) \, \me^{- \mi t H_N} \;=\; \sum_{k = 0}^\infty \widehat{\A}_N \pb{G^{(k,0)}_t(a^{(p)})} + L_N(t)\,,
\end{equation}
where $L_N(t)$, corresponding to the sum of all ``loop terms'', satisfies the estimate 
\begin{equation} \label{estimate on loop terms}
\norm{L_N(t)} \;\leq\; C(p,\kappa, t) \norm{a^{(p)}} N^{-1}\,,
\end{equation}
for small times $t$.

\subsection{Convergence of the Hartree-Fock time evolution to the tree terms}
We now give the main argument of our proof. We show that the Hartree-Fock time evolution is asymptotically ($N \to \infty$) given by the tree terms (i.e. the terms $l = 0$) of the Schwinger-Dyson series \eqref{loop expansion}. This result is summarised in Lemma \ref{lemma HF and tree terms} below.

Let $\Psi = (\psi_i)_{i = 1}^\infty \in \mathcal{K}$, and denote by $\Psi(t)$ the solution of the Hartree-Fock equation \eqref{integral Hartree-Fock} with initial data $\Psi$. Let $\gamma(t) = \gamma_{\Psi(t)}$ be the associated one-particle density matrix. 

By choosing $A = \A(\tilde{a}^{(p)})$, $a^{(p)} \in \mathcal{B}(\mathcal{H}^{(p)}_-)$, in Lemma \ref{classical Schwinger-Dyson expansion} and mimicking the proof of Lemma \ref{Hartree-Fock solves density matrix Hartree-Fock} one finds that
\begin{equation} \label{Hartree-Fock on tensor power, step 1}
\tr \pb{a^{(p)} \, \gamma(t)^{\otimes p}} \;=\; \tr \pb{a^{(p)}_t \gamma^{\otimes p}} - \ii \int_0^t \dd s \sum_{i = 1}^p \tr \pB{a^{(p)}_{t-s} \tr_{p+1} \qb{\mathcal{W}_{i \, p+1}, \gamma(s)^{\otimes (p+1)}}}
\end{equation}
It is convenient to use the representation $\tilde{\gamma}(t)$ defined in \eqref{interaction picture density matrix}. Using the substitution $a^{(p)} \mapsto a^{(p)}_{-t}$ in \eqref{Hartree-Fock on tensor power, step 1} we get
\begin{equation*}
\tr \pb{a^{(p)} \, \tilde{\gamma}(t)^{\otimes p}} \;=\; \tr \pb{a^{(p)} \gamma^{\otimes p}} - \ii \int_0^t \dd s \sum_{i = 1}^p \tr \pB{a^{(p)} \tr_{p+1} \qb{\mathcal{W}_{i \, p+1, s}, \tilde{\gamma}(s)^{\otimes (p+1)}}}
\end{equation*}
Recall that $\mathcal{W}_{ij} = W_{ij} (\umat - E_{ij})$. Also, $E_{ij}$ commutes with $W_{ij}$ and with $\tilde{\gamma}(s)^{\otimes (p+1)}$. Thus, $\Sigma^{(p)}_- a^{(p)} = p! \, a^{(p)}$ implies
\begin{equation} \label{Hartree-Fock on tensor power}
\tilde{\gamma}(t)^{\otimes p} \, \Sigma_-^{(p)} \;=\; \gamma^{\otimes p} \, \Sigma_-^{(p)} - \mi \int_0^t \dd s \; \tr_{p+1} \q{\sum_{i = 1}^p W_{i p+1, s}, \tilde{\gamma}(s)^{\otimes (p+1)} (\umat - E_{i p+1})} \Sigma_-^{(p)}.
\end{equation}

On the other hand, using
\begin{equation*}
\Sigma^{(p+1)}_- \;=\; \pbb{\umat - \sum_{j = 1}^p E_{j p+1}}\Sigma_-^{(p)}
\end{equation*}
we find
\begin{align*}
\tr_{p+1} \q{\sum_{i = 1}^p W_{i p+1, s}, \tilde{\gamma}(s)^{\otimes (p+1)} \, \Sigma_-^{(p+1)}}
&\;=\; \tr_{p+1} \q{\sum_{i = 1}^p W_{i p+1, s}, \tilde{\gamma}(s)^{\otimes (p+1)} \pbb{\umat - \sum_{j = 1}^p E_{j p+1}}\Sigma_-^{(p)}}
\\
&\;=\; \tr_{p+1} \q{\sum_{i = 1}^p W_{i p+1, s}, \tilde{\gamma}(s)^{\otimes (p+1)} \pbb{\umat - \sum_{j = 1}^p E_{j p+1}}} \Sigma_-^{(p)}\,.
\end{align*}
Together with \eqref{Hartree-Fock on tensor power} this yields
\begin{equation} \label{HF evolution of a product state}
\tilde{\gamma}(t)^{\otimes p} \, \Sigma_-^{(p)} \;=\; \gamma^{\otimes p} \, \Sigma_-^{(p)} - \mi \int_0^t \dd s \; \tr_{p+1} \q{\sum_{i = 1}^p W_{i p+1, s}, \tilde{\gamma}(s)^{\otimes (p+1)} \, \Sigma_-^{(p+1)}}
+ R_p(t)\,,
\end{equation}
with an error term
\begin{equation} \label{definition of R}
R_p(t) \;\deq\;  - \mi \sum_{1 \leq i \neq j \leq p} \int_0^t \dd s \; \tr_{p+1} \q{W_{i p+1, s}, \tilde{\gamma}(s)^{\otimes (p+1)} E_{j p+1}} \Sigma_-^{(p)}.
\end{equation}
The partial trace is most conveniently computed using operator kernels. We find
\begin{multline*}
\p{W_{i p+1, s} \tilde{\gamma}(s)^{\otimes (p+1)} E_{j p+1}} (x_1, \dots, x_{p+1}; y_1, \dots, y_{p+1})
\\
=\; \int \dd z_1 \, \dd z_2 \; \q{\prod_{r \neq i,j} \tilde{\gamma}(s)(x_r; y_r)}
\, W_s(x_i, x_{p+1}; z_1, z_2) \, \tilde{\gamma}(s)(z_1; y_i) \, \tilde{\gamma}(s)(x_j; y_{p+1}) \, \tilde{\gamma}(s)(z_2; y_j)\,,
\end{multline*}
so that
\begin{align*}
&\qquad \tr_{p+1} \p{W_{i p+1, s} \tilde{\gamma}(s)^{\otimes (p+1)} E_{j p+1}} (x_1, \dots, x_p; y_1, \dots, y_p)
\\
&=\; \int \dd z_1 \, \dd z_2 \, \dd z_3 \; \q{\prod_{r \neq i,j} \tilde{\gamma}(s)(x_r; y_r)}
\, W_s(x_i, z_3; z_1, z_2) \, \tilde{\gamma}(s)(z_1; y_i) \, \tilde{\gamma}(s)(x_j; z_3) \, \tilde{\gamma}(s)(z_2; y_j)
\\
&=\; \p{\tilde{\gamma}_j(s) W_{i j,s} \tilde{\gamma}(s)^{\otimes p}} (x_1, \dots, x_p; y_1, \dots, y_p)\,.
\end{align*}
The second term of the commutator in \eqref{definition of R} is the adjoint of the first and we get 
\begin{align}
R_p(t) &\;=\;  - \mi \sum_{1 \leq i \neq j \leq p} \int_0^t \dd s \; \pB{\tilde{\gamma}_j(s) W_{i j,s} \tilde{\gamma}(s)^{\otimes p} - \tilde{\gamma}(s)^{\otimes p} W_{i j,s} \tilde{\gamma}_j(s)} \Sigma^{(p)}_-
\notag \\
&\;=\;  - \mi \sum_{1 \leq i \neq j \leq p} \int_0^t \dd s \; \pB{\tilde{\gamma}_j(s) W_{i j,s} \tilde{\gamma}(s)^{\otimes p} \Sigma^{(p)}_- - \tilde{\gamma}(s)^{\otimes p} \Sigma^{(p)}_- W_{i j,s} \tilde{\gamma}_j(s)}\,.
\end{align}

We proceed to show that, up to an error term, the expansion of the Hartree-Fock time-evolution is equal to the tree terms of the microscopic quantum-mechanical evolution. Let $a^{(p)} \in \mathcal{B}(\mathcal{H}^{(p)}_-)$. Using \eqref{HF evolution of a product state} we find
\begin{align*}
\tr \p{a^{(p)} \, \gamma(t)^{\otimes p} \, \Sigma_-^{(p)}} &\;=\; \tr \p{a_t^{(p)} \, \tilde{\gamma}(t)^{\otimes p} \, \Sigma^{(p)}_-}
\\
&\;=\; \tr \p{a_t^{(p)} \, \gamma^{\otimes p} \, \Sigma^{(p)}_-} + \ii \int_0^t \dd s \, \sum_{i = 1}^p \tr \p{\q{W_{i \, p+1,s}, a^{(p)}_t \otimes \umat} \, \tilde{\gamma}(s)^{\otimes (p+1)} \, \Sigma^{(p+1)}_-} 
\\
&\qquad{}+{} \tr \p{a^{(p)}_t \, R_p(t)}\,.
\end{align*}
Iterating this $K$ times yields our main expansion
\begin{align}
\tr \p{a_t^{(p)} \, \tilde{\gamma}(t)^{\otimes p} \, \Sigma_-^{(p)}} &\;=\; \sum_{k = 0}^{K-1} \int_{\Delta^k(t)} \dd \ul{t} \; \tr \p{G_{t, \ul{t}}^{(k,0)}(a^{(p)}) \, \gamma^{\otimes (p+k)} \Sigma_-^{(p+k)}}
\notag \\
&\qquad {}+{} \int_{\Delta^K(t)} \dd \ul{t} \; \tr \p{G_{t, \ul{t}}^{(K,0)}(a^{(p)}) \, \tilde{\gamma}(t_K)^{\otimes (p+K)} \Sigma_-^{(p+K)}}
\notag \\ \label{main expansion}
&\qquad {}+{} \sum_{k = 0}^{K-1} \sum_{1 \leq i \neq j \leq p+k} R^{k}_{ij}(t)\,,
\end{align}
where
\begin{multline*}
R^{k}_{ij}(t) \;\deq\; - \ii \, \int_{\Delta^{k+1}(t)} \dd \ul{t} \; \tr \Bigl(G^{(k,0)}_{t, t_1, \dots, t_k}(a^{(p)}) \, \tilde{\gamma}_j(t_{k+1}) \, W_{ij, t_{k + 1}} \, \tilde{\gamma}(t_{k+1})^{\otimes (p+k)} \, \Sigma_-^{(p+k)} \\
- W_{ij, t_{k+1}} \, \tilde{\gamma}_j(t_{k+1}) \, G^{(k,0)}_{t, t_1, \dots, t_k}(a^{(p)}) \, \tilde{\gamma}(t_{k+1})^{\otimes (p+k)} \, \Sigma_-^{(p+k)}
\Bigr)\,.
\end{multline*}

We now derive a bound on $R^{k}_{ij}(t)$. Let us concentrate on the first term, which we rewrite using the renaming $t_{k+1} \to s$ as
\begin{equation} \label{half of the error term}
\int_{\Delta^k(t)} \dd \ul{t} \int_0^{\wedge \ul{t}} \dd s \; \tr \p{G^{(k,0)}_{t, t_1, \dots, t_k}(a^{(p)}) \, \tilde{\gamma}_j(s) \, W_{ij, s} \, \tilde{\gamma}(s)^{\otimes (p+k)} \, \Sigma_-^{(p+k)}}\,,
\end{equation}
where $\wedge \ul{t} \deq \min\{t_1, \dots, t_k\}$. The idea is to use a tree expansion on $\tilde{\gamma}(s)$.
\begin{lemma}
Let $a^{(p)} \in \mathcal{B}(\mathcal{H}_-^{(p)})$. For small times we have the tree expansion
\begin{equation} \label{tree expansion of p-particle observable}
\tr \pb{a^{(p)} \tilde{\gamma}(t)^{\otimes p} \Sigma_-^{(p)}} \;=\;
\sum_{k = 0}^\infty \int_{\Delta^k(t)} \dd \ul{t} \; \tr\pb{T^{(k)}_{\ul{t}}(a^{(p)}) \gamma^{\otimes (p+k)} \Sigma_-^{(p)}}\,,
\end{equation}
where $T^{(k)}_{\ul{t}}$ is the linear operator defined by $T^{(0)}(a^{(p)}) \deq a^{(p)}$ and
\begin{equation*}
T^{(k)}_{t_1 \dots t_k}(a^{(p)}) \;=\; \ii \sum_{i = 1}^{p+k-1} \q{\mathcal{W}_{i \, p + k, t_k}, T^{(k-1)}_{t_1 \dots t_{k-1}}(a^{(p)}) \otimes \umat}\,.
\end{equation*}
\end{lemma}
\begin{proof}
Lemma \ref{classical Schwinger-Dyson expansion} applied to $A = \A(\tilde{a}^{(p)})$ yields
\begin{equation*}
\tr \pb{a^{(p)} \tilde{\gamma}(t)^{\otimes p}} \;=\;
\sum_{k = 0}^\infty \int_{\Delta^k(t)} \dd \ul{t} \; \tr\pb{T^{(k)}_{\ul{t}}(a^{(p)}) \gamma^{\otimes (p+k)}}\,.
\end{equation*}
The claim then follows by noting that $\Sigma^{(p)}_-a^{(p)} = p! a^{(p)}$ and that $\sum_{i = 1}^{p+k-1} \mathcal{W}_{i\, p+1, t_k}$ commutes with $\Sigma_-^{(p)}$.
The convergence of the series is shown below.
\end{proof}

Thus \eqref{half of the error term} is equal to
\begin{equation*}
\sum_{k' = 0}^\infty \int_{\Delta^k(t)} \dd \ul{t} \int_0^{\wedge \ul{t}} \dd s \int_{\Delta^{k'}(s)} \dd \ul{t}' \;
\tr\hB{T_{\ul{t}'}^{(k')} \pB{G_{t,\ul{t}}^{(k,0)}(a^{(p)}) \tilde{\gamma}_j(s) W_{ij, s}} \gamma^{\otimes (p+k+k')} \Sigma_-^{(p+k)}}\,.
\end{equation*}

Next, we recall from \eqref{graph expansion} that $G^{(k,0)}_{t, t_1, \dots, t_k}(a^{(p)})$ can be written as a sum over tree graphs $\mathcal{Q} \in \mathscr{Q}(p,k,0)$ of elementary terms of the form \eqref{elementary term}. Also, since the definition of $T^{(k)}_{t_1, \dots, t_k}(a^{(p)})$ is the same as the definition of $G^{(k,0)}_{0, t_1, \dots, t_k}(a^{(p)})$ with $W$ replaced by $\mathcal{W}$, we immediately get that $T^{(k)}_{t_1, \dots, t_k}(a^{(p)})$ is equal to a sum over tree graphs $\mathcal{Q} \in \mathscr{Q}(p,k,0)$ of elementary terms of the form
\begin{equation*}
P_- \, \mathcal{W}_{i_1 j_1, t_{v_1}} \cdots \mathcal{W}_{i_r j_r, t_{v_r}} \,\pb{a_t^{(p)} \otimes \umat^{(k-l)}} \, \mathcal{W}_{i_{r+1} j_{r+1}, t_{v_{r+1}}} \cdots \mathcal{W}_{i_k j_k, t_{v_k}} P_- \,.
\end{equation*}
This implies that the series \eqref{tree expansion of p-particle observable} converges for small times.

Applying the tree expansion to both $G^{(k,0)}_{t, t_1, \dots, t_k}(a^{(p)})$ and $\tilde{\gamma}(s)^{\otimes (p+k)}$ in \eqref{half of the error term}, we see that \eqref{half of the error term} is equal to 
\begin{multline*}
\sum_{k' = 0}^\infty \frac{\ii^{k+k'}}{2^{k+k'}}\sum_{\mathcal{Q} \in \mathscr{Q}(p,k,0)} \sum_{\mathcal{Q}' \in \mathscr{Q}(p+k, k',0)} i_{\mathcal{Q}} i_{\mathcal{Q}'} \int_{\Delta^k_{\mathcal{Q}}(t)} \dd \ul{t} \int_0^{\wedge \ul{t}} \dd s \int_{\Delta^{k'}_{\mathcal{Q}'}(s)} \dd \ul{t}'
\\
\tr\hB{A \, a^{(p)}_{1\dots p,t} \, B P_-^{(p+k)} 
\tilde{\gamma}_j(s) W_{ij, s} \, C \, \gamma^{\otimes (p+k+k')} \Sigma_-^{(p+k)}}\,,
\end{multline*}
where $A, B, C$ are operators that depend on $(\mathcal{Q}, \mathcal{Q}', k, k', \ul{t}, \ul{t})$. $A$, $B$ and $C$ are each a product of operators of the form $W_{i' j', r}$, or $\mathcal{W}_{i' j', r}$, where $r$ stands for a time variable in $\{t_1, \dots, t_k, t_1, \dots, t_{k'}\}$. Moreover, the product $ABC$ contains $k$ $W$'s and $k'$ $\mathcal{W}$'s. Finally, each time variable in $t_1, \dots, t_k, t_1, \dots, t_{k'}$ appears exactly once in the product $ABC$.

Let $\varphi \in \mathcal{H}^{\otimes (p+k+k')}$ and estimate
\begin{multline*}
I \;\deq\; \normbb{\sum_{\mathcal{Q} \in \mathscr{Q}(p,k,0)} \sum_{\mathcal{Q}' \in \mathscr{Q}(p+k, k',0)} i_{\mathcal{Q}} i_{\mathcal{Q}'} \int_{\Delta^k_{\mathcal{Q}}(t)} \dd \ul{t} \int_0^{\wedge \ul{t}} \dd s \int_{\Delta^{k'}_{\mathcal{Q}'}(s)} \dd \ul{t}' \;
A \, a^{(p)}_{1\dots p,t} \, B P_-^{(p+k)} 
\tilde{\gamma}_j(s) W_{ij, s} \, C \, \varphi}
\\
\leq\; 
\sum_{\mathcal{Q} \in \mathscr{Q}(p,k,0)} \sum_{\mathcal{Q}' \in \mathscr{Q}(p+k, k',0)} \int_{[0,t]^k} \dd \ul{t} \int_0^t \dd s \int_{[0,t]^{k'}} \dd \ul{t}' \;
\normB{A \, a^{(p)}_{1\dots p, t} \, B P_-^{(p+k)} 
\tilde{\gamma}_j(s) W_{ij, s} \, C \, \varphi}
\end{multline*}
We now perform all time integrations, starting from the left, and using at each step the estimate \eqref{Kato smoothing} as well as
\begin{equation*}
\int_0^t \dd r \; \normb{\mathcal{W}_{i' j', r} \varphi} \;\leq\; \sqrt{2 \pi \kappa^2 t} \, \norm{\varphi}\,,
\end{equation*}
which follows trivially from \eqref{Kato smoothing}.
Also, Lemma \ref{lemma: Zagatti} implies that $\norm{\tilde{\gamma}(s)} = \norm{\gamma}$. Thus we find that  
\begin{equation*}
I \;\leq\;
\sum_{\mathcal{Q} \in \mathscr{Q}(p,k,0)} \sum_{\mathcal{Q}' \in \mathscr{Q}(p+k, k',0)} \pbb{\frac{\pi \kappa^2 t}{2}}^{(k+1)/2} \p{2\pi \kappa^2 t}^{k'/2}
\norm{a^{(p)}} \norm{\gamma} \norm{\varphi}
\end{equation*}
Using the bound
\begin{equation*}
\abs{\mathscr{Q}(p,k,0)} \;\leq\; 4^p\,32^k\,,
\end{equation*}
which can be inferred from \eqref{bound on number of graphs}, we find
\begin{align*}
I &\;\leq\; 4^p \, 32^k \, 4^{p+k} 32^{k'} \pbb{\frac{\pi \kappa^2 t}{2}}^{(k+1)/2} \p{2\pi \kappa^2 t}^{k'/2}
\norm{a^{(p)}} \norm{\gamma} \norm{\varphi}
\\
&\;\leq\;
16^p \sqrt{2 \pi \kappa^2 t} \, \pb{64 \sqrt{2 \pi \kappa^2 t}}^k \, \pb{32 \sqrt{2 \pi \kappa^2 t}}^{k'} \, \norm{a^{(p)}} \norm{\gamma} \norm{\varphi}\,.
\end{align*}
Let $t < (2^{11} \pi \kappa^2)^{-1}$. Now Lemma \ref{lemma: bound on gamma_p} implies that $\norm{\gamma^{\otimes (p+k+k')} \Sigma_-^{(p+k)}}_1 \leq 1$. Using the inequality $\tr(A \Gamma) \leq \norm{A} \norm{\Gamma}_1$ we therefore find that \eqref{half of the error term} is bounded by
\begin{equation*}
16^p \sum_{k' = 0}^\infty \pb{32 \sqrt{2 \pi \kappa^2 t}}^k \, \pb{16 \sqrt{2 \pi \kappa^2 t}}^{k'} \, \norm{a^{(p)}} \norm{\gamma} 
\;=\; 16^p \frac{\pb{32 \sqrt{2 \pi \kappa^2 t}}^k}{1 - 16 \sqrt{2 \pi \kappa^2 t}} \, \norm{a^{(p)}} \norm{\gamma}\,.
\end{equation*}

The second term of $R_{ij}^k(t)$ is equal to the complex conjugate of the first. We thus arrive at the desired bound
\begin{equation} \label{bound on R}
\abs{R_{ij}^k(t)} \;\leq\; 2 \cdot 16^p \frac{\pb{32 \sqrt{2 \pi \kappa^2 t}}^k}{1 - 16 \sqrt{2 \pi \kappa^2 t}} \, \norm{a^{(p)}} \norm{\gamma}\,.
\end{equation}
Therefore the last line of \eqref{main expansion} is bounded by
\begin{multline*}
2 \cdot 16^p \frac{1}{1 - 16 \sqrt{2 \pi \kappa^2 t}} \, \norm{a^{(p)}} \norm{\gamma} \sum_{k = 0}^\infty (p+k)^2 \pb{32 \sqrt{2 \pi \kappa^2 t}}^k
\\
\leq\;
4 \cdot 16^p \, \ee^p \frac{1}{1 - 16 \sqrt{2 \pi \kappa^2 t}} \frac{1}{\pb{1 - 32 \sqrt{2 \pi \kappa^2 t}}^3} \, \norm{a^{(p)}} \norm{\gamma}\,,
\end{multline*}
where we used the estimate $\sum_{k = 0}^\infty (p+k)^L \, x^k \leq \frac{\me^p \, L!}{(1 - x)^{L+1}}$.

Next, we note that the second line of \eqref{main expansion}, i.e.\ the rest term, vanishes in the limit $K \to \infty$. The procedure is almost identical to (in fact easier than) the above estimation of $\abs{R_{ij}^k(t)}$. The result is
\begin{equation*}
\absbb{\int_{\Delta^K(t)} \dd \ul{t} \; \tr \p{G_{t, \ul{t}}^{(K,0)}(a^{(p)}) \, \tilde{\gamma}(t_K)^{\otimes (p+K)} \Sigma_-^{(p+K)}}}
\;\leq\; 2 \cdot 16^p \frac{\pb{32 \sqrt{2 \pi \kappa^2 t}}^K}{1 - 16 \sqrt{2 \pi \kappa^2 t}} \, \norm{a^{(p)}} \;\longrightarrow\; 0\,,
\end{equation*}
as $K \to \infty$.

Summarising, we have proven:
\begin{lemma} \label{lemma HF and tree terms}
Let $a^{(p)} \in \mathcal{B}(\mathcal{H}_-^{(p)})$. Then, for small times, 
\begin{equation*}
\absbb{\tr \pB{a_t^{(p)} \, \tilde{\gamma}(t)^{\otimes p} \, \Sigma_-^{(p)}} - \sum_{k = 0}^\infty \tr \pB{G_{t}^{(k,0)}(a^{(p)}) \, \gamma^{\otimes (p+k)} \Sigma_-^{(p+k)}}} \;\leq\; \norm{a^{(p)}} \, \norm{\gamma} \, C(p, \kappa,t)\,,
\end{equation*}
for some constant $C(p, \kappa, t)$.
\end{lemma}

\subsection{The mean-field limit}
We now have all the necessary ingredients to prove our main result.
Take some fixed orthonormal sequence $\Phi = (\varphi_i)_{i \in \N}$. Denote by $\Phi_N$ the truncated sequences $\Phi_N = (\varphi_1, \dots, \varphi_N)$, and let $\Phi_N(t)$ be the solution of the Hartree-Fock equation \eqref{rescaled integral Hartree-Fock} with initial data $\Phi_N$. 
The $N$-particle density matrix evolved with the Hartree-Fock dynamics is
\begin{equation*}
\tilde{\Gamma}_N(t) \;\deq\; \ketb{S(\Phi_N(t))} \brab{S(\Phi_N(t))}\,.
\end{equation*}
The $N$-particle density matrix evolved with the microscopic dynamics is
\begin{equation*}
\Gamma_N(t) \;\deq\; \ee^{- \ii t H_N} \ketb{S(\Phi_N)} \brab{S(\Phi_N)} \ee^{\ii t H_N}\,.
\end{equation*}
The $p$-particle marginals are defined by
\begin{equation*}
\Gamma_N^{(p)}(t) \;\deq\; \tr_{p+1 \dots N} \Gamma_N(t) \,,\qquad \tilde{\Gamma}_N^{(p)}(t) \;\deq\; \tr_{p+1 \dots N} \tilde{\Gamma}_N(t)\,.
\end{equation*}
The one-particle density matrix satisfying \eqref{integral Hartree-Fock for density matrices} is
\begin{equation*}
\gamma_N(t) \;\deq\; \tilde{\Gamma}^{(1)}_N(t)\,.
\end{equation*}
The quantities $\gamma_N^{(p)}(t) = \Sigma_-^{(p)} \gamma_N(t)^{\otimes p}$ and $\tilde{\Gamma}_N^{(p)}(t)$ are asymptotically equal:
\eqref{relation between marginals} implies that
\begin{equation} \label{asymptotic equality}
\normb{\gamma_N^{(p)}(t) - \tilde{\Gamma}_N^{(p)}(t)}_1 \;\leq\; \frac{p^2}{N}\,.
\end{equation}

Next, let $a^{(p)} \in \mathcal{B}(\mathcal{H}^{(p)}_-)$. From \eqref{expectation of second-quantised operator} and \eqref{one loop expansion} we get, for small times $t$, 
\begin{align*}
\tr \p{\widehat{\A}_N(a^{(p)}) \, \Gamma_N(t)}
&\;=\; \scalarb{S(\Phi_N)}{\me^{\mi t H_N} \, \widehat{\A}_N (a^{(p)}) \, \me^{- \mi t H_N} \, S(\Phi_N)}
\\
&\;=\; \sum_{k = 0}^\infty \scalarb{S(\Phi_N)}{\widehat{\A}_N \pb{G^{(k,0)}_t(a^{(p)})} \, S(\Phi_N)} + \scalarb{S(\Phi_N)}{L_N(t) \, S(\Phi_N)}\,,
\\
&\;=\;
\sum_{k = 0}^\infty \tr \p{G^{(k,0)}_t(a^{(p)}) \, \gamma_N^{(p+k)}}
 + \scalarb{S(\Phi_N)}{L_N(t) \, S(\Phi_N)}\,.
\end{align*}
Using Lemma \ref{lemma HF and tree terms} and \eqref{operator norm of gamma} we therefore get that, for small times,
\begin{equation*}
\absB{\tr \p{a_t^{(p)} \, \tilde{\gamma}_N^{(p)}(t)} - \tr \p{\widehat{\A}_N(a^{(p)}) \, \Gamma_N(t)}} \;\leq\; \frac{C(p,\kappa, t)}{N} \, \norm{a^{(p)}}\,.
\end{equation*}
Since the quantum-mechanical and the Hartree-Fock time-evolutions preserve the trace norm, we may iterate the above result, like in \cite{bosons}, to get: For all times $t \in \R$ we have that
\begin{equation*}
\absB{\tr \p{a^{(p)} \, \gamma_N^{(p)}(t)} - \tr \p{\widehat{\A}_N(a^{(p)}) \, \Gamma_N(t)}} \;\leq\; C(p, \kappa, t) \, \norm{a^{(p)}} \, f(N)\,,
\end{equation*}
with $\lim_{N \to \infty} f(N) = 0$.
Thus, by using the duality $\mathcal{B} \;=\; (\mathcal{L}^1)^*$ we find
\begin{theorem} \label{theorem: main result}
Let $p \in \N$ and $t \in \R$. Then
\begin{equation*}
\normb{\tilde{\Gamma}_N^{(p)}(t) - \Gamma_N^{(p)}(t)}_1 \;\longrightarrow\; 0\,, \qquad N \to \infty\,.
\end{equation*}
In particular, if $a^{(p)} \in \mathcal{B}(\mathcal{H}_-^{(p)})$, we have
\begin{equation*}
\scalarb{\ee^{-\ii t H_N} S(\Phi_N)}{a^{(p)} \otimes \umat^{(N-p)} \, \ee^{-\ii t H_N} S(\Phi_N)} - \scalarb{S(\Phi_N(t))}{a^{(p)} \otimes \umat^{(N-p)} \, S(\Phi_N(t))} \;\longrightarrow\; 0\,,
\end{equation*}
as $N \to \infty$.
\end{theorem}

\begin{remarks}
1.\ The limit $N \to \infty$ of $\Gamma^{(p)}_N(t)$ does not exist in $\norm{\cdot}_1$. Indeed, $\lim_{N \to \infty} \norm{\Gamma^{(p)}_N(t)} = 0$ but $\tr \Gamma^{(p)}_N(t) = 1$ (similarly for $\tilde{\Gamma}_N^{(p)}(t)$).
\\
2.\ As in \cite{bosons}, one can show that the function $f$ is a power law: $f(N) \sim N^{-\beta(t)}$, with $\beta(t) > 0$ for all $t$. However, $\beta(t) \to 0$ as $t \to \infty$. Our bound on the rate of convergence is therefore far from the expected optimal rate $\beta(t) = 1$, which we only obtain for short times.
\\
3.\ Although the exchange term $- \frac{1}{N} \sum_{j = 1}^N \p{w * (\varphi_i \bar{\varphi}_j)} \varphi_j$ is essential for our proof, it is not clear from our analysis whether it is of leading order as $N \to \infty$. The exchange term is known to be of subleading order in the scaling of \cite{NarnhoferSewell}, and hence in that case it does not play a role in the limiting dynamics (see \cite{ElgartErdosSchleinYau}).
\end{remarks}


\section{A Egorov-type result for small times} \label{section: Egorov}
In this section we describe how the many-body dynamics of fermions may be seen as the quantisation of a classical ``superhamiltonian'' system, whose dynamics is approximately described by the Hartree-Fock equation.

\subsection{A graded algebra of observables}
We start by defining a Grassmann algebra of anticommuting variables over the one-particle space $\mathcal{H} = \LL^2(\R^3)$, and equip it with a suitable norm. Formally, we consider the infinite-dimensional Grassmann algebra generated by $\{\ol{\psi(x)}, \psi(x)\}_{x \in \R^3}$. As it turns out, this algebra can be made into a Banach algebra under a natural choice of norm. This norm is most conveniently formulated by identifying elements of the Grassmann algebra with bounded operators between $\LL^2$-spaces.

Let
\begin{equation}
a \;=\; (a^{(p,q)})_{p,q \in \N} \,, \qquad a^{(p,q)} \in \mathcal{B}(\mathcal{H}^{(q)}_-; \mathcal{H}^{(p)}_-),
\end{equation}
be a family of bounded operators. Such objects will play the role of observables in the following.
By a slight abuse of notation we identify $a^{(p,q)}$ with the family obtained by adjoining zeroes to it.

Define
\begin{equation*}
\mathfrak{B}^G \;\deq\; \{a = (a^{(p,q)}) \,:\, a^{(p,q)} = 0 \text{ for all but finitely many } (p,q)\}\,.
\end{equation*}
We introduce a norm on $\mathfrak{B}^G$ through
\begin{equation}
\norm{a}_{\mathfrak{B}^G} \;\deq\; \sum_{p,q \in \N} \norm{a^{(p,q)}}\,,
\end{equation}
and define $\ol{\mathfrak{B}}^G$ as the completion of $\mathfrak{B}^G$.

We also introduce a multiplication on $\ol{\mathfrak{B}}^G$ defined by
\begin{equation}
(ab)^{(p,q)} \;\deq\; \sum_{\substack{p_1 + p_2 = p \\ q_1 + q_2 = q}} (-1)^{p_2(p_1 + q_1)} \, P_- (a^{(p_1, q_1)} \otimes b^{(p_2, q_2)}) P_-\,.
\end{equation}
The seemingly odd choice of sign will soon become clear. It is now easy to check that $\ol{\mathfrak{B}}^G$ is an associative Banach algebra with identity
\begin{equation*}
\umat^{(p,q)} \;=\; \delta_{p0} \, \delta_{q0}\,.
\end{equation*}
Note that $\ol{\mathfrak{B}}^G$ bears a $\Z$-grading, with degree map 
\begin{equation*}
\deg a^{(p,q)} \;\deq\; p-q\,. 
\end{equation*}
An observable is gauge invariant when its degree is equal to $0$.
One readily sees that
\begin{equation*}
ab = (-1)^{\deg a \, \deg b} ba\,.
\end{equation*}


We now identify $\ol{\mathfrak{B}}^G$ with a Grassmann algebra of anticommuting variables. For $f \in \mathcal{H}$ define $\psi(f) \in \mathcal{B}(\mathcal{H}; \C) \subset \ol{\mathfrak{B}}^G$ through
\begin{equation}
\psi(f) \, g \;\deq\; \scalar{f}{g}
\end{equation}
and $\ol{\psi}(f) \in \mathcal{B}(\C; \mathcal{H}) \subset \ol{\mathfrak{B}}^G$ through
\begin{equation}
\ol{\psi}(f) z \;\deq\; f z\,.
\end{equation}
We may now consider arbitrary polynomials in the variables $\{\ol{\psi}(f), \psi(f) \,:\, f \in \mathcal{H}\}$. It is a simple matter to check that
\begin{equation*}
\psi(f) \, \psi(g) + \psi(g) \, \psi(f) \;=\; \psi(f) \ol{\psi}(g) + \ol{\psi}(g) \, \psi(f) \;=\; \ol{\psi}(f) \, \ol{\psi}(g) + \ol{\psi}(g) \, \ol{\psi}(f) \;=\; 0\,, 
\end{equation*}
for all $f,g \in \mathcal{H}$. Furthermore, we have that
\begin{equation} \label{polynomial in psis}
\ol{\psi}(f_p) \cdots \ol{\psi}(f_1) \psi(g_1) \cdots \psi(g_q) \;=\; P_-^{(p)} \, \ket{f_1 \otimes \cdots \otimes f_p} \bra{g_1 \otimes \cdots \otimes g_q} \, P_-^{(q)}\,. 
\end{equation}
Linear combinations of expressions of the form \eqref{polynomial in psis} are dense in $\ol{\mathfrak{B}}^G$.
It is often convenient to write a family $a$ of bounded operators using the ``Grassmann generators'' $\{\ol{\psi}, \psi\}$. To this end we set
\begin{equation*}
\psi(x) \;\deq\; \psi(\delta_x) \,,\qquad \ol{\psi}(x) \;\deq\; \ol{\psi}(\delta_x)\,,
\end{equation*}
where $\delta_x$ is Dirac's delta mass at $x$.
Expressions of the form \eqref{polynomial in psis} are now understood as densely defined quadratic forms. One immediately finds
\begin{multline}
a \;\equiv\; \A^{\! G}(a) \;\deq\; \sum_{p,q} \int \md x_1 \dots \md x_p \, \md y_1 \dots \md y_q
\\
{}\times{} \ol{\psi}(x_p) \cdots \ol{\psi}(x_1) \, a^{(p,q)}(x_1, \dots, x_p; y_1, \dots, y_q) \, \psi(y_1) \cdots \psi(y_q)\,.
\end{multline}
We use the notation $\A^{\! G}(a)$ to emphasize that the family $a$ is represented using Grassmann generators.

\subsection{A graded Poisson bracket}
Next, we note that $\mathfrak{B}^G$ carries the graded Poisson bracket
\begin{equation}
\{a, b\} \;\deq\; \mi \int \md x \; \qbb{a \, \frac{\overset{\leftarrow}{\delta}}{\delta \ol{\psi}(x)}
\frac{\overset{\rightarrow}{\delta}}{\delta \psi(x)} \, b
+
a \, \frac{\overset{\leftarrow}{\delta}}{\delta \psi(x)}
\frac{\overset{\rightarrow}{\delta}}{\delta \ol{\psi}(x)} \, b}\,,
\end{equation}
where $a,b \in \mathfrak{B}^G$. Here we use the usual conventions for derivatives with respect to Grassmann variables (see e.g.\ \cite{Salmhofer}, Appendix B).
In terms of kernels the graded Poisson bracket can be expressed as
\begin{equation}
\{\psi(x), \ol{\psi}(y)\} \;=\; \mi \delta(x-y) \, \qquad \{\psi(x), \psi(y)\} \;=\; \{\ol{\psi}(x), \ol{\psi}(y)\} \;=\; 0\,.
\end{equation}
We now list the important properties of the graded Poisson bracket.
\begin{itemize}
\item[(i)]
$\{a, b\} \;=\; (-1)^{1 + \deg a \, \deg b} \{b,a\}$\,.
\item[(ii)]
$(-1)^{\deg b \, (\deg a +  \deg c)} \{a,\{b, c\}\} + \text{cyclic permutations} \;=\; 0$\,.
\item[(iii)]
$\{a, bc\} \;=\; \{a, b\} c + (-1)^{\deg a \, \deg b} b \{a,c\}$\,.
\end{itemize}
\begin{proof}
Let us start with (i):
\begin{align*}
\{a,b\} &\;=\; \mi \int \md x \, \qBB{(-1)^{\deg a + 1}\frac{\delta a}{\delta \ol{\psi}(x)} \frac{\delta b}{\delta \psi(x)} + (-1)^{\deg a + 1} \frac{\delta a}{\delta \psi(x)} \frac{\delta b}{\delta \ol{\psi}(x)}}
\\
&\;=\; \mi \int \md x \, \qBB{(-1)^{\deg a \, \deg b + \deg b} \frac{\delta b}{\delta \psi(x)} \frac{\delta a}{\delta \ol{\psi}(x)} + (-1)^{\deg a \, \deg b + \deg b} \frac{\delta b}{\delta \ol{\psi}(x)} \frac{\delta a}{\delta \psi(x)}}
\\
&\;=\; (-1)^{1 + \deg a\, \deg b} \{b,a\}\,.
\end{align*}
In order to show (ii), we note that the left-hand side can be written as a sum of three terms, the first of which contains second derivatives of $a$, the second second derivatives of $b$ and the third second derivatives of $c$. Let us consider the third one. It is equal to the terms containing second derivatives of $c$ of
\begin{multline*}
(-1)^{\deg b (\deg a + \deg c)} \{a, \{b,c\}\} + (-1)^{\deg c (\deg b + \deg a)} \{b, \{c,a\}\}
\\
=\; (-1)^{\deg b (\deg a + \deg c)} \{a, \{b,c\}\} + (-1)^{\deg c (\deg b + \deg a) + 1 + \deg c \, \deg a} \{b, \{a,c\}\}\,,
\end{multline*}
where (i) was used. Define the derivation $L_a b \deq \{a,b\}$\,. Thus we need to compute the terms containing second derivatives of $c$ of
\begin{equation*}
(-1)^{\deg a \, \deg b + \deg b \, \deg c} L_a L_b c - (-1)^{\deg b \, \deg c} L_b L_a c\,.
\end{equation*}
Since we are only considering terms containing second derivatives of $c$, both derivations $L_a$ and $L_b$ must act only on $c$, and one finds
\begin{equation*}
(-1)^{\deg a \, \deg b + \deg b \, \deg c} L_a L_b c - (-1)^{\deg a \, \deg b + \deg b \, \deg c} L_a L_b c \;=\; 0\,.
\end{equation*}
We omit the straightforward proof of (iii).
\end{proof}
Furthermore, one finds by explicit calculation
\begin{multline} \label{graded Poisson bracket computed}
\{\A^{\! G}(a^{(p_1,q_1)}), \A^{\! G}(b^{(p_2,q_2)})\} 
\\
\;=\; \mi \, (-1)^{(p_2 + 1)(p_1 + q_1)}  q_1 p_2 \, \A^{\! G} \qB{\pb{a^{(p_1, q_1)} \otimes \umat^{(p_2 - 1)}} \, \pb{\umat^{(q_1 - 1)} \otimes b^{(p_2, q_2)}}} 
\\
{}-{}\mi \, (-1)^{(q_1 + 1)(p_2 + q_2)} p_1 q_2 \, \A^{\! G} \qB{\pb{b^{(p_2, q_2)} \otimes \umat^{(p_1 - 1)}} \, \pb{\umat^{(q_2 - 1)} \otimes a^{(p_1, q_1)}}}\,.
\end{multline}
%Now \eqref{graded Poisson bracket computed} implies
%\begin{equation}
%\{\A^{\! G}(a^{(p)}), \A^{\! G}(b^{(q)})\} \;=\; \mi pq \A^{\! G} \pb{[a^{(p)}, b^{(q)}]_1}\,.
%\end{equation}

\subsection{States}
To the algebra of observables $(\ol{\mathfrak{B}}^G, \norm{\cdot}_{\mathfrak{B}^G})$ is associated the space of states
\begin{equation*}
\mathfrak{R} \;\deq\; \pb{\ol{\mathfrak{B}}^G}^*\,.
\end{equation*}
By using the standard argument of the proof that $(l^1)^* = l^\infty$ (see e.g.\ \cite{ReedSimonI}), one finds that
\begin{equation*}
\mathfrak{R} \;=\; \hb{\rho = (\rho^{(p,q)})_{p,q \in \N} \,:\, \rho^{(p,q)} \in \mathcal{B}(\mathcal{H}_-^{(q)}; \mathcal{H}_-^{(p)})\,, \; \norm{\rho}_{\mathfrak{R}} < \infty}\,,
\end{equation*}
where
\begin{equation*}
\norm{\rho}_{\mathfrak{R}} \;\deq\; \sup_{p,q \in \N} \norm{\rho^{(p,q)}}_1
\end{equation*}
and
\begin{equation*}
\norm{\rho^{(p,q)}}_1 \;\deq\; \sup \hb{\absb{\tr (\rho^{(p,q)} a^{(q,p)})} \,:\, a^{(q,p)} \in \mathcal{B}(\mathcal{H}_-^{(p)}, \mathcal{H}_-^{(q)}) \,,\; \norm{a^{(q,p)}} \leq 1}\,.
\end{equation*}
Note that if $p = q$ then $\norm{\cdot}_1$ is the usual trace norm. The dual action is given by
\begin{equation*}
\scalar{\rho}{a} \;\deq\; \sum_{p,q \in \N} \tr (\rho^{(p,q)} a^{(q,p)})\,.
\end{equation*}
We abbreviate $\rho^{(p,p)} \equiv \rho^{(p)}$ in the case of gauge invariant states. Next, we note that \eqref{polynomial in psis} implies that the operator kernel of $\rho^{(p,q)}$ may be expressed as
\begin{equation} \label{operator kernel of rho}
\rho^{(p,q)}(x_1, \dots, x_p;y_1, \dots, y_q) \;=\; \scalarb{\rho}{\ol{\psi}(y_q) \cdots \ol{\psi}(y_1) \psi(x_1) \cdots \psi(x_p)}\,.
\end{equation}

There is a particular subset of gauge invariant states that is of interest for studying the Hartree-Fock dynamics. Let $\gamma \in \mathcal{D}$ be a one-particle density matrix. Define the state $\rho_\gamma$ through $\rho_\gamma^{(p,q)} = 0$ if $p \neq q$ and
\begin{equation} \label{general state from density matrix}
\rho_\gamma^{(p)} \;\deq\; \gamma^{(p)}\,,
\end{equation}
where $\gamma^{(p)}$ is defined in \eqref{p-particle action of quasi-free state}. One immediately finds $\norm{\rho_\gamma}_1 = \norm{\gamma}_1$. %The set $\{\rho_\gamma \,:\, \gamma \in \mathcal{D}\} \;\subset\; \mathfrak{R}$ is a cone of positive elements in the sense that $\scalar{\rho_\gamma}{a} \geq 0$ whenever $a^{(p)} \geq 0$ for all $p$.

\subsection{Hamilton function and dynamics}
Let $h$ be the one-particle Hamiltonian and $w$ the two-body interaction potential. We define a Hamilton function on (a dense %\footnote{For instance: Define $S \deq (\umat - \Delta)^{1/2}$, and consider the space $\{\rho \,:\, \norm{S \rho^{(1)} S}_1 < \infty\,,\; \norm{S_1 S_2 \rho^{(2)} S_2 S_1}_1 < \infty\}$. The one sees immediately that $\abs{\tr(\rho^{(1)}h)} < \infty$. Furthermore,
%\begin{equation*}
%\abs{\tr(\rho^{(2)} w)} \;=\; \absb{\tr\pb{S_1 S_2 \rho^{(2)} S_2 S_1 S_1^{-1} S_2^{-1} w S_2^{-1} S_1^{-1}}}
%\;\leq\; \normb{S_1 S_2 \rho^{(2)} S_2 S_1}_1 \, \normb{S_1^{-1} S_2^{-1} w S_2^{-1} S_1^{-1}}\,.
%\end{equation*}
%Then $S_1 S_2 \geq (\umat - \Delta_{12})^{1/2}$ and the Hardy-Littlewood-Sobolev and Sobolev inequalities imply that $ \normb{S_1^{-1} S_2^{-1} w S_2^{-1} S_1^{-1}} < \infty$ potentials $w$ with Coulomb-type singularities.}
subset of) the phase space $\mathfrak{R}$ through
\begin{align}
H &\;\deq\; \A^{\! G}(h) + \frac{1}{2} \A^{\! G}(w) 
\notag \\ \label{Grassmann Hamiltonian}
&\;=\; \int \md x \, \md y \; \ol{\psi}(x) \, h(x;y) \, \psi(y) + \frac{1}{2} \int \md x \, \md y \; \ol{\psi}(y) \ol{\psi}(x) \, w(w - y) \, \psi(x) \psi(y)\,.
\end{align}

The Hamiltonian equations of motion read
\begin{equation*}
\dot{a} \;=\; \{H, a\}\,,
\end{equation*}
where $a \in \mathfrak{B}^G$.
Instead of the ``Heisenberg'' evolution of $a$ we consider the dual ``Schr\"odinger'' evolution of states:
\begin{equation*}
\scalar{\rho(t)}{a} \;\deq\; \scalar{\rho}{a(t)}\,.
\end{equation*}
The equation of motion for states reads
\begin{multline} \label{equation of motion for grassmann states}
\mi \partial_t \rho^{(p,q)}(x_1, \dots, x_p; y_1, \dots, y_q) \;=\; \pBB{\sum_{i = 1}^p h_{x_i} - \sum_{i = 1}^q h_{y_i}} \rho^{(p,q)}(x_1, \dots, x_p; y_1, \dots, y_q)
\\
+ \int \md u \; \pBB{\sum_{i = 1}^p w(u - x_i) - \sum_{i = 1}^q w(u - y_i)} \rho^{(p+1,q+1)}(x_1, \dots, x_p, u; y_1, \dots, y_q, u)\,.
\end{multline}
This has the form of an infinite hierarchy, which decouples over subspaces of different degree.
In order to show \eqref{equation of motion for grassmann states} we compute 
\begin{multline*}
\mi \hb{H, \ol{\psi}(y_q) \cdots \ol{\psi}(y_1) \psi(x_1) \cdots \psi(x_p)} \;=\; \pBB{\sum_{i = 1}^p h_{x_i} - \sum_{i = 1}^q h_{y_i}} \, \ol{\psi}(y_q) \cdots \ol{\psi}(y_1) \psi(x_1) \cdots \psi(x_p)
\\
{}+{} \int \md u \pBB{\sum_{i = 1}^p w(u - x_i) - \sum_{i = 1}^q w(u - y_i)}
\ol{\psi}(u) \ol{\psi}(y_q) \cdots \ol{\psi}(y_1) \psi(x_1) \cdots \psi(x_p) \psi(u)\,.
\end{multline*}
Then \eqref{equation of motion for grassmann states} follows from \eqref{operator kernel of rho} and
\begin{align*}
\mi \partial_t \rho^{(p,q)}(x_1, \dots, x_p;y_1, \dots, y_q) &\;=\; \mi \partial_t \scalarb{\rho}{\ol{\psi}(y_q) \cdots \ol{\psi}(y_1) \psi(x_1) \cdots \psi(x_p)}
\\
&\;=\; \scalarb{\rho}{\mi \hb{H, \ol{\psi}(y_q) \cdots \ol{\psi}(y_1) \psi(x_1) \cdots \psi(x_p)}}\,.
\end{align*}

Next, we outline how to solve the equation of motion \eqref{equation of motion for grassmann states}. Let us first rewrite it as
\begin{multline*}
\ii \partial_t \rho^{(p,q)} \;=\; \sum_{i = 1}^p h_i \rho^{(p,q)} - \sum_{i = 1}^q \rho^{(p,q)} h_i
\\
{}+{} \sum_{i = 1}^p \tr_{p+1, q+1} \pb{W_{i \, p+1} \rho^{(p+1, q+1)}} - \sum_{i = 1}^q \tr_{p+1, q+1} \pb{\rho^{(p+1, q+1)} W_{i \, q+1}}\,.
\end{multline*}
We may now proceed exactly as with the density matrix Hartree-Fock equation \eqref{Hartree-Fock for density matrices}, i.e.\ express it as an integral equation in the interaction picture. This yields a tree expansion for the quantity $\tr\pb{\rho^{(p,q)}(t) \, a^{(q,p)}}$, where $\rho(0) \in \mathfrak{R}$. We omit the uninteresting details. As above, the tree expansion converges if $t < T$, where
\begin{equation} \label{definition of T}
T \;\deq\; (2^{11} \pi \kappa^2)^{-1}\,.
\end{equation}
Unfortunately, the time evolution \eqref{equation of motion for grassmann states} does not preserve the norm of $\rho$, which means that we cannot iterate the short-time result.

From now on, we only consider gauge invariant invariant quantities. Take some gauge invariant state $\rho = (\rho^{(p)})_{p \in \N} \in \mathfrak{R}$. For simplicity, we assume that the sequence $\rho$ is finite (as is the case if $\rho$ is defined by a Slater determinant, see below). Let us denote the Hamiltonian flow on $\mathfrak{R}$ by $\phi^t$. We have seen that $\phi^t$ is well-defined by its tree expansion for $t < T$. The solution of \eqref{equation of motion for grassmann states} with initial data $\rho$, $\rho(t) = \phi^t(\rho)$, satisfies the equation
\begin{equation} \label{grassman state duhamel}
\tilde{\rho}^{(p)}(t) \;=\; \rho^{(p)} - \ii \int_0^t \dd s \; \sum_{i = 1}^p\tr_{p+1} \qb{W_{i\,p+1,s}, \tilde{\rho}^{(p+1)}(s)}\,,
\end{equation}
where $\tilde{\rho}^{(p)}(t) \deq \ee^{\ii \sum_i h_i t} \, \rho^{(p)}(t) \, \ee^{-\ii \sum_i h_i t}$. Let us take a gauge invariant observable $a^{(p,p)} \equiv a^{(p)} \in \mathfrak{A}^G$, where
\begin{equation*}
\mathfrak{A}^G \;\deq\; \{a \in \mathfrak{B}^G \,:\, a^{(p,q)} = 0 \text{ if } p \neq q\}
\end{equation*}
is the set of gauge invariant observables. Then \eqref{grassman state duhamel} implies
\begin{multline*}
\tr \pb{a^{(p)} \rho^{(p)}(t)} \;=\; \tr \pb{a^{(p)}_t \tilde{\rho}^{(p)}(t)} 
\\
=\; \tr \pb{a^{(p)}_t \rho^{(p)}} + \ii \int_0^t \dd s \; \sum_{i = 1}^p \tr \pB{\qB{W_{i\,p+1, s}, a^{(p)}_t \otimes \umat} \tilde{\rho}^{(p+1)}(s)}\,.
\end{multline*}
Iteration of this identity gives
\begin{equation*}
\tr \pb{a^{(p)} \rho^{(p)}(t)} \;=\; \sum_{k = 0}^\infty \tr \pB{G_t^{(k,0)}(a^{(p)}) \, \rho^{(p+k)}}\,.
\end{equation*}
Summarising:
\begin{multline*}
\scalarb{a^{(p)} \circ \phi^t}{\rho} \;=\; \scalarb{a^{(p)}}{\rho(t)} \;=\; \tr \pb{a^{(p)} \rho^{(p)}(t)} 
\\
=\; \sum_{k = 0}^\infty \tr \pB{G_t^{(k,0)}(a^{(p)}) \, \rho^{(p+k)}} \;=\; \scalarBB{\sum_{k = 0}^\infty G_t^{(k,0)}(a^{(p)})}{\rho}\,.
\end{multline*}
This series converges for $t < T$, uniformly for bounded $\norm{a^{(p)}}_{\mathfrak{B}^G}$ and $\norm{\rho}_{\mathfrak{R}}$. Therefore we get the norm-convergent series
\begin{equation}  \label{tree expansion for superhamiltonian flow}
\A^G(a^{(p)}) \circ \phi^t \;=\; \sum_{k = 0}^\infty \A^G \pB{G_t^{(k,0)}(a^{(p)})}\,,
\end{equation}
provided that $t < T$.

Finally, we discuss the relationship between the Hartree-Fock dynamics and the dynamics generated by \eqref{equation of motion for grassmann states}. Take a density matrix $\gamma \in \mathcal{D}$ and consider the state $\rho = \rho_\gamma$ defined in \eqref{general state from density matrix}. If one chooses a sequence $\gamma_N$ such that $\norm{\gamma_N} \to 0$ as $N \to \infty$ (e.g.\ a sequence of Slater determinants), then Lemma \ref{lemma HF and tree terms} implies that \eqref{equation of motion for grassmann states} and the Hartree-Fock equation describe the same dynamics for large $N$.

\subsection{Quantisation and a Egorov-type theorem}
In this final section we introduce a Wick quantisation of the above ``superhamiltonian'' system and formulate the mean-field limit as a Egorov-type theorem. From now on, we use $n$ to denote the number of particles, and $N$ to denote a free parameter (the inverse ``deformation parameter'' of the quantisation). Ultimately, we shall choose $n = N$. The underlying ideas were described in detail in \cite{bosons}; here we merely state how they apply in the current setting.

Consider the fermionic Fock space
\begin{equation*}
\mathcal{F}_- \;=\; \bigoplus_{n \in \N} \mathcal{H}^{(n)}_-\,,
\end{equation*}
on which act the creation and annihilation operators $a^*(x), a(x)$ (see e.g.\ \cite{BratteliRobinson}). We define the rescaled operators as $\widehat{\psi}^\#_N(x) \deq a^\#(x)/\sqrt{N}$, where $\# \in \{*, \text{nothing}\}$. They satisfy the canonical anticommutation relations
\begin{equation*}
\qB{\widehat{\psi}_N(x), \widehat{\psi}_N(y)}_+ \;=\; \qB{\widehat{\psi}_N^*(x), \widehat{\psi}_N^*(y)}_+ \;=\; 0\,, \qquad 
\qB{\widehat{\psi}_N(x), \widehat{\psi}_N^*(y)}_+ \;=\; \frac{1}{N} \delta(x-y)\,,
\end{equation*}
where $[A,B]_+ \deq AB + BA$ is the anticommutator.

Let $a^{(p)} \in \mathcal{B}(\mathcal{H}^{(p)}_-)$ with distribution kernel $a^{(p)}(x_1, \dots, x_p;y_1, \dots, y_p)$. The second quantisation of $a^{(p)}$ is defined by
\begin{multline*}
\widehat{\A}^G_N(a^{(p)}) \;\deq\; \int \md x_1 \cdots \md x_p \, \md y_1 \cdots \md y_p 
\\
\widehat{\psi}_N^*(x_p) \cdots \widehat{\psi}_N^*(x_1) \, a^{(p)}(x_1, \dots, x_p;y_1, \dots, y_p) \, 
\widehat{\psi}_N(y_1) \cdots \widehat{\psi}_N(y_p)\,.
\end{multline*}
This defines a closed operator on $\mathcal{F}_-$. Moreover, restricted to the subspace $\mathcal{H}^{(N)}_-$, $\widehat{\A}^G_N(a^{(p)})$ is equal to $\widehat{\A}_N(a^{(p)})$ defined in \eqref{second quantisation of a}. Additional properties of the operation $\widehat{\A}^G_N$ are listed in \cite{bosons}.

Let $\widehat{\mathfrak{A}}^G$ be the linear hull of
\begin{equation*}
\hb{\widehat{\A}^G_N(a^{(p)}) \,:\, a^{(p)} \in \mathcal{B}^{(p)}\,,\, p \in \N}\,,
\end{equation*}
We define \emph{quantisation} as the linear map $\widehat{(\cdot)}_N \,:\, \mathfrak{A}^G \to \widehat{\mathfrak{A}}^G$ defined by the formal replacement $\psi(x) \mapsto \widehat{\psi}_N(x)$ and $\bar{\psi}(x) \mapsto \widehat{\psi}^*_N(x)$ followed by Wick ordering. In other words,
\begin{equation*}
\widehat{(\cdot)}_N \,:\, \A^{\! G}(a^{(p)}) \;\mapsto\; \widehat{\A}^G_N(a^{(p)})\,.
\end{equation*}
Using \eqref{graded Poisson bracket computed}, it is easy to see that, for $A,B \in \mathfrak{A}^G$,
\begin{equation*}
\qb{\widehat{A}_N, \widehat{B}_N}_+ \;=\; \frac{N^{-1}}{\ii} \widehat{\h{A, B}}_N + O(N^{-2})\,.
\end{equation*}
This identifies $N^{-1}$ as the deformation parameter of $\widehat{(\cdot)}_N$.

Extending the definition of $\widehat{(\cdot)}_N$ to unbounded operators in the obvious way, we define a Hamiltonian $\widehat{H}_N$ on $\mathcal{F}_-$ as the quantisation of the Hamilton function $H$ defined in \eqref{Grassmann Hamiltonian}. When restricted to $\mathcal{H}^{(N)}_-$, $N \widehat{H}_N$ is equal to the Hamiltonian with mean-field scaling \eqref{mean-field Hamiltonian}.

Now \eqref{one loop expansion}, \eqref{estimate on loop terms} and \eqref{tree expansion for superhamiltonian flow} yield the following Egorov-type theorem.
\begin{theorem}
Let $A \in \mathfrak{A}^G$ and $t < T$, with $T$ defined in \eqref{definition of T}. Then
\begin{equation*}
\normbb{\pbb{\ee^{\ii t N \widehat{H}_N} \, \widehat{A}_N\, \ee^{-\ii t N \widehat{H}_N}  - \widehat{(A \circ \phi^t)}_N}\biggr|_{\mathcal{H}^{(N)}_-}} \;\leq\; \frac{C}{N}\,,
\end{equation*}
for some $C > 0$.
\end{theorem}

\appendix
\section{Hamiltonian formulation for density matrices}
In Section \ref{section: Hartree-Fock equation}, we chose a Hamiltonian formulation of the Hartree-Fock equation \eqref{Hartree-Fock} in terms of sequences of orbitals. Alternatively, we could just as well have used a Hamiltonian formulation in terms of density matrices. To see how the density matrix Hartree-Fock equation \eqref{Hartree-Fock for density matrices} can be written as a Hamiltonian equation of motion of a classical Hamiltonian system, consider the Hilbert space
\begin{equation*}
\widehat{\mathcal{H}} \;=\; \mathcal{L}^2(\mathcal{H})\,,
\end{equation*}
the space of Hilbert-Schmidt operators, with scalar product
\begin{equation*}
\scalar{\kappa}{\rho} \;\deq\; \tr (\kappa^* \rho)\,.
\end{equation*}
We write the density matrix $\gamma \in \mathcal{D}$ as $\gamma = \kappa \kappa^*$, where $\kappa \in \widehat{\mathcal{H}}$.
The classical phase space is then given by a Sobolev-type space of Hilbert-Schmidt operators
\begin{equation*}
\widehat{\Gamma} \;\deq\; \hb{\kappa \in \widehat{\mathcal{H}} \,:\, \tr (\kappa^* (\umat - \Delta) \kappa) < \infty}\,.
\end{equation*}
We define polynomial functions on $\widehat{\Gamma}$ through
\begin{equation*}
\B(a^{(p)})(\kappa) \;\deq\; \scalarb{\kappa^{\otimes p}}{a^{(p)} \kappa^{\otimes p}}\,,
\end{equation*}
where $a^{(p)} \in \mathcal{B}(\mathcal{H}^{\otimes p})$. 

The affine space $\widehat{\Gamma}$ carries a Symplectic form defined by
\begin{equation*}
\omega \;=\; -\ii \int \dd x \, \dd y \; \dd \bar{\kappa}(x,y) \wedge \dd \kappa(x,y)\,,
\end{equation*}
where $\kappa(x,y) \in \LL^2(\R^3 \times \R^3)$ is the operator kernel of $\kappa$. The Poisson bracket then reads
\begin{align*}
\h{\kappa^\#(x,y), \kappa^\#(x',y')} &\;=\; 0
\\
\h{\kappa(x,y), \bar{\kappa}(x',y')} &\;=\; - \ii \delta(x - x') \delta(y - y')\,.
\end{align*}
The Hamilton function is defined by
\begin{equation*}
H \;\deq\; \B(h) + \frac{1}{2} \B(\mathcal{W})\,.
\end{equation*}
By using Sobolev-type inequalities one readily sees that $H$ is well-defined on $\widehat{\Gamma}$.
After a short computation, one finds that the Hamiltonian equation of motion,
\begin{equation*}
\ii \partial_t \kappa(x,y) \;=\; \frac{\delta H}{\delta \ol{\kappa}(x,y)} \;=\; \ii \, \{H, \kappa(x,y)\}\,,
\end{equation*}
reads
\begin{equation*}
\ii \partial_t \kappa \;=\; h \kappa + \tr_2 \pb{\mathcal{W} \, \kappa \otimes (\kappa \kappa^*)}\,.
\end{equation*}
It follows that $\gamma = \kappa \kappa^*$ satisfies \eqref{Hartree-Fock for density matrices}.

\begin{thebibliography}{99}

\bibitem{Bach}
V.\ Bach, \emph{Error bound for the Hartree-Fock energy of atoms and molecules}, Commun.\ Math.\ Phys. \textbf{147} (1992), no.\ 3, 527--548. 

\bibitem{Bove}
A.\ Bove, G.\ Da Prato and G.\ Fano, \emph{An existence proof for the Hartree-Fock time-dependent problem with bounded two-body interaction}, Commun.\ Math.\ Phys.\ \textbf{37} (1974), 183--191.

\bibitem{BardosGolse}
C.\ Bardos, F.\ Golse, A.\ D.\ Gottlieb and J.\ Mauser, \emph{Mean-field dynamics of fermions and the time-dependent Hartree-Fock equation},
J.\ Math.\ Pures et Appl.\ \textbf{82} (2003), no.\ 6, 665--683.

\bibitem{BratteliRobinson}
O.~Bratteli and D.~W.~Robinson, \emph{Operator Algebras and Quantum Statistical Mechanics 2}, Springer, 2002.

\bibitem{Chadam}
J.~M.~Chadam, \emph{The Time Dependent Hartree-Fock Equations with Coulomb Two-Body Interaction}, Commun.\ Math.\ Phys.\ \textbf{46} (1976), 99--104.

\bibitem{ElgartErdosSchleinYau}
A.\ Elgart, L.\ Erd\H{o}s, B.\ Schlein and H.-T.\ Yau, \emph{Nonlinear Hartree equation as the mean fiald limit of weakly coupled fermions}, J.\ Math.\ Pures et Appl.\ \textbf{83} (2004), 1241--1273.

\bibitem{Fefferman1}
C.\ Fefferman and L.\ Seco, \emph{On the energy of a large atom}, Bull.\ Amer.\ Soc.\ (N.S.) \textbf{23} (1990), no.\ 2, 525--530.

\bibitem{Fefferman2}
C.\ Fefferman and L.\ Seco, \emph{On the Dirac and Schwinger corrections to the ground energy of an atom}, Adv.\ Math.\  \textbf{107} (1994), no.\ 1, 1--185.

\bibitem{bosons}
J.\ Fr\"ohlich, A.\ Knowles and S.\ Schwarz, \emph{On the mean-field limit of bosons with Coulomb two-body interaction}, math-ph/0805.4299v1.

\bibitem{GrafSolovej}
G.\ M.\ Graf and J.\ P.\ Solovej, \emph{A correlation estimate with applications to quantum systems with Coulomb interactions}, Rev.\ Math.\ Phys.\ \textbf{6} (1994), no.\ 5a, 977-997.

\bibitem{LiebLoss}
E.\ H.\ Lieb and M.\ Loss, \emph{Analysis}, American Mathematical Society, 2001.

\bibitem{LiebSimon1}
E.\ H.\ Lieb and B.\ Simon, \emph{The Hartree-Fock theory for Coulomb systems}, Commun.\ Math.\ Phys.\ \textbf{53} (1977), no.\ 3, 185--194.

\bibitem{LiebSimon2}
E.\ H.\ Lieb and B.\ Simon, \emph{The Thomas-Fermi theory of atoms, molecules and solids}, Adv.\ in Math.\ \textbf{23} (1977), no.\ 1, 22--116.

\bibitem{NarnhoferSewell}
H.\ Narnhofer and G.\ L.\ Sewell, \emph{Vlasov hydrodynamics of a quantum mechanical model}, Commun.\ Math.\ Phys.\ \textbf{79} (1981), 9--24.

\bibitem{ReedSimonI}
M.~Reed and B.~Simon, \emph{Methods of Modern Mathematical Physics I: Functional Analysis}, Academic Press, 1980.

\bibitem{ReedSimonII}
M.~Reed and B.~Simon, \emph{Methods of Modern Mathematical Physics II: Fourier Analysis, Self-Adjointness}, Academic Press, 1975.

\bibitem{Salmhofer}
M.\ Salmhofer, \emph{Renormalization, An Introduction}, Springer, 1999.

\bibitem{Spohn}
H.\ Spohn, \emph{On the Vlasov hierarchy}, Math.\ Meth.\ Appl.\ Sci.\ \textbf{3} (1981), no.\ 4, 445--455.

\bibitem{Zagatti}
S.~Zagatti, \emph{The Cauchy problem for Hartree-Fock time-dependent equations}, Ann.\ Inst.\ Henri~Poincar\'e (A) \textbf{56} (1992), 357--374.

\end{thebibliography}
\end{document}
---------------0810240337487--
