\documentclass[reqno]{amsart} \usepackage{graphicx} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2009(2009), No. 12, pp. 1--22.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2009 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2009/12\hfil Fully discrete Galerkin schemes] {Fully discrete Galerkin schemes for the nonlinear and nonlocal Hartree equation} \author[W. H. Aschbacher\hfil EJDE-2009/12\hfilneg] {Walter H. Aschbacher} \address{Walter H. Aschbacher \newline Technische Universit\"at M\"unchen \\ Zentrum Mathematik, M5\\ 85747 Garching, Germany} \email{aschbacher@ma.tum.de} \thanks{Submitted September 26, 2008. Published January 12, 2009.} \subjclass[2000]{35Q40, 35Q35, 35J60, 65M60, 65N30} \keywords{Hartree equation; quantum many-body system; \hfill\break\indent weakly nonlinear dispersive waves; Newtonian gravity; Galerkin theory; \hfill\break\indent finite element methods; discretization accuracy} \begin{abstract} We study the time dependent Hartree equation in the continuum, the semidiscrete, and the fully discrete setting. We prove existence-uniqueness, regularity, and approximation properties for the respective schemes, and set the stage for a controlled numerical computation of delicate nonlinear and nonlocal features of the Hartree dynamics in various physical applications. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{proposition}[theorem]{Proposition} \newtheorem{lemma}[theorem]{Lemma} \newtheorem{definition}[theorem]{Definition} \newtheorem{remark}[theorem]{Remark} \newtheorem{assumption}[theorem]{Assumption} \section{Introduction} In this paper, we study the nonlinear and nonlocal Hartree initial-boundary value problem for the (wave) function $\psi(x,t)$ being defined by \begin{equation}\label{def:cp} \begin{gathered} \mathop{\rm i}\dot\psi =(-\Delta+v+\lambda V\ast |\psi|^2)\,\psi, \quad \text{if }(x,t)\in \overline{\Omega}\times [0,T),\\ \psi=0, \quad\text{if } (x,t)\in \partial\Omega\times [0,T),\\ \psi=\psi_0, \quad\text{if } (x,t)\in \overline{\Omega}\times\{0\}, \end{gathered} \end{equation} where $\Omega$ is some domain in $\mathbb{R}^d$ with boundary $\partial\Omega$, and $T>0$ is the upper limit of the time interval on which we want to study the time evolution of $\psi$ (here, ${\dot \psi}$ is the derivative of $\psi$ with respect to the time variable $t$). Moreover, $v$ stands for an external potential, $\lambda$ denotes the coupling strength, and $V$ is the interaction potential responsible for the nonlinear and nonlocal interaction generated by the convolution term (to be made precise below). The system \eqref{def:cp} has many physical applications, in particular for the case $\Omega=\mathbb{R}^d$. As a first application, we mention the appearance of \eqref{def:cp} within the context of the quantum mechanical description of large systems of nonrelativistic bosons in their so-called mean field limit. For the case of a local nonlinearity, i.e. for $V=\delta$, an important application of equation \eqref{def:cp} lies in the domain of Bose-Einstein condensation for repulsive interatomic forces where it governs the condensate wave function and is called the Gross-Pitaevskii equation. This dynamical equation, and its corresponding energy functional in the stationary case, have been derived rigorously (see for example \cite{Spohn, Erdos} and \cite{Lieb}, respectively; moreover, the nonlocal Hartree equation has been derived for weakly coupled fermions in \cite{Elgart}). In \cite{Aschbacher}, minimizers of this nonlocal Hartree functional have been studied in the attractive case, and symmetry breaking has been established for sufficiently large coupling. A large coupling phase segregation phenomenon has also been rigorously derived for a system of two coupled Hartree equations which are used to describe interacting Bose-Einstein condensates (see \cite{Squassina, Aschbacher08} and references therein). Such coupled systems also appear in the description of crossing sea states of weakly nonlinear dispersive surface water waves in hydrodynamics (with local nonlinearity, see for example \cite{Onorato,Shukla}), of electromagnetic waves in a Kerr medium in nonlinear optics, and in nonlinear plasma physics. Furthermore, we would like to mention that equation \eqref{def:cp} with attractive interaction potential possesses a so-called point particle limit. Consider the situation where the initial condition is composed of several interacting Hartree minimizers sitting in an external potential which varies slowly on the length scale defined by the extension of the minimizers. It turns out that, in a time regime inversely related to this scale, the center of mass of each minimizer follows a trajectory which is governed, up to a small friction term, by Newton's equation of motion for interacting point particles in the slowly varying external potential. Hence, in this limit, the system can be interpreted as the motion of interacting extended particles in a shallow external potential and weakly coupled to a dispersive environment with which mass and energy can be exchanged through the friction term. This allows to describe, and hence to numerically compute, some type of structure formation in Newtonian gravity (see \cite{Froehlich1, Froehlich2, diss}). The main content of the present paper consists in setting up the framework for the numerical analysis on bounded $\Omega$ which will be used in \cite{A09} for the study through numerical computation of such phenomena, like, for example, the dissipation through radiation for a Hartree minimizer oscillating in an external confining potential (see also \cite{diss}). In Section \ref{sec:cp}, we start off with a brief study of the Hartree initial-value boundary problem \eqref{def:cp} in the continuum setting and we discuss its existence-uniqueness and regularity properties. In Section \ref{sec:semidiscrete}, the system \eqref{def:cp} is discretized in space with the help of Galerkin theory. We derive existence-uniqueness and a bound on the $L^2$-approximation error. In the main Section \ref{sec:fd}, we proceed to the full discretization of \eqref{def:cp}, more precisely, we discretize the foregoing semidiscrete problem in time focusing on two time discretization schemes of Crank-Nicholson type. The first is the so-called one-step one-stage Gauss-Legendre Runge-Kutta method which conserves the mass of the discretized wave function under the discrete time evolution. The second one is the so-called Delfour-Fortin-Payre scheme which, besides the mass, also conserves the energy of the system. We prove existence-uniqueness using contraction methods suitable for implementation in \cite{diss, A09}. Moreover, we derive a time quadratic accuracy estimate on the $L^2$-error of these approximation schemes. In the proofs of these assertions, we write down rather explicit expressions for the bounds in order to have some qualitative idea how to achieve a good numerical control of the fully discrete approximations of the Hartree initial-value boundary problem \eqref{def:cp} for the computation of delicate nonlinear and nonlocal features of the various physical scenarios discussed above. \section{The continuum problem} \label{sec:cp} As discussed in the Introduction, we start off by briefly studying the Hartree initial-boundary value problem \eqref{def:cp} in a suitable continuum setting. For this purpose, we make the following assumptions concerning the domain $\Omega$, the external potential $v$, and the interaction potential $V$, a choice which is motivated by the perspective of the fully discrete problem and the numerical analysis dealt with in Section \ref{sec:fd} and the numerical computations in \cite{A09} (some Hartree-dynamical computations have already been performed in \cite{diss}). \begin{assumption} \label{ass:Omega} \rm $\Omega\subset\mathbb{R}^d$ is a bounded domain with smooth boundary $\partial\Omega$. \end{assumption} \begin{assumption} \label{ass:Vv} \rm Assumption \ref{ass:Omega} holds, $v\in C^\infty_0(\Omega,\mathbb{R})$, and $V\in C^\infty_0(\mathbb{R}^d,\mathbb{R})$. \end{assumption} \begin{assumption} \label{ass:even} \rm Assumption \ref{ass:Vv} holds, and $V(-x)=V(x)$ for all $x\in\mathbb{R}^d$. \end{assumption} The Lebesgue and Sobolev spaces used in the following are always defined over the domain $\Omega$ from Assumption \ref{ass:Omega} unless something else is stated explicitly. Thus, we suppress $\Omega$ in the notation of these spaces. Moreover, under Assumption \ref{ass:Vv}, let the Hilbert space $\mathcal{H}$, the linear operator $A$ on $\mathcal{H}$ with domain of definition $D(A)$, and the nonlinear mapping $J$ on $\mathcal{H}$ be given by \begin{gather*} \mathcal{H}:=L^2,\\ D(A):=H^2\cap H^1_0,\\ A:=-\Delta,\\ J[\psi]:=v\psi+f[\psi], \end{gather*} where the nonlinear mapping $f$ is defined by \begin{gather*} f[\psi]:=\lambda g_V[|\psi|^2]\psi,\\ \lambda \in C^\infty_0(\Omega,\mathbb{R}),\\ g_V[\varphi](x) :=\int_\Omega{\rm d}y\,V(x-y) \varphi(y). \end{gather*} \begin{remark} \label{rmk4} \rm The function $\lambda$ stands for some space depending coupling function which can be chosen to be a smooth characteristic function of the domain $\Omega$. Such a choice, on one hand, insures that all derivatives of $f[\psi]$ vanish at the boundary $\partial\Omega$, and, on the other hand, switches the nonlocal interaction off in some neighborhood of $\partial\Omega$ where, in the numerical computation, transparent boundary conditions have to be matched with the outgoing flow of $\psi$ (see \cite{diss, A09}). \end{remark} \begin{remark} \label{rmk5} \rm Under Assumption \ref{ass:Vv}, we have $g_V[\varphi]\in C^\infty_0(\mathbb{R}^d,\mathbb{C})$ for all $\varphi\in L^1$, and $f[\psi]\in L^2$ for all $\psi\in L^2$, see estimate \eqref{est} below. \end{remark} We now make the following definition. \begin{definition} \label{def:cs} \rm Let Assumption \ref{ass:Vv} hold, and let $T\in (0,\infty)\cup \{\infty\}$. We call a differentiable function $\psi:[0,T)\to L^2$ a continuum solution of the Hartree initial-value problem \eqref{def:cp} with initial condition $\psi_0\in D(A)$ if \begin{equation} \label{cp} \begin{gathered} \mathop{\rm i}\dot\psi(t) =A\psi(t)+J[\psi(t)],\quad \forall t\in[0,T),\\ \psi(0) =\psi_0. \end{gathered} \end{equation} If $T<\infty$, the solution is called local, and if $T=\infty$, it is called global. \end{definition} We make use of the following theorem to prove that there exists a unique global solution of the Hartree initial-value problem \eqref{def:cp} in the sense of Definition \ref{def:cs}. In addition, this solution has higher regularity properties in time which are required for the bounds on the constants appearing in the $L^2$-error estimates in the fully discrete setting of Section \ref{sec:fd}. \begin{theorem}[{\cite[p.301]{RS2}}]\label{thm:cs} Let $\mathcal{H}$ be a Hilbert space and $A$ a linear operator on $\mathcal{H}$ with domain of definition $D(A)$, and $A^\ast=A$. \noindent (a) Let $n\in\mathbb{N}$, and let $J$ be a mapping which satisfies the following conditions for all $\psi,\xi\in D(A^k)$, \begin{gather} \label{cond_a} JD(A^k) \subseteq D(A^k),\quad \forall k=1,\dots ,n,\\ \label{cond_b} \|J[\psi]\| \leq C(\|\psi\|) \|\psi\|,\\ \label{cond_c} \|A^k J[\psi]\| \leq C(\|\psi\|,\dots ,\|A^{k-1}\psi\|) \|A^k\psi\|,\quad \forall k=1,\dots ,n,\\ \label{cond_d} \|A^k(J[\psi]-J[\xi])\| \leq C(\|\psi\|,\|\xi\|,\dots ,\|A^k\psi\|,\|A^k\xi\|) \|A^k\psi-A^k\xi\|,\quad \forall k=0,\dots ,n, \end{gather} where each constant $C$ is a monotone increasing and everywhere finite function of all its variables. Then, for each $\psi_0\in D(A^n)$, there exists a unique local continuum solution in the sense of Definition \ref{def:cs} with $\psi(t)\in D(A^n)$ for all $t\in[0,T)$. \noindent (b) Moreover, this solution is global if \begin{equation} \label{boundedness} \|\psi(t)\| \le C,\quad \forall t\in[0,T). \end{equation} \noindent (c) In addition to the conditions in (a), let $J$ satisfy the following conditions for all $1\le k< n$: if $\psi\in C^k([0,\infty),\mathcal{H})$ with $\frac{{\rm d}^l}{{\rm d} t^l}\psi(t) \in D(A^{n-l})$ for all $0\le l\le k$ and $A^{n-l}\frac{{\rm d}^l}{{\rm d} t^l}\psi(t)\in C([0,\infty),\mathcal{H})$ for all $0\le l\le k$, then \begin{gather} \label{cond_e} J[\psi(t)] \text{ is $k$ times differentiable},\\ \label{cond_f} \frac{{\rm d}^k}{{\rm d} t^k}J[\psi(t)] \in D(A^{n-k-1}),\\ \label{cond_g} A^{n-k-1}\frac{{\rm d}^k}{{\rm d} t^k}J[\psi(t)] \in C([0,\infty),\mathcal{H}). \end{gather} If this condition holds, then the local solution from (a) is $n$ times differentiable and $\frac{{\rm d}^k}{{\rm d} t^k}\psi(t) \in D(A^{n-k})$ for all $1\le k\le n$. \end{theorem} \begin{remark} \label{rmk8} \rm With the help of Theorem \ref{thm:cs} (and its proof) the constants in the estimates \eqref{cfd:c4} and \eqref{ifd:c4} below on the $L^2$-error of the fully discrete approximations are finite and can be estimated explicitly. \end{remark} \begin{proof}[Proof] To prove the assertions of Theorem \ref{thm:cs}, we verify its assumptions for the situation specified above. Let us start off by checking condition \eqref{cond_b}. \\ \emph{Condition} \eqref{cond_b}: ${\|J[\psi]\|}_{L^2}\le C({\|\psi\|}_{L^2}) {\|\psi\|}_{L^2}$\\ Using Assumption \ref{ass:Vv} and the estimate \begin{equation} \label{est} {\|\lambda g_V[\varphi\psi]\chi\|}_{L^2} \le {\|\lambda\|}_{L^\infty} {\|V\|}_{L^\infty(\mathbb{R}^d)} {\|\varphi\|}_{L^2}{\|\psi\|}_{L^2}{\|\chi\|}_{L^2}, \end{equation} we immediately get \[ {\|J[\psi]\|}_{L^2} \le ({\|v\|}_{L^\infty} +{\|\lambda\|}_{L^\infty} {\|V\|}_{L^\infty(\mathbb{R}^d)}{\|\psi\|}_{L^2}^2) {\|\psi\|}_{L^2}, \] and the prefactor $C(\|\psi\|)$ from \eqref{cond_b} is a monotone increasing function of ${\|\psi\|}_{L^2}$. Next, let us check condition \eqref{cond_c}. \emph{Condition} \eqref{cond_c}: ${\|\Delta^k J[\psi]\|}_{L^2} \le C({\|\psi\|}_{L^2}){\|\Delta^k\psi\|}_{L^2}$ for all $k\in\mathbb{N}$\\ To show \eqref{cond_c}, we have to control the $L^2$-norm of $\Delta^k(v\psi)$ and of $\Delta^k(\lambda g_V[|\psi|^2]\psi)$. Hence, we write the powers of the Dirichlet-Laplacian as follows, \begin{equation} \label{Deltak1} \begin{aligned} &\Delta^k(\varphi\psi)\\ &=\sum_{\substack{j_1,\dots ,j_k=1,\dots ,d\\ \alpha_1,\beta_1,\dots ,\alpha_k,\beta_k=0,1}} \hspace{-5mm}c_{\alpha_1\beta_1\dots \alpha_k\beta_k}\, (\partial_{j_1}^{\alpha_1}\partial_{j_1}^{\beta_1} \dots \partial_{j_k}^{\alpha_k}\partial_{j_k}^{\beta_k}\varphi) (\partial_{j_1}^{1-\alpha_1}\partial_{j_1}^{1-\beta_1} \dots \partial_{j_k}^{1-\alpha_k}\partial_{j_k}^{1-\beta_k}\psi), \end{aligned} \end{equation} \begin{equation} \label{Deltak2} \begin{aligned} &\Delta^k(g_V[\varphi]\psi)\\ &= \sum_{\substack{j_1,\dots ,j_k=1,\dots ,d\\ \alpha_1,\beta_1,\dots ,\alpha_k,\beta_k=0,1}} \hspace{-5mm}c_{\alpha_1\beta_1\dots \alpha_k\beta_k}\, g_{\partial_{j_1}^{\alpha_1}\partial_{j_1}^{\beta_1} \dots \partial_{j_k}^{\alpha_k}\partial_{j_k}^{\beta_k}V}[\varphi]\, (\partial_{j_1}^{1-\alpha_1}\partial_{j_1}^{1-\beta_1} \dots \partial_{j_k}^{1-\alpha_k}\partial_{j_k}^{1-\beta_k}\psi), \end{aligned} \end{equation} where $c_{\alpha_1\beta_1\dots \alpha_k\beta_k}$ denote some combinatorial constants. Hence, using \eqref{est} and the following Schauder type estimate\footnote{I.e., the elliptic regularity estimate of generalized solutions up to the boundary, see for example \cite[p.383]{Zeidler}.} \[ {\|\psi\|}_{H^{2+m}} \le C {\|\Delta\psi\|}_{H^m},\quad \forall \psi\in H^{2+m}\cap H^1_0,\quad \forall m\in\mathbb{N}_0, \] we get the following bounds on \eqref{Deltak1} and \eqref{Deltak2}, \begin{gather}\label{est1} {\|\Delta^k(\varphi\psi)\|}_{L^2} \le C {\|\varphi\|}_{W^{2k,\infty}} {\|\psi\|}_{H^{2k}} \le C {\|\varphi\|}_{W^{2k,\infty}} {\|\Delta^k\psi\|}_{L^2},\\ \label{est2} {\|\Delta^k(g_V[\varphi]\psi)\|}_{L^2} \le C {\|V\|}_{W^{2k,\infty}(\mathbb{R}^d)} {\|\varphi\|}_{L^1} {\|\Delta^k\psi\|}_{L^2}. \end{gather} Using \eqref{est1}, and \eqref{est2} twice, we arrive at \begin{equation} \label{cond_cc} {\|\Delta^k J[\psi]\|}_{L^2} \le C ({\|v\|}_{W^{2k,\infty}} +{\|\lambda\|}_{W^{2k,\infty}} {\|V\|}_{W^{2k,\infty}(\mathbb{R}^d)} {\|\psi\|}_{L^2}^2) {\|\Delta^k\psi\|}_{L^2}, \end{equation} where the prefactor $C(\|\psi\|, \|A\psi\|,\dots ,\|A^{k-1}\psi\|)$ from \eqref{cond_c} depends on ${\|\psi\|}_{L^2}$ only and is monotone increasing in ${\|\psi\|}_{L^2}$. Next, we check condition \eqref{cond_a}. \emph{Condition} \eqref{cond_a}: $JD(\Delta^k)\subseteq D(\Delta^k)$ for all $k\in\mathbb{N}$\\ Due to \eqref{cond_cc} and since $D(\Delta^k)=\{\psi\in H^2\cap H^1_0\,|\, \Delta\psi\in H^2\cap H^1_0,\dots ,\Delta^{k-1}\psi\in H^2\cap H^1_0\}$, it remains to show that $\Delta^{j}J(\psi)\in H^1_0$ for all $\psi\in D(\Delta^k)$ and for all $j=0,\dots ,k$. But this follows since $v,\lambda\in C^\infty_0$ and from the fact that $C^\infty(\overline{\Omega})$ is dense in $H^m$ with respect to the $H^m$-norm for all $m\in\mathbb{N}_0$. \emph{Condition} \eqref{cond_d}: \[ {\|\Delta^k(J[\psi]-J[\xi])\|}_{L^2} \le C({\|\psi\|}_{L^2},{\|\xi\|}_{L^2}, {\|\Delta^k\psi\|}_{L^2}, {\|\Delta^k\xi\|}_{L^2}) {\|\Delta^k\psi-\Delta^k\xi\|}_{L^2} \] To show \eqref{cond_d}, we write the difference with the help of the decomposition \begin{equation}\label{dc1} \begin{aligned} &g_V[\varphi_1\psi_1]\chi_1-g_V[\varphi_2\psi_2]\chi_2\\ &=g_V[\varphi_1(\psi_1-\psi_2)]\chi_1 +\,g_V[(\varphi_1-\varphi_2)\psi_2]\chi_1 +\,g_V[\varphi_2\psi_2](\chi_1-\chi_2). \end{aligned} \end{equation} Each term on the right-hand side of \eqref{dc1} can then be estimated with the help of \eqref{est2}. Hence, as in \eqref{cond_cc}, we get \begin{equation} \label{cond_dd} \begin{aligned} {\|\Delta^k(J[\psi]-J[\xi])\|}_{L^2} &\leq C ({\|v\|}_{W^{2k,\infty}} +{\|\lambda\|}_{W^{2k,\infty}} {\|V\|}_{W^{2k,\infty}(\mathbb{R}^d)}\\ &\quad\times ({\|\psi\|}_{L^2}^2 +({\|\psi\|}_{L^2}+{\|\xi\|}_{L^2}){\|\Delta^k\psi\|}_{L^2})) {\|\Delta^k\psi-\Delta^k\xi\|}_{L^2}, \end{aligned} \end{equation} where the prefactor $C(\|\psi\|,\|\xi\|,\dots ,\|A^k\psi\|,\|A^k\xi\|)$ from \eqref{cond_d} depends on the lowest and the highest power of $\Delta$ only, and it is monotone increasing in all its variables. \emph{Condition} \eqref{boundedness}: ${\|\psi(t)\|}_{L^2}={\|\psi_0\|}_{L^2}$ for all $t\in [0,T)$\\ This condition is satisfied due to Proposition \ref{prop:conservation} below. \emph{Condition} \eqref{cond_e}:\quad $J[\psi(t)]$ is $k$ times differentiable in $L^2$ with respect to time $t$ Let $n\in\mathbb{N}$ be fixed, and let $k=1$. Then, it is shown in \cite[p.299]{RS2} that part (a) and conditions \eqref{cond_e}, \eqref{cond_f}, and \eqref{cond_g} for $k=1$ imply that $\psi\in C^2([0,\infty),L^2)$ with $\frac{{\rm d}^2}{{\rm d} t^2}\psi(t) \in D(A^{n-2})$ and $A^{n-2}\frac{{\rm d}^2}{{\rm d} t^2}\psi(t)\in C([0,\infty),L^2)$. Then, using conditions \eqref{cond_e}, \eqref{cond_f}, and \eqref{cond_g} for subsequent $1\le k0$ whose closure is the union of the $(n-1)^2$ congruent closed subsquares generated by dividing each side of $\Omega$ equidistantly into $n-1$ intervals. Let us denote by $N_h=(n-2)^2$ the total number of interior vertices of this lattice and by $h=D/(n-1)$ the lattice spacing.\footnote{\label{foot:bijection}As bijection from the one-dimensional to the two-dimensional lattice numbering, we may use the mapping $\tau:\{0,\dots ,m-1\}^{\times 2} \to \{0,\dots ,m^2-1\}$ with $j=\tau(m_1,m_2):=m_1+m_2m$.} Moreover, let us choose the Galerkin space $S_h$ to be spanned by the bilinear Lagrange rectangle finite elements $\varphi_j\in C(\overline{\Omega})$ whose reference basis function $\varphi_0:\overline{\Omega}\to [0,\infty)$ is defined on its support $[0,2h]^{\times 2}$ by \begin{equation} \label{def:blfe} \varphi_0(x,y):=\frac{1}{h^2} \begin{cases} xy, & \text{if } (x,y)\in [0,h]^{\times 2},\\ (2h-x)y, & \text{if } (x,y)\in [h,2h]\times [0,h],\\ (2h-x)(2h-y), & \text{if } (x,y)\in [h,2h]^{\times 2},\\ x(2h-y), & \text{if } (x,y)\in [0,h]\times [h,2h], \end{cases} \end{equation} see Figure \ref{fig:refel}. The functions $\varphi_j$ are then defined to be of the form \eqref{def:blfe} having their support translated by $(m_1 h,m_2 h)$ with $m_1,m_2=0,\dots ,m-3$. Hence, with this choice, we have $S_h\subset C(\overline{\Omega})\cap H_0^1(\Omega)$ and $\dim S_h=N_h$. \end{remark} \begin{figure}[ht] \begin{center} \includegraphics[width=4cm,height=4.5cm]{fig1} \caption{$\varphi_0(x,y)$ on its support $[0,2h]^{\times 2}$ with maximum at vertex $(h,h)$.} \label{fig:refel} \end{center} \end{figure} Motivated by the weak formulation \eqref{def:wcp}, we make the following definition. \begin{definition} \label{def:sdp} \rm Let Assumptions \ref{ass:Vv} and \ref{ass:Sh} hold. We call $\psi_h: [0,T)\to S_h$ with $\psi_h,\dot\psi_h\in L^2(0,T;S_h)$ a semidiscrete solution of the Hartree initial- boundary value problem \eqref{def:cp} with initial condition $\psi_{0h}\in S_h$ if \begin{equation} \label{sdp} \begin{gathered} \mathop{\rm i}\, \frac{{\rm d}}{{\rm d}t}\,{(\varphi,\psi_h)}_{L^2} ={(\nabla\varphi,\nabla\psi_h)}_{L^2} +{(\varphi,v\psi_h)}_{L^2} +{(\varphi,f[\psi_h])}_{L^2},\quad \forall\varphi\in S_h,\;t\in [0,T),\\ \psi_h(0)=\psi_{0h}. \end{gathered} \end{equation} \end{definition} \begin{remark} \label{rem:Gelfand} \rm In general, the weak problem \eqref{def:wcp} is set up using the Gelfand evolution triple $H_0^1\subset L^2\subset {(H_0^1)}^\ast=H^{-1}$. One then looks for weak solutions $\psi\in W_2^1(0,T; H_0^1, L^2)\subset C([0,T),L^2)$ motivating Definition \ref{def:sdp}. \end{remark} We assume the Galerkin subspace $S_h$ from Assumption \ref{ass:Sh} to satisfy the following additional approximation and inverse inequalities. \begin{assumption} \label{ass:AccuracyOrder} \rm Let Assumption \ref{ass:Sh} hold. Then, there exists a constant $C_A>0$ such that \[ \inf_{\varphi\in S_h}\big({\|\psi-\varphi\|}_{L^2} +h{\|\psi-\varphi\|}_{H^1}\big) \le C_A h^2 {\|\psi\|}_{H^2}, \quad\forall \psi\in H^2\cap H^1_0. \] \end{assumption} \begin{remark} \label{rmk16} \rm For an order of accuracy $r\ge 2$ of the family $\{S_h\}_{h\in (0,1)}$, the usual assumption replaces the right-hand side by $C_A h^s {\|\psi\|}_{H^s}$ for $1\le s\le r$ and is asked to hold for all $\psi\in H^s\cap H^1_0$. For simplicity, we stick to Assumption \ref{ass:AccuracyOrder}. \end{remark} \begin{assumption} \label{ass:inverse} \rm Let Assumption \ref{ass:Sh} hold. Then, there exists a constant $C_B>0$ such that \[ {\|\varphi\|}_{H^1}\le C_B h^{-1} {\|\varphi\|}_{L^2},\quad\forall \varphi\in S_h. \] \end{assumption} \begin{remark} \label{rmk18} \rm For the two-dimensional bilinear Lagrange finite element setting of Remark \ref{rem:Lagrange}, both Assumption \ref{ass:AccuracyOrder} and Assumption \ref{ass:inverse} hold (see for example \cite[p.109,111]{Brenner}). \end{remark} Furthermore, we make an assumption on the approximation quality of the initial condition $\psi_{0h}\in S_h$ of the semidiscrete problem \eqref{sdp} compared to the initial condition $\psi_0\in H^2\cap H^1_0$ of the continuum problem \eqref{cp}. \begin{assumption} \label{ass:initial} \rm Let Assumption \ref{ass:Sh} hold. Then, there exists a constant $C_0>0$ such that \begin{equation} \label{def:initial} {\|\psi_0-\psi_{0h}\|}_{L^2}\le C_0 h^2. \end{equation} \end{assumption} The semidiscrete scheme has the following conservation properties. \begin{proposition} \label{sdp:ExUnique} Let Assumptions \ref{ass:even} and \ref{ass:Sh} hold, and let $\psi_h$ be a semidiscrete solution of the Hartree initial-boundary value problem \eqref{def:cp} in the sense of Definition \ref{def:sdp}. Then, the mass and energy of $\psi_h$ are conserved under the time evolution, \begin{equation} \label{sdp:conservation} \begin{gathered} \mathcal{M}[\psi_h(t)]=\mathcal{M}[\psi_{0h}],\quad \forall t\in[0,T),\\ \mathcal{H}[\psi_h(t)]=\mathcal{H}[\psi_{0h}],\quad \forall t\in[0,T). \end{gathered} \end{equation} \end{proposition} \begin{proof} If we plug $\varphi=\psi_h(t)$ into \eqref{sdp} and take the imaginary part of the resulting equation, we get the conservation of the mass. If we plug $\varphi=\dot\psi_h(t)$ into \eqref{sdp} and take the real part of the resulting equation, we get the conservation of the energy using Assumption \ref{ass:even} (in the sense of Remark \ref{rem:Gelfand}). \end{proof} Existence-uniqueness is addressed in the following theorem. \begin{theorem} \label{sdp:ExUniqueCons} Let Assumptions \ref{ass:even} and \ref{ass:Sh} hold. Then, there exists a unique global semidiscrete solution $\psi_h$ of the Hartree initial-boundary value problem \eqref{def:cp} in the sense of Definition \ref{def:sdp}. \end{theorem} \begin{proof} Let $\{\varphi_j\}_{j=1}^{N_h}$ be a basis of the Galerkin space $S_h$, and let us write \begin{equation} \label{expansion1} \psi_h(t) =\sum_{j=1}^{N_h}z_j(t)\,\varphi_j. \end{equation} Plugging \eqref{expansion1} into the semidiscrete system \eqref{sdp}, for $z(t):=(z_1(t),\dots ,z_{N_h}(t))\in\mathbb{C}^{N_h}$, we obtain \begin{equation}\label{sdp2} \begin{gathered} \mathop{\rm i}\dot z(t) =A^{-1}(B+Y)z(t) +A^{-1}H[z(t)]z(t), \quad\forall t\in[0,T),\\ z(0)=z_0, \end{gathered} \end{equation} where $\psi_{0h}=\sum_{j=1}^{N_h}(z_0)_j\,\varphi_j$ and the matrices $A,B\in\mathbb{C}^{N_h\times N_h}$ are the positive definite mass and stiffness matrices, respectively, \[ A_{ij}:={(\varphi_i,\varphi_j)}_{L^2},\quad B_{ij}:={(\nabla\varphi_i,\nabla\varphi_j)}_{L^2}. \] Moreover, $Y\in\mathbb{C}^{N_h\times N_h}$ is the external potential matrix, \[ Y_{ij}:={(\varphi_i,v\varphi_j)}_{L^2}, \] and the matrix-valued function $H:\mathbb{C}^{N_h}\to\mathbb{C}^{N_h\times N_h}$ is defined by \[ {H[z]}_{ij} :=\sum_{k,l=1}^{N_h}\bar z_k z_l\, {(\varphi_i,\lambda g[\bar\varphi_k\varphi_l]\varphi_j)}_{L^2}. \] Since the function $\mathbb{C}^{N_h}\ni z\mapsto A^{-1}(B+Y)z+A^{-1}H[z]z\in\mathbb{C}^{N_h}$ is locally Lipschitz continuous analogously to the continuum case, the Picard-Lindel\"of theory for ordinary differential equations implies local existence and uniqueness of the initial-value problem \eqref{sdp2}. Moreover, this local solution is a global solution if it remains restricted to a compact subset of $\mathbb{C}^{N_h}$. But this is the case due to the mass conservation from \eqref{sdp:conservation}. \end{proof} We next turn to the $L^2$-error estimate of the semidiscretization. For that purpose, we introduce the Ritz projection (also called elliptic projection). \begin{definition} \label{def:Ritz} \rm Let Assumption \ref{ass:Sh} hold. Then, the Ritz projection $R_h: H^1_0\to S_h$ is defined to be the orthogonal projection from $H^1_0$ onto $S_h$ with respect to the Dirichlet scalar product ${(\nabla\cdot,\nabla\cdot)}_{L^2}$ on $H^1_0$; i.e., \[ \label{Ritz} {(\nabla\varphi,\nabla R_h\psi)}_{L^2} ={(\nabla\varphi,\nabla\psi)}_{L^2}, \quad\forall \varphi\in S_h. \] \end{definition} The Ritz projection satisfies the following error estimate. \begin{lemma}[{\cite[p.8]{VT97}}] Let Assumptions \ref{ass:Sh} and \ref{ass:AccuracyOrder} hold. Then, there exists a constant $C_R>0$ such that \begin{equation} \label{RitzEstimate} {\|(1-R_h)\psi\|}_{L^2}+h{\|(1-R_h)\psi\|}_{H^1} \le C_R h^2 {\|\psi\|}_{H^2}, \quad\forall \psi\in H^2\cap H^1_0. \end{equation} \end{lemma} The next theorem is the main assertion of this section. \begin{theorem} \label{thm:ErrorSD} Let Assumptions \ref{ass:Vv}, \ref{ass:Sh}, \ref{ass:AccuracyOrder}, and \ref{ass:initial} hold, and let $\psi$ be the solution of the continuum problem from Theorem \ref{thm:cs} and $\psi_h$ the solution of the semidiscrete problem from Theorem \ref{sdp:ExUniqueCons}. Then, for any $00$ such that \[ \max_{t\in[0,T]}{\|\psi(t)-\psi_h(t)\|}_{L^2} \le C_E h^2. \] \end{theorem} \begin{proof} We decompose the difference of $\psi$ and $\psi_h$ as \[ \psi(t)-\psi_h(t)=\rho(t)+\theta(t), \] where $\rho(t)$ and $\theta(t)$ are defined with the help of the Ritz projection $R_h$ from \eqref{Ritz} by \[ \rho(t):=(1-R_h)\psi(t),\quad \theta(t):=R_h\psi(t)-\psi_h(t). \] Making use of the schemes \eqref{def:wcp}, \eqref{sdp}, and \eqref{Ritz}, we can write \begin{equation} \label{ThetaAndRho} \begin{aligned} &\mathop{\rm i}{(\varphi,{\dot \theta}(t))}_{L^2} -{(\nabla\varphi,\nabla\theta(t))}_{L^2}\\ &= -\mathop{\rm i}{(\varphi,{\dot \rho}(t))}_{L^2} +{(\varphi,v(\psi(t)-\psi_h(t)))}_{L^2} + {(\varphi,f[\psi(t)]-f[\psi_h(t)])}_{L^2}. \end{aligned} \end{equation} Plugging $\varphi=\theta(t)\in S_h$ into \eqref{ThetaAndRho} and taking the imaginary part of the resulting equation, we get the differential inequality \begin{equation} \label{DiffInequality1} \begin{aligned} &\frac12\,\frac{{\rm d}}{{\rm d}t}\,{\|\theta(t)\|}_{L^2}^2\\ &\le \left({\|{\dot \rho}(t)\|}_{L^2} +{\|v(\psi(t)-\psi_h(t))\|}_{L^2} +{\|f[\psi(t)]-f[\psi_h(t)]\|}_{L^2}\right) {\|\theta(t)\|}_{L^2}\\ &\le \left({\|{\dot \rho}(t)\|}_{L^2} +c_1\left({\|\rho(t)\|}_{L^2}+{\|\theta(t)\|}_{L^2}\right) \right) {\|\theta(t)\|}_{L^2}, \end{aligned} \end{equation} where we used the conservation laws \eqref{cp:conservation} and \eqref{sdp:conservation}, and \eqref{cond_dd} to define the constant \[ c_1 :={\|v\|}_{L^\infty}+2{\|\lambda\|}_{L^\infty} {\|V\|}_{L^\infty(\mathbb{R}^d)}\left(\mathcal{M}[\psi_0]+\mathcal{M}[\psi_{0h}]\right). \] Using $\epsilon>0$ to regularize the time derivative of ${\|\theta\|}_{L^2}$ at $\theta=0$ by rewriting the left-hand side of \eqref{DiffInequality1} as $\frac12\frac{{\rm d}}{{\rm d}t}{\|\theta\|}_{L^2}^2 =\frac12\frac{{\rm d}}{{\rm d}t}({\|\theta\|}_{L^2}^2+\epsilon^2)$, we get \begin{equation} \label{DiffInequality2} \frac{{\rm d}}{{\rm d}t}\left({\|\theta\|}_{L^2}^2+\epsilon^2\right)^{1/2} \le {\|{\dot \rho}(t)\|}_{L^2} +c_1\left({\|\rho(t)\|}_{L^2}+{\|\theta(t)\|}_{L^2}\right), \end{equation} where we used ${\|\theta(t)\|}_{L^2}\le \big({\|\theta\|}_{L^2}^2+\epsilon^2\big)^{1/2}$. Integrating \eqref{DiffInequality2} from $0$ to $t$, letting $\epsilon\to 0$, and applying Gr\"onwall's lemma to the resulting inequality, we find \begin{equation}\label{Groenwall} \begin{aligned} {\|\theta(t)\|}_{L^2} &\le {\|\theta(0)\|}_{L^2} +\int_0^t {\rm d}s\,\big({\|{\dot \rho}(s)\|}_{L^2} +c_1{\|\rho(s)\|}_{L^2}\big)\\ &\quad +\,c_1\int_0^t {\rm d}s\,\Big( {\|\theta(0)\|}_{L^2} +\int_0^s {\rm d}u\,\big({\|{\dot \rho}(u)\|}_{L^2} +c_1{\|\rho(u)\|}_{L^2}\big)\Big)\, {\rm e}^{c_1(t-s)}. \end{aligned} \end{equation} To extract the factor $h^2$, we apply \eqref{RitzEstimate} and Assumption \ref{ass:initial} to obtain \begin{gather} \label{rhoEstimate} {\|\rho(t)\|}_{L^2} \le C_R{\|\psi(t)\|}_{H^2} h^2,\\ \label{rhotEstimate} {\|{\dot \rho}(t)\|}_{L^2} \le C_R{\|{\dot \psi}(t)\|}_{H^2} h^2,\\ \label{theta0Estimate} {\|\theta(0)\|}_{L^2} \le \big(C_0+C_R{\|\psi_0\|}_{H^2}\big) h^2. \end{gather} Plugging \eqref{rhoEstimate}, \eqref{rhotEstimate}, and \eqref{theta0Estimate} into \eqref{Groenwall}, we finally arrive at \[ {\|\psi(t)-\psi_h(t)\|}_{L^2} \le {\|\theta(t)\|}_{L^2}+{\|\rho(t)\|}_{L^2} \le c_2(t) h^2, \] where the time dependent prefactor is defined by \begin{align*} c_2(t) &:=\big(C_0+C_R{\|\psi_0\|}_{H^2}\big) {\rm e}^{c_1 t} +C_R\int_0^t {\rm d}s \,\big( {\|{\dot \psi}(s)\|}_{H^2}+ c_1{\|\psi(s)\|}_{H^2}\big)\\ &\quad +c_1 C_R \int_0^t {\rm d}s\int_0^s {\rm d}u\, \big({\|{\dot \psi}(u)\|}_{H^2}+c_1{\|\psi(u)\|}_{H^2}\big) \,{\rm e}^{c_1(t-s)} +C_R {\|\psi(t)\|}_{H^2}. \end{align*} Setting $C_E:=\max_{t\in[0,T]}c_2(t)$ brings the proof of Theorem \ref{thm:ErrorSD} to an end. \end{proof} \begin{remark} \label{rmk25} \rm For the local case $f[\psi]=|\psi|^2\psi$, one replaces the original locally Lipschitz nonlinearity $f$ by a globally Lipschitz continuous nonlinearity which coincides with $f$ in a given neighborhood of the solution $\psi$ of the continuum problem. One then first shows that the semidiscrete solution of the modified problem satisfies the desired $L^2$-error bound, and, second, that for $h$ sufficiently small, the modified solution lies in the given neighborhood of $\psi$. But for such $h$, the solution of the modified problem coincides with the solution of the original problem, and, hence, the solution of the original problem satisfies the desired $L^2$-error bound, too (see for example \cite{ADK91}). \end{remark} \section{The fully discrete approximation} \label{sec:fd} In this section, we discretize the semidiscrete problem \eqref{sdp} in time. To this end, let us denote by $N\in\mathbb{N}$ the desired fineness of the time discretization with time discretization scale $\tau$ and its multiples $t_n$ for all $n=0,1,2,\dots ,N$, \begin{equation} \label{def:tau} \tau:=\frac{T}{N},\quad t_n:=n\tau. \end{equation} As mentioned in the Introduction, we will use two different time discretization schemes of Crank-Nicholson type to approximate the semidiscrete solution $\psi_h$ of Theorem \ref{sdp:ExUniqueCons} at time $t_n$ by $\Psi_n\in\Psi$, where \begin{equation} \label{def:Psi} \Psi:=(\Psi_0,\Psi_1,\dots ,\Psi_N)\in S_h^{\times(N{+}1)}. \end{equation} These two schemes differ in the way of approximating the nonlinear term $g_V[|\psi|^2]\psi$ as follows. Let $\mathcal{N}:=\{1,2,\dots ,N\}$ and $\mathcal{N}_0:=\mathcal{N}\cup\{0\}$, and set \begin{equation} \label{def:average} \Psi_{n-1/2}:=\frac{1}{2}\left(\Psi_{n}+\Psi_{n-1}\right), \quad\forall n\in\mathcal{N}. \end{equation} The first scheme implements the one-step one-stage Gauss-Legendre Runge-Kutta method in which the nonlinear term is discretized by \begin{equation} \label{gCoherent} g_V[|\Psi_{n-1/2}|^2]\Psi_{n-1/2}. \end{equation} In this method, the mass $\mathcal{M}[\Psi_n]$ is conserved under the discrete time evolution. The second scheme, introduced in \cite{DFP81} and applied in \cite{ADK91}, discretizes the nonlinear term by \begin{equation} \label{gIncoherent} g_V[\tfrac12(|\Psi_{n}|^2+|\Psi_{n-1}|^2)]\Psi_{n-1/2}. \end{equation} This method, in addition to the mass, also conserves the energy $\mathcal{H}[\Psi_n]$ of the system. In the following, for convenience, we will call the first scheme \emph{coherent} and the second one \emph{incoherent}. \subsection{Coherent scheme} \label{sec:cfd} To define what we mean by a coherent solution of the Hartree initial-boundary value problem \eqref{def:cp}, we set \[ {\dot \Psi}_n :=\frac{1}{\tau}\left(\Psi_{n}-\Psi_{n-1}\right), \quad\forall n\in\mathcal{N}. \] \begin{definition} \label{def:cfd} \rm Let Assumption \ref{ass:Vv} and \ref{ass:Sh} hold. We call $\Psi\in S_h^{\times(N{+}1)}$ a coherent fully discrete solution of the Hartree initial-boundary value problem \eqref{def:cp} with initial condition $\psi_{0h}\in S_h$ if \begin{equation} \label{cfd} \begin{gathered} \mathop{\rm i}{(\varphi,{\dot \Psi}_n)}_{L^2} ={(\nabla\varphi,\nabla\Psi_{n-1/2})}_{L^2} +\,{(\varphi,v\Psi_{n-1/2})}_{L^2} +{(\varphi,f[\Psi_{n-1/2}])}_{L^2},\\ \forall\varphi\in S_h,\, \forall n\in\mathcal{N},\\ \Psi_0=\psi_{0h}. \end{gathered} \end{equation} \end{definition} The coherent solution has the following conservation property. \begin{proposition} \label{prop:cfdMassConservation} Let $\Psi\in S_h^{\times(N{+}1)}$ be a coherent fully discrete solution of the Hartree initial-boundary value problem \eqref{def:cp} in the sense of Definition \ref{def:cfd}. Then, the mass of $\Psi$ is conserved under the discrete time evolution, \begin{equation} \label{cfdMassConservation} \mathcal{M}[\Psi_n]=\mathcal{M}[\psi_{0h}], \quad\forall n\in\mathcal{N}_0. \end{equation} \end{proposition} \begin{proof} If we plug $\varphi=\Psi_{n-1/2}$ into \eqref{cfd} and take the imaginary part of the resulting equation, we get \[ 0=\mathop{\rm Im}\mathop{\rm i}{(\Psi_{n-1/2},{\dot \Psi}_n)}_{L^2} =\frac{1}{2\tau}\left(\mathcal{M}[\Psi_{n}]-\mathcal{M}[\Psi_{n-1}]\right). \] \end{proof} \begin{remark} \label{rmk28} \rm The energy $\mathcal{H}[\Psi_n]$ of the coherent solution \eqref{cfd} is not conserved under the discrete time evolution (see \cite{ADK91} and references therein, in particular \cite{Sanz-Serna} and \cite{Tourigny} for the local case with $d=1$). \end{remark} The question of existence and uniqueness of a coherent solution is addressed in the following theorem. \begin{theorem} \label{prop:ExistenceCFD} Let Assumptions \ref{ass:Vv}, \ref{ass:Sh}, and \ref{ass:inverse} hold, and let the time discretization scale $\tau$ be sufficiently small. Then, there exists a unique coherent fully discrete solution of the Hartree initial-boundary value problem \eqref{def:cp} in the sense of Definition \ref{def:cfd}. \end{theorem} \begin{proof} Let $\phi\in S_h$ be given, and define the mapping $F_\phi: S_h\to S_h$ by \begin{equation} \label{def:Fn} {(\varphi,F_\phi[\psi])}_{L^2} :={(\varphi,\phi)}_{L^2} -\frac{\mathop{\rm i}\tau}{2}\, \big({(\nabla\varphi,\nabla\psi)}_{L^2} +{(\varphi,v\psi)}_{L^2} +{(\varphi,f[\psi])}_{L^2}\big),\quad \forall \varphi\in S_h. \end{equation} For some $n\in\mathcal{N}$, let the $n$-th component $\Psi_{n-1}$ of $\Psi\in S_h^{\times(N{+}1)}$ from \eqref{def:Psi} be given. Adding $2\mathop{\rm i}\,(\varphi,\Psi_{n-1})/\tau$ on both sides of \eqref{cfd}, we can rewrite \eqref{cfd} with the help of \eqref{def:Fn} in the form of a fixed point equation for $\Psi_{n-1/2}$, \begin{equation} \label{cfdFixedPoint} \Psi_{n-1/2}=F_{\Psi_{n-1}}[\Psi_{n-1/2}], \end{equation} from which we retrieve the unknown component $\Psi_n$ by \eqref{def:average}. In order to construct the unique solution of \eqref{cfdFixedPoint}, we make use of Banach's fixed point theorem on the compact ball $\mathcal{B}_{n-1}:=\{\psi\in S_h\,|\, {\|\psi\|}_{L^2}\le \mathcal{M}[\Psi_{n-1}]^{1/2}+1\}$ in $S_h$. Using Assumption \ref{ass:inverse} and \eqref{cond_dd}, we get, for $\psi,\xi\in S_h$, \begin{equation} \label{banach1} \begin{aligned} &|{(\varphi,F_{\Psi_{n-1}}[\psi]-F_{\Psi_{n-1}}[\xi])}_{L^2}|\\ &\le\frac\tau 2 \Big(C_B^2 h^{-2}+{\|v\|}_{L^\infty} +2{\|\lambda\|}_{L^\infty}{\|V\|}_{L^\infty(\mathbb{R}^d)} \big({\|\psi\|}_{L^2}^2+{\|\xi\|}_{L^2}^2\big)\Big) {\|\psi-\xi\|}_{L^2} {\|\varphi\|}_{L^2}. \end{aligned} \end{equation} Plugging $\varphi=F_{\Psi_{n-1}}[\psi]-F_{\Psi_{n-1}}[\xi]$ into \eqref{banach1} and picking $\psi$ and $\xi$ from $\mathcal{B}_{n-1}$, we find \begin{equation} \label{banach2} {\|F_{\Psi_{n-1}}[\psi]-F_{\Psi_{n-1}}[\xi]\|}_{L^2} \le \frac{\alpha_{n-1}}{2}\,\tau {\|\psi-\xi\|}_{L^2}, \end{equation} where the constant $\alpha_{n-1}$ is defined, for all $n\in\mathcal{N}$, by \begin{equation} \label{taun} \alpha_{n-1}:=C_B^2 h^{-2}+{\|v\|}_{L^\infty} +4{\|\lambda\|}_{L^\infty}{\|V\|}_{L^\infty(\mathbb{R}^d)}(\mathcal{M}[\Psi_{n-1}]^{1/2}+1)^2. \end{equation} Let now $\alpha_{n-1}(\mathcal{M}[\Psi_{n-1}]^{1/2}+1)\tau\le 1$. Then, it follows from \eqref{banach2} and \eqref{taun} that $F_{\Psi_{n-1}}$ maps $\mathcal{B}_{n-1}$ into $\mathcal{B}_{n-1}$ (set $\xi=0$ in \eqref{banach2}) and that $F_{\Psi_{n-1}}$ is a strict contraction on $\mathcal{B}_{n-1}$. Therefore, for such $\tau$, Banach's fixed point theorem implies the existence of a unique solution $\Psi_{n-1/2}\in\mathcal{B}_{n-1}$ of the fixed point equation \eqref{cfdFixedPoint}. Moreover, due to the mass conservation \eqref{cfdMassConservation}, there exists no solution $\Psi_{n-1/2}$ of \eqref{cfdFixedPoint} with $\Psi_{n-1/2}\in S_h\setminus \mathcal{B}_{n-1}$. Hence, the component $\Psi_n$ of the coherent solution exists and is unique for such $\tau$. Starting at $\Psi_0=\psi_{0h}$ and proceeding iteratively, we get all $n+1$ components of the coherent solution $\Psi\in S_h^{\times(N{+}1)}$. Moreover, again due to \eqref{cfdMassConservation}, we get a uniform bound on the size of the time discretization scale $\tau$, e.g. \[ \alpha_0(\mathcal{M}[\psi_{0h}]^{1/2}+1)\tau\le 1. \] \end{proof} \begin{remark} \label{rm30}\rm Since $\alpha_0\ge C_B^2h^{-2}$, we have that $\tau\le C_B^{-2}h^2$, where $C_B$ stems from Assumption \ref{ass:inverse}. \end{remark} We next turn to the first of the two main assertions of the present paper which is the time quadratic accuracy estimate on the $L^2$-error of the coherent solution. \begin{theorem} \label{thm:cfderror} Let Assumptions \ref{ass:Vv}, \ref{ass:Sh}, \ref{ass:AccuracyOrder}, and \ref{ass:initial} hold, and let $\Psi\in S_h^{\times(N{+}1)}$ be the coherent solution from Theorem \ref{prop:ExistenceCFD}. Then, there exists a constant $C_K>0$ such that \[ \max_{n\in\mathcal{N}_0}{\|\psi(t_n)-\Psi_n\|}_{L^2} \le C_K (\tau^2+h^2). \] \end{theorem} \begin{remark} \label{rmk32} \rm The constant $C_K$ depends on higher Sobolev norms of the continuum solution $\psi$. These norms exist due to the regularity assertion in Theorem \ref{thm:cs}(c). \end{remark} \begin{proof} Let $n\in\mathcal{N}$ be fixed and define $\psi_n:=\psi(t_n)$ with $t_n$ from \eqref{def:tau}. As in the proof of Theorem \ref{thm:ErrorSD}, we decompose the difference to be estimated as \begin{equation} \label{decomposition} \psi_n-\Psi_n=\rho_n+\theta_n, \end{equation} where $\rho_n$ and $\theta_n$ are again defined with the help of the Ritz projection from \eqref{Ritz} by \begin{equation} \label{def:RhoTheta} \rho_n:=(1-R_h)\psi_n,\quad \theta_n:=R_h\psi_n-\Psi_n. \end{equation} Using Taylor's theorem in order to expand $\psi_n$ around $t=0$ up to zeroth order in $t_n$ and the estimate on the Ritz projection \eqref{RitzEstimate}, we immediately get \begin{equation} \label{RhoEstimate} {\|\rho_n\|}_{L^2} \le C_Rh^2\Big({\|\psi_0\|}_{H^2} +\int_0^{t_n} {\rm d}t\,{\|{\dot \psi}(t)\|}_{H^2}\Big). \end{equation} To estimate $\theta_n$, we want to extract suitable small differences from the expression \begin{equation} \label{def:Ln} L_{n,\varphi}:= \frac{\mathop{\rm i}}{\tau}{(\varphi,\theta_{n}-\theta_{n-1})}_{L^2} -\frac12 {(\nabla\varphi,\nabla (\theta_{n}+\theta_{n-1}))}_{L^2} -\frac12 {(\varphi,v (\theta_{n}+\theta_{n-1}))}_{L^2} \end{equation} which contains all the linear terms in \eqref{cfd} moved to the left-hand side with $\Psi_n$ replaced by $\theta_n$. For this purpose, we first plug the definition of $\theta_n$ into \eqref{def:Ln}, and then use the definition of the Ritz projection \eqref{Ritz} and the scheme \eqref{cfd} to get \begin{equation} \label{Ln1} \begin{aligned} L_{n,\varphi}&= \frac{\mathop{\rm i}}{\tau}\,{(\varphi,R_h(\psi_{n}-\psi_{n-1}))}_{L^2} -\frac12 {(\nabla\varphi,\nabla (\psi_{n}+\psi_{n-1}))}_{L^2}\\ &\quad -\frac12 {(\varphi,v R_h(\psi_{n}+\psi_{n-1}))}_{L^2} -{(\varphi,f[\Psi_{n-1/2}])}_{L^2}. \end{aligned} \end{equation} Rewriting the first term on the right-hand side of \eqref{Ln1} with the help of the continuum solution satisfying the weak formulation \eqref{def:wcp}, we have \begin{equation} \label{Rh} \begin{aligned} &\frac{\mathop{\rm i}}{\tau}\,{(\varphi,R_h(\psi_{n}-\psi_{n-1}))}_{L^2}\\ &= \frac{\mathop{\rm i}}{\tau}\,{(\varphi,(R_h-1)(\psi_{n}-\psi_{n-1}))}_{L^2} +\mathop{\rm i}\,{(\varphi,\tfrac{1}{\tau}(\psi_{n}-\psi_{n-1})-{\dot \psi}_{n-1/2})}_{L^2} \\ &\quad +\,{(\nabla\varphi,\nabla\psi_{n-1/2})}_{L^2} +{(\varphi,v\psi_{n-1/2})}_{L^2} +{(\varphi,f[\psi_{n-1/2}])}_{L^2}, \end{aligned} \end{equation} where we used $\psi_{n-1/2}:=\psi(t_n-\tau/2)$ and ${\dot \psi}_{n-1/2}:={\dot \psi}(t_n-\tau/2)$. Plugging \eqref{Rh} into \eqref{Ln1}, we can express $L_{n,\varphi}$ in the form \begin{equation} \label{def:omegas} L_{n,\varphi}=\sum_{j=1}^6\,{(\varphi,\omega_n^{(j)})}_{L^2}, \end{equation} where the functions $\omega_n^{(j)}$ with $j=1,\dots ,6$ are defined by \begin{gather*} \omega_n^{(1)} :=\frac{\mathop{\rm i}}{\tau}(R_h-1)(\psi_{n}-\psi_{n-1}),\quad \omega_n^{(2)} :=\mathop{\rm i}\big(\tfrac{1}{\tau}(\psi_{n}-\psi_{n-1}) -{\dot \psi}_{n-1/2}\big),\\ \omega_n^{(3)}:=\Delta\big(\tfrac12\left(\psi_{n}+\psi_{n-1}\big) -\psi_{n-1/2}\right),\quad \omega_n^{(4)} :=v\left(\psi_{n-1/2}-\tfrac12\left(\psi_{n} +\psi_{n-1}\right)\right),\\ \omega_n^{(5)}:=\frac12 \,v\left(1-R_h\right)\left(\psi_{n}+\psi_{n-1}\right), \quad \omega_n^{(6)}:= f[\psi_{n-1/2}]-f[\Psi_{n-1/2}], \end{gather*} and, for $\omega_n^{(3)}$, we used again ${(\nabla\varphi,\nabla\psi)}_{L^2}={(\varphi,-\Delta\psi)}_{L^2}$ for all $\psi\in H^2\cap H^1_0$. Plugging $\varphi=(\theta_{n}+\theta_{n-1})/2$ into \eqref{def:Ln} and \eqref{def:omegas}, and taking the imaginary part of the resulting equation, we get (if $\theta_{n}=0$, we are left with \eqref{RhoEstimate}) \begin{equation} \label{thetaEstimate1} {\|\theta_{n}\|}_{L^2} \le {\|\theta_{n-1}\|}_{L^2} +\tau\sum_{j=1}^6 {\|\omega_n^{(j)}\|}_{L^2}. \end{equation} Let us next estimate the terms ${\|\omega_n^{(j)}\|}_{L^2}$ for all $j=1,\dots ,6$. For $\omega_n^{(1)}$, we expand $\psi_{n}$ around $t=t_{n-1}$ up to zeroth order in $\tau$ and use \eqref{RitzEstimate} such that \begin{equation} \label{omega1Estimate} {\|\omega_n^{(1)}\|}_{L^2} \le {C_Rh^2}\tau^{-1}\int_{t_{n-1}}^{t_{n}} {\rm d}t\, {\|{\dot \psi}(t)\|}_{H^2}. \end{equation} For $\omega_n^{(2)}$, we expand $\psi_{n-1}$ and $\psi_{n}$ around $t=t_n-\tau/2$ up to second order in $\tau/2$, \begin{equation} \label{omega2Estimate} \begin{aligned} {\|\omega_n^{(2)}\|}_{L^2} &\le \frac{1}{2\tau}\Big(\int_{t_{n-1}}^{t_n-\tau/2} {\rm d}t\, (t_{n-1}-t)^2{\|{\dddot \psi}(t)\|}_{L^2} +\int_{t_n-\tau/2}^{t_{n}} {\rm d}t\, (t_{n}-t)^2{\|{\dddot \psi}(t)\|}_{L^2}\Big)\\ &\le \frac{\tau}{8}\int_{t_{n-1}}^{t_{n}} {\rm d}t\, {\|{\dddot \psi}(t)\|}_{L^2}. \end{aligned} \end{equation} Analogously, for $\omega_n^{(3)}$ and $\omega_n^{(4)}$, we expand $\psi_{n-1}$ and $\psi_{n}$ around $t=t_n-\tau/2$ up to first order in $\tau/2$, \begin{gather} \label{omega3Estimate} {\|\omega_n^{(3)}\|}_{L^2} \le \frac{\tau}{4}\int_{t_{n-1}}^{t_{n}} {\rm d}t\, {\|\Delta{\ddot \psi}(t)\|}_{L^2},\\ \label{omega4Estimate} {\|\omega_n^{(4)}\|}_{L^2} \le \frac{\tau}{4}\,{\|v\|}_{L^\infty}\int_{t_{n-1}}^{t_{n}} {\rm d}t\, {\|{\ddot \psi}(t)\|}_{L^2}. \end{gather} For $\omega_n^{(5)}$, expanding $\psi_{n-1}$ and $\psi_{n}$ around $t=0$ up to zeroth order in time, we get, analogously to the estimate of $\omega_n^{(1)}$, \begin{equation} \label{omega5Estimate} {\|\omega_n^{(5)}\|}_{L^2} \le C_R h^2 {\|v\|}_{L^\infty} \Big({\|\psi_0\|}_{H^2} +\int_{0}^{t_{n}} {\rm d}t\,{\|{\dot \psi}(t)\|}_{H^2}\Big). \end{equation} Finally, for $\omega_n^{(6)}$, we apply the local Lipschitz continuity \eqref{cond_dd} to get \begin{equation} \label{omega6Estimate1} {\|\omega_n^{(6)}\|}_{L^2} \le c_1\,{\|\psi_{n-1/2}-\Psi_{n-1/2}\|}_{L^2}, \end{equation} where we used the continuum mass conservation \eqref{cp:conservation} and the coherent fully discrete mass conservation \eqref{cfdMassConservation} to define the constant \begin{equation} \label{def:c1} c_1 :=2 {\|\lambda\|}_{L^\infty}{\|V\|}_{L^\infty(\mathbb{R}^d)} \left(\mathcal{M}[\psi_0]+\mathcal{M}[\psi_{0h}]\right). \end{equation} Since we want to reinsert the decomposition \eqref{decomposition} into the right-hand side of \eqref{omega6Estimate1}, we write \begin{equation} \label{insert} \begin{aligned} {\|\psi_{n-1/2}-\Psi_{n-1/2}\|}_{L^2} &\le {\|\psi_{n-1/2}-\tfrac12(\psi_{n}+\psi_{n-1})\|}_{L^2}\\ &\quad +\,\frac12 \left({\|\rho_{n-1}\|}_{L^2}+{\|\rho_{n}\|}_{L^2} +{\|\theta_{n-1}\|}_{L^2}+{\|\theta_{n}\|}_{L^2}\right). \end{aligned} \end{equation} Plugging the estimates \eqref{RhoEstimate}, \eqref{omega1Estimate} to \eqref{omega6Estimate1}, and \eqref{insert} into \eqref{thetaEstimate1}, we find \begin{equation} \label{ThetaClosed} {\|\theta_{n}\|}_{L^2} {\le \|\theta_{n-1}\|}_{L^2} +\big(A^{(1)}_n+\tau A^{(2)}_n\big)h^2 +A^{(3)}_n\tau^2 +\frac{c_1}{2}\,\tau \left({\|\theta_{n-1}\|}_{L^2} +{\|\theta_{n}\|}_{L^2}\right), \end{equation} where the first term on the right-hand side of \eqref{insert} was estimated as in $\omega_n^{(3)}$ or $\omega_n^{(4)}$, and \begin{gather} \label{def:A1} A^{(1)}_n :=C_R\int_{t_{n-1}}^{t_{n}} {\rm d}t\,{\|\dot \psi(t)\|}_{H^2},\\ \label{def:A2} A^{(2)}_n := C_R \left(c_1+{\|v\|}_{L^\infty}\right) \Big({\|\psi_0\|}_{H^2} +\int_{0}^{t_{n}} {\rm d}t\,{\|\dot\psi(t)\|}_{H^2}\Big),\\ \label{def:A3} \begin{aligned} A^{(3)}_n &:= \frac14\left(c_1+{\|v\|}_{L^\infty}\right) \int_{t_{n-1}}^{t_{n}} {\rm d}t\, {\|\ddot\psi(t)\|}_{L^2} +\frac18\int_{t_{n-1}}^{t_{n}} {\rm d}t\, {\|\dddot\psi(t)\|}_{L^2}\\ &\quad +\frac14\int_{t_{n-1}}^{t_{n}} {\rm d}t\, {\|\Delta\ddot\psi(t)\|}_{L^2}. \end{aligned} \end{gather} If we choose the time discretization scale $\tau$ to be small enough, e.g. $c_1\tau\le 1$, we can construct the following recursive bound on ${\|\theta_{n}\|}_{L^2}$ from inequality \eqref{ThetaClosed}, \begin{equation} \label{theta_n} {\|\theta_{n}\|}_{L^2} \le B^{(1)}{\|\theta_{n-1}\|}_{L^2}+B^{(2)}_n, \end{equation} where we used that $1/(1-c_1\tau/2)\le 1+c_1\tau$ if $c_1\tau\le 1$ to define \begin{gather*} B^{(1)} :=(1+c_1\tau)^2,\\ B^{(2)}_n :=(1+c_1\tau)\big( \big(A^{(1)}_n +\tau A^{(2)}_n\big)h^2+A^{(3)}_n\tau^2\big). \end{gather*} Therefore, if we iterate the bound \eqref{theta_n} until we arrive at ${\|\theta_0\|}_{L^2}$, we get \begin{equation} \label{theta_n2} \begin{aligned} {\|\theta_n\|}_{L^2} &\le \left(B^{(1)}\right)^n{\|\theta_0\|}_{L^2} +\sum_{k=1}^{n}\left(B^{(1)}\right)^{n-k} B^{(2)}_ {k}\\ &\le c_2\Big({\|\theta_0\|}_{L^2} +\sum_{k=1}^{n}\big(\big(A^{(1)}_{k} +\tau A^{(2)}_{k}\big)h^2 +A^{(3)}_{k}\tau^2 \big)\Big), \end{aligned} \end{equation} where, on the second line of \eqref{theta_n2}, we first extract the global factor $(B^{(1)})^n$ which can then be estimated as $(B^{(1)})^n\le (1+c_1 T/N)^{2N}\le c_2$ with the definition \begin{equation} \label{def:c2} c_2 :={\rm e}^{2c_1 T}. \end{equation} It remains to estimate ${\|\theta_0\|}_{L^2}$ on the right-hand side of \eqref{theta_n2}. This is again done by using the estimate on the Ritz projection \eqref{RitzEstimate}, \begin{equation} \label{theta_0} {\|\theta_0\|}_{L^2} \le {\|\psi_0-\psi_{0h}\|}_{L^2} +{\|(R_h-1)\psi_0\|}_{L^2} \le {\|\psi_0-\psi_{0h}\|}_{L^2} +C_Rh^2{\|\psi_0\|}_{H^2}. \end{equation} Hence, with estimate \eqref{RhoEstimate} on ${\|\rho_n\|}_{L^2}$ and the estimates \eqref{theta_n2} and \eqref{theta_0} on ${\|\theta_n\|}_{L^2}$ in the decomposition \eqref{def:RhoTheta}, taking the maximum over all times, we finally arrive at \begin{equation} \label{cfd:max} \max_{n\in\mathcal{N}_0}{\|\psi_n-\Psi_n\|}_{L^2} \le c_2{\|\psi_0-\psi_{0h}\|}_{L^2}+c_3 h^2+c_4\tau^2, \end{equation} where the constants $c_3$ and $c_4$ are defined by \begin{gather} \label{cfd:c3} c_3 := C_R\left(1+c_2\left(1+T \left(c_1 +{\|v\|}_{L^\infty}\right)\right)\right) \Big({\|\psi_0\|}_{H^2}+\int_0^T {\rm d}t\,{\|\dot\psi(t)\|}_{H^2}\Big), \\ \label{cfd:c4} c_4 :=\frac{c_2}{4}\,\Big((c_1+{\|v\|}_{L^\infty}) \int_0^T {\rm d}t\,{\|\ddot\psi(t)\|}_{L^2} +\frac12 \int_0^T {\rm d}t\,{\|\dddot\psi(t)\|}_{L^2} +\int_0^T {\rm d}t\,{\|\Delta\ddot\psi(t)\|}_{L^2}\Big). \end{gather} The constants $c_1$ and $c_2$ are given in \eqref{def:c1} and \eqref{def:c2}, respectively. Using Assumption \ref{ass:initial} and setting $C_K:=\max\{c_2C_0+c_3,c_4\}$ brings the proof of Theorem \ref{thm:cfderror} to an end. \end{proof} \begin{remark} \label{rem:bound} \rm We can compute an explicit bound on the integrands in \eqref{cfd:c3} and \eqref{cfd:c4} on any finite time interval. As an example, for the last term in the constant $c_4$ from \eqref{cfd:c4}, we have (see Theorem \ref{thm:cs}(c) and \cite[p.299]{RS2}) \[ {\|\Delta\ddot\psi(t)\|}_{L^2} \le {\|\Delta^3\psi(t)\|}_{L^2} +{\|\Delta^2J[\psi(t)]\|}_{L^2} +{\|\Delta \tfrac{{\rm d}}{{\rm d}t}J[\psi(t)]\|}_{L^2}. \] The first term is exponentially bounded in time using the conditions from Theorem \ref{thm:cs} (a) and (b) and Gr\"onwall's lemma on the Duhamel integral form of the differential equation \eqref{cp}, \[ {\|\Delta^3\psi(t)\|}_{L^2} \le {\|\Delta^3\psi_0\|}_{L^2}\, {\rm e}^{C(\mathcal{M}[\psi_0])t}, \] where the constant $C$ stems from \eqref{cond_cc} (see also \cite[p.300]{RS2}). In contradistinction to the general case from Theorem \ref{thm:cs}(a), the growth rate $C$ from \eqref{cond_cc} only depends on the mass of the initial condition (and on $v$, $V$, and $\lambda$, of course). The second term is again bounded due to \eqref{cond_c}. Finally, the third term is bounded due to equation \eqref{deriv} for the time derivative of the nonlinear term $J[\psi(t)]$ and the corresponding estimates \eqref{est1} and \eqref{est2}. \end{remark} \subsection{Incoherent scheme} As described at the beginning of Section \ref{sec:fd}, we also study a second discretization scheme which approximizes the nonlinear term $g_V[|\psi|^2]\psi$ not by \eqref{gCoherent} but rather by the expression \eqref{gIncoherent}. \begin{definition} \label{def:ifd} \rm Let Assumptions \ref{ass:Vv} and \ref{ass:Sh} hold. We call $\Psi\in S_h^{\times (N+1)}$ an incoherent fully discrete solution of the Hartree initial-boundary value problem \eqref{def:cp} with initial condition $\psi_{0h}\in S_h$ if \begin{equation} \label{ifd} \begin{gathered} \begin{aligned} \mathop{\rm i}\,{(\varphi,\dot\Psi_n)}_{L^2} &={(\nabla\varphi,\nabla\Psi_{n-1/2})}_{L^2} +{(\varphi,v\Psi_{n-1/2})}_{L^2}\\ &\quad + \big(\varphi,\tfrac12\lambda g_V[|\Psi_n|^2 +|\Psi_{n-1}|^2]\Psi_{n-1/2}\big)_{L^2}, \quad \forall\varphi\in S_h,\; n\in\mathcal{N}, \end{aligned}\\ \Psi_0=\psi_{0h}. \end{gathered} \end{equation} \end{definition} The incoherent solution has the following conservation properties. \begin{proposition} \label{prop:ifdMassConservation} Let $\Psi\in S_h^{\times (N+1)}$ be an incoherent fully discrete solution of the Hartree initial-boundary value problem \eqref{def:cp} in the sense of Definition \ref{def:ifd}, and let Assumption \ref{ass:even} hold. Then, the mass and the energy of $\Psi$ are conserved under the discrete time evolution, \begin{equation} \label{ifd:conservation} \begin{gathered} \mathcal{M}[\Psi_n] = \mathcal{M}[\psi_{0h}],\quad \forall n\in\mathcal{N}_0,\\ \mathcal{H}[\Psi_n] = \mathcal{H}[\psi_{0h}],\quad \forall n\in\mathcal{N}_0. \end{gathered} \end{equation} \end{proposition} \begin{proof} Plugging $\varphi=\Psi_{n-1/2}$ into \eqref{ifd} and taking the imaginary part of the resulting equation leads to the mass conservation as in the proof of Proposition \ref{prop:cfdMassConservation}. In order to prove the energy conservation, we plug $\varphi=\dot\Psi_n$ into \eqref{ifd} and take the real part of the resulting equation. Using that ${(\Psi_n,\lambda g_V[|\Psi_{n-1}|^2]\Psi_n)}_{L^2} ={(\Psi_{n-1},\lambda g_V[|\Psi_n|^2]\Psi_{n-1})}_{L^2}$ due to Assumption \ref{ass:even}, we get \[ 0=\mathop{\rm Re}\mathop{\rm i}{(\dot \Psi_n,\dot \Psi_n)}_{L^2} =\frac{1}{2\tau}\left(\mathcal{H}[\Psi_n]-\mathcal{H}[\Psi_{n-1}]\right). \] \end{proof} We next turn to the proof of existence-uniqueness of the incoherent solution. \begin{theorem} \label{thm:ifdEx} Let Assumptions \ref{ass:Vv}, \ref{ass:Sh}, and \ref{ass:inverse} hold, and let the time discretization scale $\tau$ be sufficiently small. Then, there exists a unique incoherent fully discrete solution of the Hartree initial-boundary value problem \eqref{def:cp} in the sense of Definition \ref{def:ifd}. \end{theorem} \begin{proof} The proof for the incoherent solution is analogous to the proof of the coherent solution. Let $\phi\in S_h$ be given, and define the mapping $G_\phi: S_h\to S_h$ by \begin{equation} \label{def:Gphi} \begin{aligned} {(\varphi,G_\phi[\psi])}_{L^2} &:= {(\varphi,\phi)}_{L^2} -\frac{\mathop{\rm i}\tau}{2}\Big( {(\nabla\varphi,\nabla\psi)}_{L^2} +{(\varphi,v\psi)}_{L^2}\\ &\quad +{(\varphi,\tfrac12\lambda g_V[|2\psi-\phi|^2+|\phi|^2]\psi)}_{L^2} \Big),\quad \forall\varphi\in S_h. \end{aligned} \end{equation} For some $n\in\mathcal{N}$, let the $n$-th component $\Psi_{n-1}$ of $\Psi\in S_h^{\times (N+1)}$ from \eqref{def:Psi} be given. Adding $2\mathop{\rm i}(\varphi,\Psi_{n-1})/\tau$ on both sides of \eqref{ifd}, we rewrite \eqref{ifd} with the help of \eqref{def:Gphi} in the form of a fixed point equation for $\Psi_{n-1/2}$, \[ \Psi_{n-1/2} =G_{\Psi_{n-1}}[\Psi_{n-1/2}]. \] To use the Banach's fixed point theorem as in the proof of Theorem \ref{prop:ExistenceCFD}, we show that $G_{\Psi_{n-1}}$ maps the compact Ball $\mathcal{B}_{n-1}:=\{\psi\in S_h\,|\, {\|\psi\|}_{L^2}\le \mathcal{M}[\Psi_{n-1}]^{1/2}+1\}$ into itself and that $G_{\Psi_{n-1}}$ is a strict contraction on $\mathcal{B}_{n-1}$. To this end, we write \begin{equation} \label{DiffGn1} \begin{aligned} &|{(\varphi,G_{\Psi_{n-1}}[\psi]-G_{\Psi_{n-1}}[\xi])}_{L^2}|\\ &\le\frac{\tau}{2}\Big( {\|\nabla(\psi-\xi)\|}_{L^2} {\|\nabla\varphi\|}_{L^2} +{\|v\|}_{L^\infty} {\|\psi-\xi\|}_{L^2} {\|\varphi\|}_{L^2} +\,\frac12\,A{\|\varphi\|}_{L^2} \Big), \end{aligned} \end{equation} where, with the help of \eqref{est}, \eqref{dc1}, and $||z|^2-|w|^2|\le |z+w||z-w|$ for all $z,w\in\mathbb{C}$, the third term $A$ on the right-hand side of \eqref{DiffGn1} can be estimated as \begin{align*} A &:= {\|\lambda g_V[|2\psi-\Psi_{n-1}|^2+|\Psi_{n-1}|^2]\psi -\lambda g_V[|2\xi-\Psi_{n-1}|^2+|\Psi_{n-1}|^2]\xi\|}_{L^2}\\ &\le 8{\|\lambda\|}_{L^\infty}{\|V\|}_{L^\infty(\mathbb{R}^d)}\, \big({\|\psi\|}_{L^2}^2 +{\|\xi\|}_{L^2}^2 +{\|\Psi_{n-1}\|}_{L^2}^2\big) {\|\psi-\xi\|}_{L^2}. \end{align*} Hence, plugging $\varphi=G_{\Psi_{n-1}}[\psi]-G_{\Psi_{n-1}}[\xi]$ into \eqref{DiffGn1}, we get for $\psi,\xi\in\mathcal{B}_{n-1}$ using Assumption \ref{ass:inverse}, \[ {\|G_{\Psi_{n-1}}[\psi]-G_{\Psi_{n-1}}[\xi]\|}_{L^2} \le \frac{\alpha_{n-1}}{2}\,\tau\,{\|\psi-\xi\|}_{L^2}, \] where \[ \alpha_{n-1}:=C_B^2h^{-2} +{\|v\|}_{L^\infty} +12{\|\lambda\|}_{L^\infty}{\|V\|}_{L^\infty(\mathbb{R}^d)}\, (\mathcal{M}[\Psi_{n-1}]^{1/2}+1)^2 \] like in the coherent scheme \eqref{taun}. Therefore, we arrive at the claim as in the proof of Proposition \ref{prop:ExistenceCFD} using the mass conservation from \eqref{ifd:conservation}, i.e., the incoherent solution exists and is unique if the time discretization scale $\tau$ is sufficiently small, e.g. $\alpha_0(\mathcal{M}[\psi_{0h}]^{1/2}+1)\tau\le 1$. \end{proof} Finally, we also provide a time quadratic accuracy estimate on the $L^2$-error of the incoherent solution. Again, the proof is analogous to the corresponding proof for the coherent solution from Theorem \ref{thm:cfderror}. \begin{theorem} \label{thm:ifderror} Let Assumptions \ref{ass:Vv}, \ref{ass:Sh}, \ref{ass:AccuracyOrder}, and \ref{ass:initial} hold, and let $\Psi\in S_h^{\times(N{+}1)}$ be the incoherent solution from Theorem \ref{thm:ifdEx}. Then, there exists a constant $C_I>0$ such that \[ \max_{n\in \mathcal{N}_0}{\|\psi(t_n)-\Psi_n\|}_{L^2} \le C_I (\tau^2+h^2). \] \end{theorem} \begin{proof} As in \eqref{decomposition} and \eqref{def:RhoTheta}, we make use of the decomposition $\psi_n-\Psi_n=\rho_n+\theta_n$, and we estimate $\rho_n$ again by \eqref{RhoEstimate}. In order to estimate ${\|\theta_n\|}_{L^2}$, we define $L_{n,\varphi}$ as in \eqref{def:Ln}. Then, everything in equation \eqref{Ln1} remains unchanged up to the last term which is replaced by the expression $-(\varphi,\tfrac12 \lambda g_V[|\Psi_{n}|^2+|\Psi_{n-1}|^2]\Psi_{n-1/2})_{L^2}$. Using again \eqref{Rh}, we can rewrite $L_{n,\varphi}$ as in equation \eqref{def:omegas}, where the terms $\omega_n^{(j)}$ remain unchanged for all $j=1,\dots ,5$, whereas the term $\omega_n^{(6)}$ now has the form \[ \omega_n^{(6)} :=f[\psi_{n-1/2}]-\lambda g_V[\tfrac12(|\Psi_n|^2+|\Psi_{n-1}|^2)] \Psi_{n-1/2} = \omega_{n,1}^{(6)} +\omega_{n,2}^{(6)} +\omega_{n,3}^{(6)}, \] where we use the same notation as introduced after \eqref{Rh} to define \begin{gather*} \omega_{n,1}^{(6)} :=\lambda g_V[|\psi_{n-1/2}|^2](\psi_{n-1/2}-\tfrac12 (\psi_n+\psi_{n-1})),\\ \omega_{n,2}^{(6)} :=\frac12\, \lambda g_V[|\psi_{n-1/2}|^2-\tfrac12(|\psi_n|^2+|\psi_{n-1}|^2)] (\psi_n+\psi_{n-1}),\\ \omega_{n,3}^{(6)} :=\frac12\, \lambda g_V[\tfrac12(|\psi_n|^2+|\psi_{n-1}|^2)] (\psi_n+\psi_{n-1}) - \lambda g_V[\tfrac12(|\Psi_n|^2+|\Psi_{n-1}|^2)]\Psi_{n-1/2}. \end{gather*} For $\omega_{n,1}^{(6)}$, expanding $\psi_{n-1}$ and $\psi_{n}$ around $t=t_n-\tau/2$ up to first order in $\tau/2$, we get , using \eqref{est} and the mass conservation \eqref{cp:conservation}, \begin{equation} \label{omega6_1} {\|\omega_{n,1}^{(6)}\|}_{L^2} \le \frac{\tau}{4}\,{\|\lambda\|}_{L^\infty}{\|V\|}_{L^\infty(\mathbb{R}^d)} \mathcal{M}[\psi_0] \int_{t_{n-1}}^{t_n} {\rm d}t\, {\|\ddot\psi(t)\|}_{L^2}. \end{equation} To estimate $\omega_{n,2}^{(6)}$, we expand $\psi_{n-1}$ and $\psi_{n}$ around $t=t_n-\tau/2$ up to first order in $\tau/2$ and get similarly \begin{equation} \label{omega6_2} \begin{aligned} {\|\omega_{n,2}^{(6)}\|}_{L^2} &\le {\|\lambda\|}_{L^\infty}{\|V\|}_{L^\infty(\mathbb{R}^d) } \mathcal{M}[\psi_0]^{1/2}{\| |\psi_{n-1/2}|^2 -\tfrac12 (|\psi_n|^2+|\psi_{n-1}|^2)\|}_{L^1}\\ &\le a^{(0)}(a_n^{(1)} \tau+a_n^{(2)} \tau^2+a_n^{(3)} \tau^3), \end{aligned} \end{equation} where we define \begin{gather*} a^{(0)} := {\|\lambda\|}_{L^\infty}{\|V\|}_{L^\infty(\mathbb{R}^d) } \mathcal{M}[\psi_0]^{1/2},\\ a_n^{(1)} := \mathcal{M}[\psi_0]\int_{t_{n-1}}^{t_n} {\rm d}t\, {\|\ddot\psi(t)\|}_{L^2}\\ a_n^{(2)} :=\frac12 \,{\|\dot\psi_{n-1/2}\|}_{L^2}\Big({\|\dot\psi_{n-1/2}\|}_{L^2} +\int_{t_{n-1}}^{t_n} {\rm d}t\, {\|\ddot\psi(t)\|}_{L^2}\Big)\\ a_n^{(3)} := \frac{1}{8}\int_{t_{n-1}}^{t_n} {\rm d}t\, {\|\ddot\psi(t)\|}_{L^2}^2. \end{gather*} For $\omega_{n,3}^{(6)}$, using in particular again $||z|^2-|w|^2|\le |z+w| |z-w|$ for all $z,w\in\mathbb{C}$ and the decomposition \eqref{decomposition}, we get \begin{equation} \label{omega6_3} \begin{aligned} {\|\omega_{n,3}^{(6)}\|}_{L^2} &\leq \frac14\, {\|\lambda g_V[|\psi_n|^2+|\psi_{n-1}|^2](\psi_n-\Psi_n +\psi_{n-1}- \Psi_{n-1})\|}_{L^2}\\ &\quad +\frac14\,{\|\lambda g_V[|\psi_n|^2-|\Psi_n|^2+|\psi_{n-1}|^2-|\Psi_{n-1}|^2] (\Psi_n+\Psi_{n-1})\|}_{L^2}\\ &\hspace{-1.4cm}\leq {\|\lambda\|}_{L^\infty} {\|V\|}_{L^\infty(\mathbb{R}^d)}(\mathcal{M}[\psi_0] +\mathcal{M}[\psi_{0h}]) ({\|\rho_n\|}_{L^2}+{\|\rho_{n-1}\|}_{L^2}+{\|\theta_n\|}_{L^2}+{\|\theta_{n-1}\|}_{L^2}). \end{aligned} \end{equation} Therefore, plugging the estimates \eqref{omega6_1}, \eqref{omega6_2}, and \eqref{omega6_3} into \eqref{thetaEstimate1}, we again find the closed inequality \eqref{ThetaClosed}, the coefficients $A_n^{(1)}$ and $A_n^{(2)}$ having the same form as in \eqref{def:A1} and \eqref{def:A2}, respectively. Using estimate \eqref{omega6_1} on $\omega_{n,1}^{(6)}$, we see that the coefficient $A_n^{(3)}$ in the incoherent case contains all the terms from \eqref{def:A3} of the coherent case with $c_1$ replaced by $a^{(0)}$, plus an additional term of the form $a^{(0)}(a_n^{(1)}+a_n^{(2)}\tau+a_n^{(3)}\tau^2)$ which is due to the estimate \eqref{omega6_2} of $\omega_{n,2}^{(6)}$. Plugging the coefficients $A_n^{(1)}$, $A_n^{(2)}$, and $A_n^{(3)}$ into the iterated bound \eqref{theta_n2} and using estimate \eqref{theta_0} on $\theta_0$, we again get an estimate of the form \eqref{cfd:max} where the constant $c_3$ has the same form as in \eqref{cfd:c3} whereas the constant $c_4$, compared to \eqref{cfd:c4}, now looks like \begin{equation}\label{ifd:c4} \begin{aligned} c_4 &:=\frac{c_2}{4} \Big((a^{(0)}+{\|v\|}_{L^\infty}) \int_0^T {\rm d}t\,{\|\ddot\psi(t)\|}_{L^2} +\frac12 \int_0^T {\rm d}t\,{\|\dddot \psi(t)\|}_{L^2} +\int_0^T {\rm d}t\,{\|\Delta\ddot\psi(t)\|}_{L^2}\Big)\\ &\quad +c_2a^{(0)}\big(\mathcal{M}[\psi_0] +\tfrac{1}{2}\,T\max_{t\in[0,T]}{\|\dot\psi(t)\|}_{L^2}\big)\int_0^T {\rm d}t\,{\|\ddot\psi(t)\|}_{L^2}\\ &\quad +\frac{c_2a^{(0)}}{2}T\max_{t\in[0,T]}{\|\dot\psi(t)\|}_{L^2}^2 +\frac{c_2a^{(0)}}{8}T^2 \int_0^T {\rm d}t\,{\|\ddot\psi(t)\|}_{L^2}^2. \end{aligned} \end{equation} Herewith, as in the proof of Theorem \ref{thm:cfderror}, we arrive at the assertion. \end{proof} \begin{remark} \label{rmk38} \rm Using estimates as in Remark \ref{rem:bound}, we can again bound the constants in \eqref{ifd:c4} explicitly. \end{remark} \begin{thebibliography}{99} \bibitem{ADK91} Akrivis, G. D.; Dougalis, V. A.; Karakashian, O. A.; 1991 On fully discrete Galerkin methods of second-order temporal accuracy for the nonlinear Schr\"odinger equation \emph{Numer. Math.} \textbf{59} 31--53 \bibitem{diss} Aschbacher, W. H.; 2001 Large systems of non-relativistic Bosons and the Hartree equation \emph{Diss. ETH No. 14135} \bibitem{A09} Aschbacher, W. H.; On radiative dissipation of Hartree solitons, in progress \bibitem{Aschbacher} Aschbacher, W. H.; Fr\"ohlich, J.; Graf, G. M.; Schnee, K.; Troyer, M.; 2002 Symmetry breaking regime in the nonlinear Hartree equation \emph{J. Math. Phys.} \textbf{43} 3879--91 \bibitem{Aschbacher08} Aschbacher, W. H.; Squassina, M.; 2008 On phase segregation in nonlocal two-particle Hartree systems, submitted, and arXiv:0809.3369 \bibitem{Brenner} Brenner, S. C.; Scott, L. R.; 1994 \emph{The mathematical theory of finite element methods} (New York: Springer) \bibitem{Brezis} Brezis, H.; Gallouet, T.; 1980 Nonlinear Schr\"odinger evolution equations \emph{Nonlinear Anal.} \textbf{4} 677--81 \bibitem{Squassina} Caliari, M.; Squassina, M.; 2008 Location and phase segregation of ground and excited states for 2D Gross-Pitaevskii systems \emph{Dyn. Partial Differential Equations} \textbf{5} 117--37 \bibitem{DFP81} Delfour, M.; Fortin, M.; Payre, G.; 1981 Finite-difference solutions of a non-linear Schr\"odinger equation \emph{J. Comput. Phys.} \textbf{44} 277--88 \bibitem{Elgart} Elgart, A.; Erd\"os, L.; Schlein, B; Yau, H. T.; 2004 Nonlinear Hartree equation as the mean field limit of weakly coupled fermions \emph{J. Math. Pures Appl.} \textbf{83}, 1241--73 \bibitem{Erdos} Erd\"os, L.; Schlein, B.; Yau, H. T.; 2007 Rigorous derivation of the Gross-Pitaevskii equation \emph{Phys. Rev. Lett.} \textbf{98} 040404 \bibitem{Froehlich1} Fr\"ohlich, J.; Tsai, T. P.; Yau, H. T.; 2000 On a classical limit of quantum theory and the nonlinear Hartree equation \emph{Geom. funct. anal.} 57--78 \bibitem{Froehlich2} Fr\"ohlich, J.; Gustafson, S.; Jonsson, B. L. G.; Sigal, I. M.; 2004 Solitary wave dynamics in an external potential \emph{Commun. Math. Phys.} \textbf{250} 613--42 \bibitem{Lieb} Lieb, E.; Seiringer, R.; Yngvason, J.; 2000 Bosons in a trap: a rigorous derivation of the Gross-Pitaevskii energy functional \emph{Phys. Rev. A} \textbf{61} 043602 \bibitem{Onorato} Onorato, M.; Osborne, A. R.; Serio, M.; 2006 Modulational instability in crossing sea states: a possible mechanism for the formation of freak waves \emph{Phys. Rev. Lett} \textbf{96} 014503 \bibitem{RS2} Reed, M.; Simon, B.; 1975 \emph{Methods of modern mathematical physics II} (New York: Academic Press) \bibitem{Sanz-Serna} Sanz-Serna, J. M.; Verwer, J. G.; 1986 Conservative and nonconservative schemes for the solution of the nonlinear Schr\"odinger equation \emph{IMA J. Numer. Anal.} \textbf{6} 25--42 \bibitem{Shukla} Shukla, P. K.; Kourakis, I.; Eliasson, B.; Marklund, M.; Stenflo, L.; 2006 Instability and evolution of nonlinearly interacting water waves \emph{Phys. Rev. Lett.} \textbf{97} 094501 \bibitem{Spohn} Spohn, H.; 1980 Kinetic equations from Hamiltonian dynamics \emph{Rev. Mod. Phys.} \textbf{52} 569--615 \bibitem{VT97} Thom\'ee, V.; 1997 \emph{Galerkin finite element methods for parabolic problems} Springer Series in Computational Mathematics, Vol. 25 (Berlin: Springer) \bibitem{Tourigny} Tourigny, Y.; Morris, J. L.; 1988 An investigation into the effect of product approximation in the numerical solution of the cubic nonlinear Schr\"odinger equation \emph{J. Comput. Phys.} \textbf{76} 103--30 \bibitem{Zeidler} Zeidler, E.; 1990 \emph{Nonlinear functional analysis and its applications II/A} (New York: Springer) \end{thebibliography} \end{document}