Separation of scales: dynamical approximations for composite quantum systems^*

The quantum system is described by a time-dependent Schrödinger equation

$$\begin{equation}\mathrm{i}{\partial }_{t}\psi =H\psi ;\quad {\psi }_{\vert t=0}={\psi }_{0},\end{equation} \tag{ 1 }$$

governed by a Hamiltonian of the form

$$\begin{equation*}H={H}_{x}+{H}_{y}+W(x,y),\end{equation*}$$

where the coupling potential W(x, y) is a smooth function, that satisfies growth estimates guaranteeing existence and uniqueness of the solution to the Schrödinger equation (1) for a rather general set of initial data, as we shall see later in section 2. The overall set of space variables is denoted as $(x,y)\in {\mathbb{R}}^{n}\times {\mathbb{R}}^{d}$ such that the total dimension of the configuration space is n + d. The wave function depends on time t ⩾ 0 and both space variables, that is, ψ = ψ(t, x, y). We suppose that initially scales are separable, that is, we work with initial data of product form

$$\begin{equation}\psi (0,x,y)={\psi }_{0}(x,y)={\varphi }_{0}^{x}(x){\varphi }_{0}^{y}(y).\end{equation} \tag{ 2 }$$

In the simple case without coupling, that is, W ≡ 0, the solution stays separated, ψ(t, x, y) = φ^x (t, x)φ^y (t, y) for all time, where

$$\begin{align*}\mathrm{i}{\partial }_{t}{\varphi }^{x}={H}_{x}{\varphi }^{x};\quad {\varphi }_{\vert t=0}^{x}={\varphi }_{0}^{x},\\ \mathrm{i}{\partial }_{t}{\varphi }^{y}={H}_{y}{\varphi }^{y};\quad {\varphi }_{\vert t=0}^{y}={\varphi }_{0}^{y},\end{align*}$$

and this is an exact formula. Here, we aim at investigating the case of an actual coupling with ∂_x ∂_y W ≢ 0 and look for approximate solutions of the form

$$\begin{equation*}{\psi }_{\mathrm{app}}(t,x,y)={\varphi }^{x}(t,x){\varphi }^{y}(t,y),\end{equation*}$$

where the individual components satisfy evolution equations that account for the coupling between the variables. The main motivation for such approximations is dimension reduction, since φ^x (t, x) and φ^y (t, y) depend on variables of lower dimension than the initial (x, y). Of crucial importance is the choice of the approximate Hamiltonian H_app = H_x + H_y + W_app(t, x, y), that governs the approximate dynamics. We consider two different approximate coupling potentials: one is the time-dependent Hartree mean-field potential, the other one, computationally less demanding, is based on a brute force single-point collocation. Time-dependent Hartree methods have been known for a long time and have earned the reputation of oversimplifying the dynamics of real molecular systems. We emphasize, that our present study does not aim at rejuvenating but at deriving rigorous mathematical error estimates, which seem to be missing in the literature. Surprisingly, our error analysis provides similar estimates for both methods, the collocation and the mean-field approach. We investigate the size of the difference between the true and the approximate solution in the L²-norm,

$$\begin{equation*}{\Vert}\psi (t)-{\psi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}=\sqrt{{\int }_{{\mathbb{R}}^{n+d}}\vert \psi (t,x,y)-{\psi }_{\mathrm{app}}(t,x,y){\vert }^{2}\enspace \mathrm{d}x\enspace \mathrm{d}y}\end{equation*}$$

and in Sobolev norms. We present error estimates that explicitly depend on derivatives of the coupling potential W(x, y) and on moments of the approximate solution. As an additional error measure we also consider the deviation of true and approximate expectation values $\left\langle \psi (t),A\psi (t)\right\rangle -\left\langle {\psi }_{\mathrm{app}}(t),A{\psi }_{\mathrm{app}}(t)\right\rangle$, for self-adjoint linear operators A. Roughly speaking, the estimates we obtain for observables depend on one more derivative of the coupling potential W(x, y) than the norm estimates. This means that in many situations expectation values are more accurately described than the wave function itself. Even though rigorous error estimates that quantify the decoupling of quantum subsystems in terms of flatness properties of the coupling potential W(x, y) are naturally important, our results here seem to be the first ones of their kind.

Interacting quantum systems have traditionally been formulated from the point of view of reduced dynamics theories, based on quantum master equations in a Markovian or non-Markovian setting [1]. More recently, also tensorized representations of the full quantum system have been considered, as for example by matrix product states [2, 3] or within a multiconfiguration time-dependent Hartree (MCTDH) approach [4–7]. Both wavefunctions (pure states) and density operators (mixed states) can be described in this framework, and wavefunction-based computations can be used to obtain density matrices [8]. In the chemical physics literature, dimension reduction for quantum systems has been proposed in the context of mean-field methods [9, 10], and the quantum–classical mean-field Ehrenfest approach [11, 12]. Also, quantum–classical formulations have been derived in Wigner phase space [13, 14] and in a quantum hydrodynamic setting [15–17]. Our present mathematical formulation circumvents formal difficulties of these approaches [18–20], by preserving a quantum wavefunction description for the entire system. Previous mathematical work we are aware of is concerned with rather specific coupling models, as for example the coupling of Hartree–Fock and classical equations in [21], or the time-dependent self-consistent field equations in [22], or with adiabatic approximations which rely on eigenfunctions for one part of the system, see for example [23, 24]. To the best of our knowledge, the rather general mathematical analysis of scale separation in quantum systems we are developing here is new.

For a first approximate potential, we consider a brute-force approach, where we collocate partially at a single point, for definiteness we choose the origin, and set

$$\begin{equation*}{W}_{\text{bf}}(x,y)=W(x,0)+W(0,y)-W(0,0).\end{equation*}$$

In comparison, following the more conventional time-dependent Hartree approach, we set

$$\begin{equation*}{W}_{\mathrm{mf}}(t,x,y)={\left\langle W\right\rangle }_{y}(t,x)+{\left\langle W\right\rangle }_{x}(t,y)-\left\langle W\right\rangle (t),\end{equation*}$$

where we perform partial and full averages of the coupling potential,

$$\begin{align*}& {\left\langle W\right\rangle }_{x}={\int }_{{\mathbb{R}}^{d}}W(x,y)\enspace \vert {\varphi }^{x}(t,x){\vert }^{2}\enspace \mathrm{d}x/{\int }_{{\mathbb{R}}^{d}}\vert {\varphi }^{x}(t,x){\vert }^{2}\enspace \mathrm{d},\\ & {\left\langle W\right\rangle }_{y}={\int }_{{\mathbb{R}}^{d}}W(x,y)\enspace \vert {\varphi }^{y}(t,y){\vert }^{2}\enspace \mathrm{d}y/{\int }_{{\mathbb{R}}^{d}}\vert {\varphi }^{y}(t,y){\vert }^{2}\enspace \mathrm{d}y,\\ & \left\langle W\right\rangle ={\int }_{{\mathbb{R}}^{n+d}}W(x,y)\enspace \vert {\varphi }^{x}(t,x){\varphi }^{y}(t,y){\vert }^{2}\enspace \mathrm{d}x\enspace \mathrm{d}y/{\int }_{{\mathbb{R}}^{n+d}}\vert {\varphi }^{x}(t,x){\varphi }^{y}(t,y){\vert }^{2}\enspace \mathrm{d}x\enspace \mathrm{d}y.\end{align*}$$

For both approximations, the brute-force and the mean-field approximation, we derive various types of estimates for the error in L²-norm. Our key finding is that both methods come with error bounds that are qualitatively the same, since they draw from either evaluations or averages of the function

$$\begin{equation*}\delta W(x,{x}^{\prime },y,{y}^{\prime })=W(x,y)-W(x,{y}^{\prime })-W({x}^{\prime },y)+W({x}^{\prime },{y}^{\prime }).\end{equation*}$$

Depending on whether one chooses to control the auxiliary function δW in terms of ∇_x W, ∇_y W or ∇_x ∇_y W, the estimate requires a balancing with corresponding moments of the approximate solution. For example, proposition 1 provides L²-norm estimates of the form

$$\begin{equation*}\hspace{25.0pt}{\Vert}\psi (t)-{\psi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}\leqslant \begin{cases}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\enspace {\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\int }_{0}^{t}{\Vert}y{\varphi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s,\quad \\ \mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\enspace {\Vert}{\nabla }_{x}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}{\int }_{0}^{t}{\Vert}x{\varphi }^{x}(s){{\Vert}}_{{L}_{x}^{2}}{\Vert}y{\varphi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s,\quad \end{cases}\end{equation*}$$

where const. ∈ {1, 2, 4}, depending on whether ψ_app(t) results from the brute-force or the mean-field approximation. Example 3 discusses important variants of this estimate using different ways of quantifying the flatness of the coupling potential. Proposition 2 gives analogous estimates in Sobolev norms. In addition, we analyze the deviation of the true and the approximate expectation values in a similar vein. For the expectation values, we again obtain qualitatively similar error estimates for W_bf and W_mf. The upper bounds differ from the norm bounds in so far as they involve one more derivative of the coupling potential W and low order Sobolev norms of the approximate solution, see proposition 3. Hence, from the perspective of approximation accuracy, the brute force and the mean-field approach differ only slightly. Therefore, other assessment criteria are needed for explaining the prevalence of the Hartree method in many applications, as we will discuss in section 3.5.

In the second part of the paper we turn to a specific case of the previous general class of coupled Hamiltonians ${H}^{\varepsilon }={H}_{x}+{H}_{y}^{\varepsilon }+W(x,y)$ and consider for one part of the system a semiclassically scaled Schrödinger operator

$$\begin{equation*}{H}_{y}^{\varepsilon }=-\frac{{\varepsilon }^{2}}{2}{{\Delta}}_{y}+{V}_{2}(y),\quad \varepsilon > 0.\end{equation*}$$

We will discuss in section 5.1 system–bath Hamiltonians that can be recast in this semiclassical format. The initial data are still a product of the form 2, but the y-factor is chosen as

$$\begin{equation*}{\varphi }_{0}^{y}(y)={\varepsilon }^{-d/4}a\left(\frac{y-{q}_{0}}{\sqrt{\varepsilon }}\right){\mathrm{e}}^{\mathrm{i}{p}_{0}\cdot (y-{q}_{0})/\varepsilon },\end{equation*}$$

that is, ${\varphi }_{0}^{y}$ is a semiclassical wave packet with a smooth and rapidly decaying amplitude function $a\in \mathcal{S}({\mathbb{R}}^{d})$, and an arbitrary phase space centre $({q}_{0},{p}_{0})\in {\mathbb{R}}^{2d}$. We will choose a semiclassical wave packet approximation for φ^y (t, y) exploring two different choices for the centre (q(t), p(t)). As a first option we consider the classical trajectory

$$\begin{equation*}\dot{q}=p,\hspace{25.0pt}\dot{p}=-\nabla {V}_{2}(q),\end{equation*}$$

and as a second option the corresponding trajectory resulting from the averaged gradient of the potential V₂,

$$\begin{equation*}\left\langle \nabla {V}_{2}\right\rangle (t)={\int }_{{\mathbb{R}}^{d}}{\nabla }_{y}{V}_{2}(y)\enspace \vert {\varphi }^{y}(t,y){\vert }^{2}\enspace \mathrm{d}y/{\int }_{{\mathbb{R}}^{d}}\vert {\varphi }^{y}(t,y){\vert }^{2}\enspace \mathrm{d}y.\end{equation*}$$

Correspondingly, the approximative factor φ^x (t, x) is evolved by the partial Hamiltonian H_x + W_eff with

$$\begin{align*}& {W}_{\text{eff}}(t,x)=W(x,q(t))\quad \mathrm{o}\mathrm{r}\\ & {W}_{\text{eff}}(t,x)={\int }_{{\mathbb{R}}^{d}}W(x,y)\enspace \vert {\varphi }^{y}(t,y){\vert }^{2}\enspace \mathrm{d}y/{\int }_{{\mathbb{R}}^{d}}\vert {\varphi }^{y}(t,y){\vert }^{2}\enspace \mathrm{d}y.\end{align*}$$

We obtain error estimates in L²-norm, see proposition 4 and expectation values, see proposition 5. These estimates are given in terms of the semiclassical parameter ɛ and derivatives of the coupling potential. Again, both choices for the effective potential differ only slightly in approximation accuracy. Measuring the coupling strength in terms of $\eta ={\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}$, we obtain two-parameter estimates of order $\sqrt{\varepsilon }+\eta /\sqrt{\varepsilon }$ in norm, while the ones for the expectation values are of order ɛ + η. Hence, again the accuracy of quadratic observables is higher than the one for wavefunctions.

We choose a very classical set of assumptions on the regularity and the growth of the potential, since our focus is more on finding appropriate ways to approximate the solution in a standard framework than on treating specific situations.

Assumption 1. All the potentials that we consider are smooth, real-valued, and at most quadratic in their variables:

$$\begin{equation*}{V}_{1}\in {C}^{\infty }({\mathbb{R}}^{n};\mathbb{R}),\enspace {V}_{2}\in {C}^{\infty }({\mathbb{R}}^{d};\mathbb{R}),\enspace W\in {C}^{\infty }({\mathbb{R}}^{n+d};\mathbb{R}),\end{equation*}$$

and, for $\alpha \in {\mathbb{N}}^{n}$, $\beta \in {\mathbb{N}}^{d}$,

$$\begin{equation*}\begin{gathered}{\partial }_{x}^{\alpha }{V}_{1}\in {L}^{\infty }\enspace \mathrm{f}\mathrm{o}\mathrm{r}\enspace \vert \alpha \vert \geqslant 2,\hspace{25.0pt}{\partial }_{y}^{\beta }{V}_{2}\in {L}^{\infty }\enspace \mathrm{f}\mathrm{o}\mathrm{r}\enspace \vert \beta \vert \geqslant 2,\\ {\partial }_{x}^{\alpha }{\partial }_{y}^{\beta }W\in {L}^{\infty }\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace \vert \alpha \vert +\vert \beta \vert \geqslant 2.\end{gathered}\end{equation*}$$

We also assume that ∇_y W ∈ L^∞, but note that this condition can easily be relaxed, see example 3. All the initial date we consider are smooth and rapidly decaying, that is, Schwartz class functions:

$$\begin{equation*}{\varphi }_{0}^{x}\in \mathcal{S}({\mathbb{R}}^{n};\mathbb{C}),\quad {\varphi }_{0}^{y}\in \mathcal{S}({\mathbb{R}}^{d};\mathbb{C})\quad (\mathrm{h}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\enspace {\psi }_{0}\in \mathcal{S}({\mathbb{R}}^{n+d};\mathbb{C})).\end{equation*}$$

Under the above assumption, it is well-known that H_x , H_y and H are essentially self-adjoint on ${L}^{2}({\mathbb{R}}^{N})$, with N = n, d and n + d, respectively (as a consequence of Faris–Lavine theorem, see e.g. [26, theorem 10.38]).

Example 1. Since assumption 1 involves similar properties for V₁, V₂ or W, we present examples for V₁ only, which can readily be adapted to V₂ and W. For instance, we can consider

$$\begin{equation*}{V}_{1}(x)=\sum\limits _{j=1}^{n}\enspace {\omega }_{j}^{2}{x}_{j}^{2}+\sum\limits _{j=1}^{n}\enspace {\beta }_{j}\mathrm{exp}(-\left\langle x,{A}_{j}x\right\rangle )+{V}_{\text{per}}(x),\end{equation*}$$

with ω_j ⩾ 0 (possibly anisotropic harmonic potential), ${\beta }_{j}\in \mathbb{R}$, ${A}_{j}\in {\mathbb{R}}^{n\times n}$ positive definite symmetric matrices (Gaussian potential), and V_per a smooth potential, periodic along some lattice in ${\mathbb{R}}^{n}$.

The assumptions on the growth and smoothness of the potentials and the regularity of the initial data call for comments.

Remark 1. Concerning the growth of V₁, V₂ and W, the assumption that they are at most quadratic concerns the behavior at infinity and could be relaxed, up to suitable sign assumptions. Local behavior is rather free, for example a local double well is allowed, as long as it is not too confining at infinity. We choose to stick to the at most quadratic case, since bounded second order derivatives simplify the presentation.

Remark 2. Concerning the smoothness of our potentials, most of our results still hold assuming only smoothness of W, as long as the operators H_x and H_y are essentially self-adjoint on an adequate domain included in ${L}^{2}({\mathbb{R}}^{N})$, with N = n, d. For example, V₁ and V₂ could both present Coulomb singularities, and the results of proposition 1 would still hold. In the semiclassical régime, we can also allow a Coulomb singularity for V₁ and prove propositions 4 and 5.

Remark 3. Concerning the smoothness and the decay of the initial data, most of our results still hold, if the initial data are contained in one of the spaces ${{\Sigma}}^{k}({\mathbb{R}}^{N})$ containing functions f whose norms

$$\begin{equation}{\Vert}f{{\Vert}}_{{{\Sigma}}^{k}}=\underset{\genfrac{}{}{0pt}{}{z\in {\mathbb{R}}^{N}}{\vert \alpha \vert +\vert \beta \vert \leqslant k}}{\mathrm{sup}}{\Vert}{z}^{\alpha }{\partial }_{z}^{\beta }f{{\Vert}}_{{L}^{2}}\end{equation} \tag{ 4 }$$

are bounded. Note that $\mathcal{S}({\mathbb{R}}^{N})={\cap }_{k\in \mathbb{N}}{{\Sigma}}^{k}$. For example, proposition 1 still holds for initial data in Σ¹, while proposition 4 requires initial data in a semiclassically scaled Σ³ space.

An important class of coupled quantum systems are described by system–bath Hamiltonians [1, 25].

$$\begin{equation*}{H}_{\text{sb}}=-\frac{1}{2}{{\Delta}}_{x}-\frac{1}{2}{{\Delta}}_{y}+{V}_{\mathrm{s}}(x)+{V}_{\mathrm{b}}(y)+{V}_{\text{sb}}(x,y).\end{equation*}$$

These are naturally given in the format required by (26). In the present discussion, we specify that the bath is described by a harmonic oscillator, ${V}_{\mathrm{b}}(y)=\frac{1}{2}{k}_{2}^{0}\vert y{\vert }^{2}$ (or a set of harmonic oscillators in more than one dimension) and the system–bath coupling V_sb(x, y) = W(x, y) is of cubic form, such that we obtain in the notation of (26),

$$\begin{equation*}\hspace{25.0pt}{H}_{x}=-\frac{1}{2}{{\Delta}}_{x}+{V}_{\mathrm{s}}(x),\hspace{25.0pt}{H}_{y}=-\frac{1}{2}{{\Delta}}_{y}+\frac{1}{2}{k}_{2}^{0}\vert y{\vert }^{2},\hspace{25.0pt}W(x,y)=\frac{1}{2}\overrightarrow{\eta }\cdot x\vert y{\vert }^{2},\end{equation*}$$

where ${k}_{2}^{0} > 0$ and $\overrightarrow{\eta }\in {\mathbb{R}}^{d}$. The cubic, anharmonic coupling W(x, y) is a non-trivial case, which is employed, e.g. in the description of vibrational dephasing [27, 28] and Fermi resonances [29]. It is natural to assume smoothness and subquadratic growth for V_s(x). However, the coupling potential W(x, y) clearly fails to satisfy the growth condition of assumption 1. Moreover, it is not clear that in such a framework the total Hamiltonian H is essentially self-adjoint. On the other hand, adding a quartic confining potential,

$$\begin{equation*}{H}_{y}=-\frac{1}{2}{{\Delta}}_{y}+\frac{1}{2}{k}_{2}^{0}\vert y{\vert }^{2}+\frac{1}{4}{k}_{4}^{0}\vert y{\vert }^{4},\quad {k}_{4}^{0} > 0,\end{equation*}$$

guarantees that H is essentially self-adjoint. Indeed, Young's inequality for products yields

$$\begin{equation*}W(x,y)\geqslant -\frac{1}{4}\left(\frac{1}{{c}_{0}}\vert \eta {\vert }^{2}\vert x{\vert }^{2}+{c}_{0}\vert y{\vert }^{4}\right),\quad \forall \enspace {c}_{0} > 0,\end{equation*}$$

so choosing ${c}_{0}={k}_{4}^{0}$, we have, for the total potential V(x, y) = V_s(x) + V_b(y) + W(x, y),

$$\begin{equation*}V(x,y)\geqslant \vert {V}_{\mathrm{s}}(x)\vert -\frac{\vert \eta {\vert }^{2}}{8{k}_{4}^{0}}\vert x{\vert }^{2}\geqslant -{C}_{1}\vert x{\vert }^{2}-{C}_{2},\end{equation*}$$

for some constants C₁, C₂ ⩾ 0. Hence, the Faris–Lavine theorem implies that H_x , H_y and H are essentially self-adjoint. In the following, we will therefore also provide slight extensions of our estimates to accommodate this specific, but interesting type of coupling (see remark 8).

We consider the uncoupled system of equations

$$\begin{equation}\begin{cases}\mathrm{i}{\partial }_{t}{\varphi }^{x}={H}_{x}{\varphi }^{x}+W(x,0){\varphi }^{x};\quad {\varphi }_{\vert t=0}^{x}={\varphi }_{0}^{x},\quad \\ \mathrm{i}{\partial }_{t}{\varphi }^{y}={H}_{y}{\varphi }^{y}+W(0,y){\varphi }^{y};\quad {\varphi }_{\vert t=0}^{y}={\varphi }_{0}^{y}.\quad \end{cases}\end{equation} \tag{ 5 }$$

In view of assumption 1, these equations have unique solutions ${\varphi }^{x}\in C(\mathbb{R};{L}_{x}^{2})$, ${\varphi }^{y}\in C(\mathbb{R};{L}_{y}^{2})$, and higher regularity is propagated, ${\varphi }^{x}\in C(\mathbb{R};{{\Sigma}}_{x}^{k})$, ${\varphi }^{y}\in C(\mathbb{R};{{\Sigma}}_{y}^{k})$, for all $k\in \mathbb{N}$, where we recall that Σ^k has been defined in (4). The plain product solves

$$\begin{equation*}\mathrm{i}{\partial }_{t}({\varphi }^{x}{\varphi }^{y})=H({\varphi }^{x}{\varphi }^{y})+\left(-W(x,y)+W(x,0)+W(0,y)\right)({\varphi }^{x}{\varphi }^{y}).\end{equation*}$$

This is not the right approximation, since the residual term is not small: even if W varies very little in y, then W(x, y) − W(x, 0) − W(0, y) ≈ W(x, 0) − W(x, 0) − W(0, 0) = −W(0, 0). This term is removed by considering instead

$$\begin{equation*}{\psi }_{\mathrm{app}}(t,x,y)=\mathrm{exp}(\mathrm{i}tW(0,0)){\varphi }^{x}(t,x){\varphi }^{y}(t,y).\end{equation*}$$

It satisfies the equation

$$\begin{equation*}\mathrm{i}{\partial }_{t}{\psi }_{\mathrm{app}}=H{\psi }_{\mathrm{app}}-\underset{=:{{\Sigma}}_{\psi }}{\underbrace{\left(W(x,y)-{W}_{\mathrm{app}}(x,y)\right){\psi }_{\mathrm{app}}}}.\end{equation*}$$

with approximative potential W_app(x, y) = W(x, 0) + W(0, y) − W(0, 0). The last term Σ_ψ controls the error ψ − ψ_app, as we will see more precisely below. Saying that the coupling potential W is flat in y means that ∇_y W is small, and we write

$$\begin{equation*}W(x,y)-{W}_{\mathrm{app}}(x,y)=\underset{\approx y\cdot {\nabla }_{y}W(x,0)}{\underbrace{\left(W(x,y)-W(x,0)\right)}}-\underset{\approx y\cdot {\nabla }_{y}W(0,0)}{\underbrace{\left(W(0,y)-W(0,0)\right)}}.\end{equation*}$$

This suggests that partial flatness of W implies smallness of the approximation error.

Remark 4. For choosing another collocation point than the origin, one might use the matrix

$$\begin{equation*}M(x,y)={\partial }_{x}{\partial }_{y}W(x,y)={\left({\partial }_{{x}_{j}}{\partial }_{{y}_{k}}W(x,y)\right)}_{1\leqslant j\leqslant n,\enspace 1\leqslant k\leqslant d}.\end{equation*}$$

The Taylor expansion

$$\begin{align}\quad W(x,y)-{W}_{\mathrm{app}}(x,y)& =W(x,y)-W(x,0)-W(0,y)+W(0,0)\\ \quad & =y\cdot {\int }_{0}^{1}\left({\partial }_{y}W(x,\eta y)-{\partial }_{y}W(0,\eta y)\right)\mathrm{d}\eta \\ & ={\int }_{0}^{1}{\int }_{0}^{1}y\cdot {\partial }_{x}{\partial }_{y}W(\theta x,\eta y)\enspace x\enspace \mathrm{d}\eta \enspace \mathrm{d}\theta .\end{align} \tag{ 6 }$$

implies for (x, y) ≈ (x₀, y₀) that

$$\begin{equation*}W(x,y)-W(x,{y}_{0})-W({x}_{0},y)+W({x}_{0},{y}_{0})\approx (x\cdot M({x}_{0},{y}_{0})y),\end{equation*}$$

which corresponds to the standard normal mode expansion. Hence, choosing (x₀, y₀) such that the maximal singular value of M(x₀, y₀) is minimized, we minimize the error of the brute-force approach.

Instead of pointwise evaluations of the coupling potential, we might also use partial averages for an approximation. We consider

$$\begin{equation}\begin{cases}\mathrm{i}{\partial }_{t}{\phi }^{x}={H}_{x}{\phi }^{x}+{\left\langle W\right\rangle }_{y}(t)\enspace {\phi }^{x};\quad {\phi }_{\vert t=0}^{x}={\varphi }_{0}^{x},\quad \\ \mathrm{i}{\partial }_{t}{\phi }^{y}={H}_{y}{\phi }^{y}+{\left\langle W\right\rangle }_{x}(t)\enspace {\phi }^{y};\quad {\phi }_{\vert t=0}^{y}={\varphi }_{0}^{y},\quad \end{cases}\end{equation} \tag{ 7 }$$

where we have denoted

$$\begin{equation}\quad {\left\langle W\right\rangle }_{y}={\left\langle W\right\rangle }_{y}(t,x)=\frac{\int W(x,y)\vert {\phi }^{y}(t,y){\vert }^{2}\mathrm{d}y}{\int \vert {\phi }^{y}(t,y){\vert }^{2}\mathrm{d}y}=\frac{\int W(x,y)\vert {\phi }^{y}(t,y){\vert }^{2}\mathrm{d}y}{\int \vert {\varphi }_{0}^{y}(y){\vert }^{2}\mathrm{d}y},\end{equation} \tag{ 8 }$$

$$\begin{equation}\quad {\left\langle W\right\rangle }_{x}={\left\langle W\right\rangle }_{x}(t,y)=\frac{\int W(x,y)\vert {\phi }^{x}(t,x){\vert }^{2}\mathrm{d}x}{\int \vert {\phi }^{x}(t,x){\vert }^{2}\mathrm{d}x}=\frac{\int W(x,y)\vert {\phi }^{x}(t,x){\vert }^{2}\mathrm{d}x}{\int \vert {\varphi }_{0}^{x}(x){\vert }^{2}\mathrm{d}x},\end{equation} \tag{ 9 }$$

where we have used the fact that the L²-norms of ϕ^x and ϕ^y are independent of time, since W is real-valued. Note that (7) is the nonlinear system of equations of the time-dependent Hartree approximation. Contrary to the brute-force approach, L² regularity does not suffice to define partial averages in general. In view of assumption 1, a fixed point argument (very similar to the proof of e.g. [30, lemma 13.10]) shows that this system has a unique solution $({\phi }^{x},{\phi }^{y})\in C(\mathbb{R};{{\Sigma}}_{x}^{1}\times {{\Sigma}}_{y}^{1})$, and higher regularity is propagated, ${\phi }^{x}\in C(\mathbb{R};{{\Sigma}}_{x}^{k})$, ${\phi }^{y}\in C(\mathbb{R};{{\Sigma}}_{y}^{k})$, for all k ⩾ 2. The approximate solution is then

$$\begin{equation*}{\phi }_{\mathrm{app}}(t,x,y)={\phi }^{x}(t,x)\enspace {\phi }^{y}(t,y)\enspace \mathrm{exp}\left(\mathrm{i}{\int }_{0}^{t}\left\langle W\right\rangle \mathrm{d}s\right),\end{equation*}$$

with the phase given by the full average

$$\begin{equation*}\left\langle W\right\rangle =\left\langle W\right\rangle (t)=\frac{\int W(x,y)\vert {\phi }^{x}(t,x){\phi }^{y}(t,y){\vert }^{2}\mathrm{d}x\enspace \mathrm{d}y}{\int \vert {\varphi }_{0}^{x}(x){\varphi }_{0}^{y}(y){\vert }^{2}\mathrm{d}x\enspace \mathrm{d}y}.\end{equation*}$$

It solves the equation

$$\begin{equation*}\mathrm{i}{\partial }_{t}{\phi }_{\mathrm{app}}=H{\phi }_{\mathrm{app}}-{{\Sigma}}_{\phi },\hspace{25.0pt}{{\Sigma}}_{\phi }{:=}\left(W-{\left\langle W\right\rangle }_{x}-{\left\langle W\right\rangle }_{y}+\left\langle W\right\rangle \right){\phi }_{\mathrm{app}}.\end{equation*}$$

Remark 5. The correcting phase $\mathrm{exp}(\mathrm{i}{\int }_{0}^{t}\left\langle W\right\rangle \mathrm{d}s)$ seems to be crucial if we want to compute the wave function. On the other hand, since it is a purely time dependent phase factor, it does not affect the usual quadratic observables. The same applies for the phase exp(itW(0, 0)) of the brute-force approximation.

We begin with an approximation result at the level of L²-norms only. For its proof, see section 4.

Proposition 1. Under assumption 1, we have the following error estimates:

Brute-force approach: for ψ_app(t, x, y) = φ^x (t, x)φ^y (t, y)exp(itW(0, 0)) defined by (5),

$$\begin{equation*}{\Vert}\psi (t)-{\psi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}\leqslant \begin{cases}2{\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\int }_{0}^{t}{\Vert}y{\varphi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s,\quad \\ {\Vert}{\nabla }_{x}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}{\int }_{0}^{t}{\Vert}x{\varphi }^{x}(s){{\Vert}}_{{L}_{x}^{2}}{\Vert}y{\varphi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s.\quad \end{cases}\end{equation*}$$

Mean-field approach: for ${\phi }_{\mathrm{app}}(t,x,y)={\phi }^{x}(t,x){\phi }^{y}(t,y)\mathrm{exp}(\mathrm{i}{\int }_{0}^{t}\left\langle W\right\rangle \mathrm{d}s)$ defined by (7),

$$\begin{equation*}{\Vert}\psi (t)-{\phi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}\leqslant \begin{cases}4{\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\int }_{0}^{t}{\Vert}y{\phi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s,\quad \\ 4{\Vert}{\nabla }_{x}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}{\int }_{0}^{t}{\Vert}x{\phi }^{x}(s){{\Vert}}_{{L}_{x}^{2}}{\Vert}y{\phi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s.\quad \end{cases}\end{equation*}$$

We see that the smallness of ${\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}$ controls the error between the exact and the approximate solution in both approaches.

Example 2. An important class of examples consists of those where W is slowly varying in y: W(x, y) = w(x, ηy) with η ≪ 1 and w bounded, as well as its derivatives. In that case ${\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}=\eta {\Vert}{\nabla }_{y}w{{\Vert}}_{{L}^{\infty }}$.

Example 3. In the case W(x, y) = W₁(x)W₂(y), the averaged potentials satisfy

$$\begin{equation*}{\left\langle W\right\rangle }_{y}(t,x)={W}_{1}(x){\left\langle {W}_{2}\right\rangle }_{y}(t),\hspace{25.0pt}{\left\langle W\right\rangle }_{x}(t,y)={\left\langle {W}_{1}\right\rangle }_{x}(t){W}_{2}(y),\end{equation*}$$

with

$$\begin{align*}\hspace{25.0pt}\left\langle {W}_{2}\right\rangle ={\left\langle {W}_{2}\right\rangle }_{y}(t)=\frac{\int {W}_{2}(y)\vert {\phi }^{y}(t,y){\vert }^{2}\mathrm{d}y}{\int \vert {\phi }^{y}(t,y){\vert }^{2}\mathrm{d}y}=\frac{\int {W}_{2}(y)\vert {\phi }^{y}(t,y){\vert }^{2}\mathrm{d}y}{\int \vert {\varphi }_{0}^{y}(y){\vert }^{2}\mathrm{d}y},\\ \hspace{25.0pt}\left\langle {W}_{1}\right\rangle ={\left\langle {W}_{1}\right\rangle }_{x}(t)=\frac{\int {W}_{1}(x)\vert {\phi }^{x}(t,x){\vert }^{2}\mathrm{d}x}{\int \vert {\phi }^{x}(t,x){\vert }^{2}\mathrm{d}x}=\frac{\int {W}_{1}(x)\vert {\phi }^{x}(t,x){\vert }^{2}\mathrm{d}x}{\int \vert {\varphi }_{0}^{x}(x){\vert }^{2}\mathrm{d}x}.\end{align*}$$

In this special product case, the crucial source term takes the form

$$\begin{equation*}{{\Sigma}}_{\phi }(t)=({W}_{1}-{\left\langle {W}_{1}\right\rangle }_{x}(t))\enspace ({W}_{2}-{\left\langle {W}_{2}\right\rangle }_{y}(t))\enspace {\phi }_{\mathrm{app}}(t),\end{equation*}$$

and proposition 1 can be augmented by the gradient-free estimate

$$\begin{align}& {\Vert}\psi (t)-{\phi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}\leqslant {\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\Vert}{\varphi }_{0}^{y}{{\Vert}}_{{L}_{y}^{2}}{\int }_{0}^{t}\sqrt{\left({\left\langle {W}_{1}^{2}\right\rangle }_{x}(s)-{\left\langle {W}_{1}\right\rangle }_{x}^{2}(s)\right)\left({\left\langle {W}_{2}^{2}\right\rangle }_{y}(s)-{\left\langle {W}_{2}\right\rangle }_{y}^{2}(s)\right)}\mathrm{d}s\end{align} \tag{ 10 }$$

The L^∞-norms, that provide the upper bounds in proposition 1, separate as

$$\begin{align*}{\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}& ={\Vert}{W}_{1}{{\Vert}}_{{L}_{x}^{\infty }}{\Vert}{\nabla }_{y}{W}_{2}{{\Vert}}_{{L}_{y}^{\infty }},\\ {\Vert}{\nabla }_{x}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}& ={\Vert}{\nabla }_{x}{W}_{1}{{\Vert}}_{{L}_{x}^{\infty }}{\Vert}{\nabla }_{y}{W}_{2}{{\Vert}}_{{L}_{y}^{\infty }},\end{align*}$$

and it is ${\Vert}{\nabla }_{y}{W}_{2}{{\Vert}}_{{L}^{\infty }}$ that controls the estimates. Suppose we have W₂(y) = η|y|² with η small: ∇W₂ is not bounded, but we can adapt the proof of proposition 1 to get

$$\begin{equation*}{\Vert}\psi (t)-{\phi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}\leqslant 16\eta {\Vert}{W}_{1}{{\Vert}}_{{L}_{x}^{\infty }}{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\int }_{0}^{t}{\Vert}\vert y{\vert }^{2}{\phi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s,\end{equation*}$$

that is, the extra power of y is transferred to the ϕ^y term.

Remark 6. In the spirit of the last observations of example 3, in terms of $\eta {:=}{\Vert}{\left\langle x\right\rangle }^{-{\sigma }_{x}}{\langle y\rangle }^{-{\sigma }_{y}}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}$, for some σ_x , σ_y ⩾ 0, we get in proposition 1:

$$\begin{equation*}{\Vert}\psi (t)-{\phi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}\leqslant 8\eta {\int }_{0}^{t}{\Vert}{\left\langle x\right\rangle }^{{\sigma }_{x}}{\phi }^{x}(s){{\Vert}}_{{L}_{x}^{2}}{\Vert}\langle {y\rangle }^{{\sigma }_{y}}\vert y\vert {\phi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s.\end{equation*}$$

See appendix B for details of the argument.

Remark 7. If V₁ is confining, V₁(x)  ≳  |x|² for |x| ⩾ R (for instance, V₁(x) = c|x|^2k , c > 0 and k a positive integer, a typical case where V₁ may be super-quadratic while H_x and H remain self-adjoint), then we can estimate ${\Vert}x{\phi }^{x}{{\Vert}}_{{L}_{x}^{2}}$ uniformly in time. If V₁ = 0, or more generally if V₁(x) → 0 as |x| → ∞, we must expect some linear growth in time ${\Vert}x{\phi }^{x}(t){{\Vert}}_{{L}_{x}^{2}}\lesssim \left\langle t\right\rangle$, and the order of magnitude in t is sharp, corresponding to a dispersive phenomenon.

Remark 8. The framework of a cubic system–bath coupling $W(x,y)=\frac{1}{2}\overrightarrow{\eta }\cdot x\vert y{\vert }^{2}$ as described in section 2.2 is recovered by taking σ_x = σ_y = 1 in example 3. In addition, in the presence of a quartic confinement with ${k}_{4}^{0} > 0$, in view of remark 7, we also know that ${\Vert}\langle y\rangle \vert y\vert {\phi }^{y}(t){{\Vert}}_{{L}_{y}^{2}}$ is bounded uniformly in t.

Adding control on the gradients of the functions φ^x (t), φ^y (t) respectively ϕ^x (t), ϕ^y (t), allows also error estimates at the level of the kinetic energy. For a proof, see appendix B.2.

Proposition 2. Under assumption 1, there exists a constant C > 0 depending on second order derivatives of the potentials such that we have the following error estimates:

Brute-force approach: for ψ_app(t, x, y) = φ^x (t, x)φ^y (t, y)exp(itW(0, 0)) defined by (5),

$$\begin{align*}& {\Vert}{\nabla }_{x}\psi (t)-{\nabla }_{x}{\psi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}+{\Vert}x\psi (t)-x{\psi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}\\ & \hspace{25.0pt}\leqslant C{\Vert}{\nabla }_{x}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}\times {\int }_{0}^{t}{\mathrm{e}}^{Cs}{\Vert}y{\varphi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\left({\Vert}x{\varphi }^{x}(s){{\Vert}}_{{L}_{x}^{2}}\right.\\ & \hspace{25.0pt}\quad \left.+{\Vert}{\nabla }_{x}{\varphi }^{x}(s){{\Vert}}_{{L}_{x}^{2}}+{\Vert}\vert x\vert {\nabla }_{x}{\varphi }^{x}(s){{\Vert}}_{{L}_{x}^{2}}\right)\mathrm{d}s,\end{align*}$$

and

$$\begin{align*}& {\Vert}{\nabla }_{y}\psi (t)-{\nabla }_{y}{\psi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}+{\Vert}y\psi (t)-y{\psi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}\\ & \hspace{25.0pt}\leqslant C{\Vert}{\nabla }_{x}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}\times {\int }_{0}^{t}{\mathrm{e}}^{Cs}{\Vert}x{\varphi }^{x}(s){{\Vert}}_{{L}_{x}^{2}}\left({\Vert}y{\varphi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\right.\\ & \hspace{25.0pt}\quad +\left.{\Vert}{\nabla }_{y}{\varphi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}+{\Vert}\vert y\vert {\nabla }_{y}{\varphi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\right)\mathrm{d}s.\end{align*}$$

Mean-field approach: for ${\phi }_{\mathrm{app}}(t,x,y)={\phi }^{x}(t,x){\phi }^{y}(t,y)\mathrm{exp}(\mathrm{i}{\int }_{0}^{t}\left\langle W\right\rangle \mathrm{d}s)$ defined by (7), analogous estimates for ψ(t) − ϕ_app(t) hold.

Remark 9. The strategy used to prove proposition 2 can be iterated to infer error estimates in Sobolev spaces of higher order, ${H}^{k}({\mathbb{R}}^{n+d})$ for k ⩾ 2, provided that we consider momenta of the same order k, which explains the interest in the functional spaces Σ^k . Error estimates in such spaces can also be obtained by first proving that ψ and the approximate solution(s) remain in Σ^k , and then interpolating with the L² error estimate from proposition 1.

For obtaining quadratic estimates, we consider observables such as the energy or the momenta, that is, operators that are differential operators of order at most 2 with bounded smooth coefficients. These differential operators have their domain in ${H}^{2}({\mathbb{R}}^{n+d})$, as the operator H. More generally, we could consider pseudo-differential operators B = op(b) associated with a smooth real-valued function b = b(Z) with $Z=(z,\zeta )\in {\mathbb{R}}^{2(n+d)}$, whose action on functions $f\in \mathcal{S}({\mathbb{R}}^{n+d})$ is given by

$$\begin{equation*}\mathrm{o}\mathrm{p}(b)f(z)={(2\pi )}^{-(n+d)}{\int }_{{\mathbb{R}}^{2(n+d)}}b\left(\frac{z+{z}^{\prime }}{2},\zeta \right)\mathrm{exp}(\mathrm{i}\zeta \cdot (z-{z}^{\prime }))f({z}^{\prime })\mathrm{d}\zeta \enspace \mathrm{d}{z}^{\prime }.\end{equation*}$$

We assume that b satisfies the Hörmander condition

$$\begin{equation}\forall \enspace \alpha ,\beta \in {\mathbb{N}}^{n+d},\enspace \enspace \exists {C}_{\alpha ,\beta } > 0,\quad \vert {\partial }_{z}^{\beta }{\partial }_{\zeta }^{\alpha }b(z,\zeta )\vert \leqslant {C}_{\alpha ,\beta }{\langle \zeta \rangle }^{2-\vert \alpha \vert },\end{equation} \tag{ 11 }$$

that is, b is a symbol of order 2, see e.g. [31, chapter I.2]. We shall also consider observables that depend only on the variable x or the variable y. The following estimates are proven in appendix B.3.

Proposition 3. Under assumption 1, for $b\in {C}^{\infty }({\mathbb{R}}^{n+d})$ satisfying 11 and B = op(b), there exists a constant C_b > 0 such that we have the following error estimates:

Brute-force approach: for ψ_app(t, x, y) = φ^x (t, x)φ^y (t, y)exp(itW(0, 0)), defined by (5), the error

$$\begin{equation*}{e}_{\psi }(t)=\left\langle \psi (t),B\psi (t)\right\rangle -\left\langle {\psi }_{\mathrm{app}}(t),B{\psi }_{\mathrm{app}}(t)\right\rangle \end{equation*}$$

satisfies

$$\begin{equation*}\left\vert {e}_{\psi }(t)\right\vert \leqslant {C}_{\mathrm{b}}\underset{\vert \beta \vert \leqslant {N}_{n+d}}{\mathrm{sup}}{\Vert}{\nabla }^{\beta }M{{\Vert}}_{{L}^{\infty }}{\Vert}{\psi }_{0}{{\Vert}}_{{L}^{2}}\left(\mathcal{N}({\psi }_{\mathrm{app}})+t{\Vert}{\psi }_{0}{{\Vert}}_{{L}^{2}}\right),\end{equation*}$$

where N_n+d > 0 depends on n + d, M(x, y) = ∂_x ∂_y W(x, y), while

$$\begin{align*}\mathcal{N}({\psi }_{\mathrm{app}})& ={\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\int }_{0}^{t}{\Vert}y{\varphi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s+{\Vert}{\varphi }_{0}^{y}{{\Vert}}_{{L}_{y}^{2}}{\int }_{0}^{t}{\Vert}x{\varphi }^{x}(s){{\Vert}}_{{L}_{x}^{2}}\mathrm{d}s\\ & \quad +{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\int }_{0}^{t}{\Vert}\nabla (y{\varphi }^{y}(s)){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s+{\Vert}{\varphi }_{0}^{y}{{\Vert}}_{{L}_{y}^{2}}{\int }_{0}^{t}{\Vert}\nabla (x{\varphi }^{x}(s)){{\Vert}}_{{L}_{x}^{2}}\mathrm{d}s\\ & \quad +{\int }_{0}^{t}{\Vert}x{\varphi }^{x}(s){{\Vert}}_{{L}_{x}^{2}}{\Vert}\nabla {\varphi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s+{\int }_{0}^{t}{\Vert}\nabla {\varphi }^{x}(s){{\Vert}}_{{L}_{x}^{2}}{\Vert}y{\varphi }^{y}(s){{\Vert}}_{{L}_{y}^{2}}\mathrm{d}s.\end{align*}$$

Mean-field approach: for ${\phi }_{\mathrm{app}}(t,x,y)={\phi }^{x}(t,x){\phi }^{y}(t,y)\mathrm{exp}(\mathrm{i}{\int }_{0}^{t}\left\langle W\right\rangle \mathrm{d}s)$ defined by (7), the error ⟨ψ(t), Bψ(t)⟩ − ⟨ϕ_app(t), Bϕ_app(t)⟩ satisfies a similar estimate.

Remark 10. The averaging process involved in the action of an observable on a wave function allows to prove estimates like the one in proposition 3, that are more precise than the standard ones stemming from norm estimates,

$$\begin{equation*}\vert {e}_{\psi }(t)\vert \leqslant {\Vert}\psi (t)-{\psi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}\left({\Vert}B\psi (t){{\Vert}}_{{L}^{2}}+{\Vert}B{\psi }_{\mathrm{app}}(t){{\Vert}}_{{L}^{2}}\right).\end{equation*}$$

Remark 11. We point out that the error is governed by derivatives of second order in W, involving a derivative in the y variable that is supposed to be small. Besides, note that the direct use of an estimate on the wave function itself would have involved H² norms of ψ_app(s), while this estimate only requires H¹ norms. This first improvement is due to the averaging process present in Egorov theorem.

The error estimates of propositions 1–3, do not allow to distinguish between the brute-force single point collocation and the mean-field Hartree approach. However, in computational practice most of the employed methods are of mean-field type. Why? Our previous analysis, that specifically addresses the coupling of quantum systems, does not allow for an answer, and we resort to a more general point of view. Both approximations, the brute-force and the mean-field one, are norm-conserving. However, the mean-field approach is energy-conserving with the same energy as 1. At first sight, this is surprising, since the mean-field Hamiltonian H_mf(t) depends on time. In a more general framework, where the time-dependent Hartree approximation is considered as application of the time-dependent Dirac–Frenkel variational principle on the manifold of product functions, energy conservation is immediate, see [32, section 3.2] or [33, chapter II.1.5].

Lemma 1. Under assumption 1 and considering the mean-field approach ${\phi }_{\mathrm{app}}(t,x,y)={\phi }^{x}(t,x){\phi }^{y}(t,y)\mathrm{exp}(\mathrm{i}{\int }_{0}^{t}\left\langle W\right\rangle (s)\mathrm{d}s)$ defined by (7), we have

$$\begin{equation*}\left\langle {\phi }_{\mathrm{app}}(t),{H}_{\mathrm{mf}}(t){\phi }_{\mathrm{app}}(t)\right\rangle =\left\langle {\psi }_{0},H{\psi }_{0}\right\rangle \quad \enspace \mathrm{f}\mathrm{o}\mathrm{r}\enspace \mathrm{a}\mathrm{l}\mathrm{l}\enspace t\geqslant 0,\end{equation*}$$

where the mean-field Hamiltonian is given by

$$\begin{equation*}{H}_{\mathrm{mf}}(t)={H}_{x}+{H}_{y}+{\left\langle W\right\rangle }_{y}(t)+{\left\langle W\right\rangle }_{x}(t)-\left\langle W\right\rangle (t).\end{equation*}$$

Below in appendix B.5 we give an elementary ad-hoc proof of lemma 1.

Remark 12. In the brute-force case, the approximate Hamiltonian H_bf = H_x + H_y + W(x, 0) + W(0, y) − W(0, 0) is time-independent, and we have

$$\begin{equation*}\left\langle {\psi }_{\mathrm{app}}(t),{H}_{\text{bf}}{\psi }_{\mathrm{app}}(t)\right\rangle =\left\langle {\psi }_{0},{H}_{\text{bf}}{\psi }_{0}\right\rangle \quad \enspace \mathrm{f}\mathrm{o}\mathrm{r}\enspace \mathrm{a}\mathrm{l}\mathrm{l}\enspace t\geqslant 0.\end{equation*}$$

However this conserved value does not correspond to the exact energy of (1), but only to an approximation of it.

We reconsider the system–bath Hamiltonian with cubic coupling of section 2.2, now formulated in physical coordinates (X, Y), that is,

$$\begin{equation*}{H}_{\text{sb}}=-\frac{{\hslash }^{2}}{2{\mu }_{1}}{{\Delta}}_{X}+{V}_{\mathrm{s}}(X)-\frac{{\hslash }^{2}}{2{\mu }_{2}}{{\Delta}}_{Y}+\frac{{\mu }_{2}{\omega }_{2}^{2}}{2}\vert Y{\vert }^{2}+\frac{1}{2}\overrightarrow{\eta }\cdot X\vert Y{\vert }^{2},\end{equation*}$$

where the coordinates X and Y of the system and the bath part are prescaled, resulting in the single mass parameters μ₁, μ₂ for each subsystem and one single harmonic frequency ω₂ for the bath (noting that, alternatively, several harmonic bath frequencies ω_2,j could be introduced, without modifying the conclusions detailed below). The corresponding time-dependent Schrödinger equation reads

$$\begin{equation*}\mathrm{i}\hslash {\partial }_{\tau }{\Psi}(\tau ,X,Y)={H}_{\text{sb}}{\Psi}(\tau ,X,Y).\end{equation*}$$

We perform a local harmonic expansion of the potential V_s(X) around the origin X = 0 and assume that it is possible to determine a dominant frequency ω₁. We then define the natural length scale of the system as

$$\begin{equation*}{a}_{1}=\sqrt{\frac{\hslash }{{\mu }_{1}{\omega }_{1}}}.\end{equation*}$$

Rescaling coordinates as $(x,y)=\frac{1}{a}(X,Y)$, we obtain

$$\begin{equation*}{H}_{\text{sb}}=\hslash {\omega }_{1}\left(-\frac{1}{2}{{\Delta}}_{x}+{V}_{1}(x)-\frac{{\varepsilon }^{2}}{2}{{\Delta}}_{y}+\frac{1}{2}\frac{{\varpi }^{2}}{{\varepsilon }^{2}}\vert y{\vert }^{2}+\frac{1}{2}{\overrightarrow{\eta }\;}^{\prime }\cdot x\vert y{\vert }^{2}\right),\end{equation*}$$

where we have introduced the dimensionless parameters

$$\begin{equation*}\varepsilon =\sqrt{\frac{{\mu }_{1}}{{\mu }_{2}}},\hspace{25.0pt}\varpi =\frac{{\omega }_{2}}{{\omega }_{1}},\end{equation*}$$

and denoted ${V}_{1}(x)=\frac{1}{\hslash {\omega }_{1}}{V}_{\mathrm{s}}({a}_{1}x)$ and ${\overrightarrow{\eta }\;}^{\prime }=\frac{{a}_{1}}{{\mu }_{1}{\omega }_{1}^{2}}\enspace \overrightarrow{\eta }$. The rescaling of the system potential V_s and the coupling vector $\overrightarrow{\eta }$ do not alter their role in the Hamiltonian, whereas the two dimensionless parameters ɛ and ϖ deserve further attention. We now consider the régime where both the mass ratio ɛ between system and bath and the frequency ratio ϖ between bath and system are small, that is, where the system is viewed as 'light' and 'fast' when compared to the 'heavy' and 'slow' bath.

Example 4. For the hydrogen molecule H₂, where the electrons are considered as the quantum subsystem while the interatomic vibration is considered as the classical subsystem, we have μ₁ = m_e and μ₂ = 918.6m_e. Further, the characteristic electronic energy is of the order of ℏω₁ = 1E_h while the first vibrational level is found at ℏω₂ = 0.020 05E_h. Hence the dimensionless parameters are both small, ɛ = 0.032 99 and ϖ = 0.020 05.

Example 5. As a second example, we consider coupled molecular vibrations, exemplified by the H₂ molecule in a 'bath' of rare-gas atoms, here chosen as krypton (Kr) atoms. The H₂ vibration is now considered as a quantum system interacting with weak intermolecular vibrations. The reduced masses are given as μ₁ (H–H) = 0.5u = 911.44m_e (where u refers to atomic mass units), μ₂ (Kr–Kr) = 41.9u = 76.379 × 10³ m_e, and μ₃ (H₂–Kr) = 1.953u = 3560.10m_e. The vibrational quanta associated with these vibrations are ℏω₁(H–H) = 4159.2 cm⁻¹ = 0.0189E_h, ℏω₂(Kr–Kr) = 21.6 cm⁻¹ = 9.82 × 10⁻⁵ E_h, and ℏω₃(H₂–Kr) = 26.8 cm⁻¹ = 1.22 × 10⁻⁴ E_h (see references [34, 35]). The resulting dimensionless mass ratios are given as ${{\epsilon}}_{12}=\sqrt{{\mu }_{1}/{\mu }_{2}}=0.109$ and ${{\epsilon}}_{13}=\sqrt{{\mu }_{1}/{\mu }_{3}}=0.51$, and the corresponding frequency ratios are ϖ₁₂ = ω₂/ω₁ = 0.005 and ϖ₁₃ = ω₃/ω₁ = 0.006. In the case of the H₂–Kr relative motion, note that the frequency ratio ϖ₁₃ is indeed small whereas the mass ratio is ₁₃ ∼ 0.5; this shows that the quantum–classical boundary is less clearly defined than in the first example of coupled electronic-nuclear motions. In such cases, different choices can be made in defining the quantum–classical partitioning.

In an idealized setting, where ɛ is considered as a small positive parameter whose size can be arbitrarily small, we would say that

$$\begin{equation*}\varpi =\mathcal{O}(\varepsilon )\quad \mathrm{a}\mathrm{s}\enspace \varepsilon \to 0,\end{equation*}$$

and view the system–bath Hamiltonian H_sb as an instance of a partially semiclassical operator

$$\begin{equation*}{H}^{\varepsilon }=-\frac{1}{2}{{\Delta}}_{x}+{V}_{1}(x)-\frac{{\varepsilon }^{2}}{2}{{\Delta}}_{y}+{V}_{2}(y)+W(x,y),\end{equation*}$$

whose potentials V₁(x) and V₂(y) are independent of the semiclassical parameter ɛ and satisfy the growth conditions of assumption 1. As emphasized in section 2.2, the cubic coupling potential W(x, y) does not satisfy the subquadratic estimate, but can be controlled by additional moments of the approximate solution. A corresponding rescaling of time, t = ɛω₁ τ, translates the time-dependent Schrödinger equation to its semiclassical counterpart

$$\begin{equation}\mathrm{i}\varepsilon {\partial }_{t}{\psi }^{\varepsilon }(t,x,y)={H}^{\varepsilon }{\psi }^{\varepsilon }(t,x,y),\end{equation} \tag{ 15 }$$

where the physical and the rescaled wave functions are related via

$$\begin{equation*}{\psi }^{\varepsilon }(t,x,y)={a}^{(n+d)/2}\enspace {\Psi}(\tau /(\varepsilon {\omega }_{1}),aX,aY).\end{equation*}$$

Remark 13. Criteria for justifying a semiclassical description are somewhat versatile in the literature. Our scaling analysis shows, that for system–bath Hamiltonians with cubic coupling the obviously small parameter ɛ, that describes a ratio of reduced masses, has to be complemented by an equally small ratio of frequencies ϖ. Otherwise, the standard form of an ɛ-scaled Hamiltonian, as it is typically assumed in the mathematical literature, does not seem appropriate.

As before, the initial data separate scales,

$$\begin{equation}{\psi }^{\varepsilon }(0,x,y)={\varphi }_{0}^{x}(x){\mathtt{\text{g}}}^{\varepsilon }(y),\end{equation} \tag{ 16 }$$

where we now assume that g ^ɛ is a semiclassically scaled wave packet,

$$\begin{equation}{\mathtt{\text{g}}}^{\varepsilon }(y)=\frac{1}{{\varepsilon }^{d/4}}a \left(\frac{y-{q}_{0}}{\sqrt{\varepsilon }}\right)\mathrm{exp}(\mathrm{i}{p}_{0}\cdot (y-{q}_{0})/\varepsilon ),\end{equation} \tag{ 17 }$$

with $({q}_{0},{p}_{0})\in {\mathbb{R}}^{2d}$, a smooth and rapidly decreasing, i.e. $a\in \mathcal{S}({\mathbb{R}}^{d};\mathbb{C})={\cap }_{k\in \mathbb{N}}{{\Sigma}}^{k}$. In the typical case, where the bath is almost structureless (say, near harmonic), the amplitude a could be chosen as a complex Gaussian, but not necessarily. We now seek an approximate solution of the form ${\psi }_{\mathrm{app}}^{\varepsilon }(t,x,y)={\psi }_{1}^{\varepsilon }(t,x){\psi }_{2}^{\varepsilon }(t,y)$, where ${\psi }_{2}^{\varepsilon }$ is a semiclassically scaled wave packet for all time,

$$\begin{equation}{\psi }_{2}^{\varepsilon }(t,y)=\frac{1}{{\varepsilon }^{d/4}}{u}_{2} \left(t,\frac{y-q(t)}{\sqrt{\varepsilon }}\right)\mathrm{exp}(\mathrm{i}p(t)\cdot (y-q(t))/\varepsilon +\mathrm{i}S(t)/\varepsilon ).\end{equation} \tag{ 18 }$$

Here, $(q(t),p(t))\in {\mathbb{R}}^{2d}$, the phase $S(t)\in \mathbb{R}$, and the amplitude ${u}_{2}(t)\in \mathcal{S}({\mathbb{R}}^{d},\mathbb{C})$ must be determined.

Remark 14. We note that our approximation ansatz differs from the adiabatic one, that would write the full Hamiltonian as ${H}^{\varepsilon }=-\frac{{\varepsilon }^{2}}{2}{{\Delta}}_{y}+{H}_{\mathrm{f}}(y)$, where

$$\begin{equation*}{H}_{\mathrm{f}}(y)=-\frac{1}{2}{{\Delta}}_{x}+{V}_{1}(x)+{V}_{2}(y)+W(x,y)\end{equation*}$$

is an operator, that parametrically depends on the 'slow' variable y and acts on the 'fast' degrees of freedom x. From the adiabatic point of view, one would then construct an approximate solution as ${\psi }_{\text{bo}}^{\varepsilon }(t,x,y)={\Phi}(x,y){\psi }_{2}^{\varepsilon }(t,y)$, where Φ(x, y) is an eigenfunction of the operator H_f(y); here, the subscript 'bo' stands for Born–Oppenheimer. The result of corollary 2 emphasizes the difference between these two points of view.

We denote by

$$\begin{equation}{u}_{\mathrm{app}}^{\varepsilon }(t,x,z)={\psi }_{1}^{\varepsilon }(t,x){u}_{2}(t,z)\quad \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\quad z=\frac{y-q(t)}{\sqrt{\varepsilon }}\end{equation} \tag{ 19 }$$

the part of the approximate solution that just contains the amplitude. With this notation,

$$\begin{equation*}{\psi }_{\mathrm{app}}^{\varepsilon }(t,x,y)=\frac{1}{{\varepsilon }^{d/4}}{u}_{\mathrm{app}}^{\varepsilon }(t,x,z){\left.\mathrm{exp}(\mathrm{i}p(t)\cdot z/\sqrt{\varepsilon }+iS(t)/\varepsilon )\right\vert }_{z=\frac{y-q(t)}{\sqrt{\varepsilon }}}.\end{equation*}$$

The analysis developed in the next two sections allows to derive two different approximations, based on ordinary differential equations governing the semiclassical wave packet part, which are justified in section 5.5 (see proposition 4).

Plugging the expression of ${\psi }_{\mathrm{app}}^{\varepsilon }(t,x,y)$ into (15) and writing $y=q(t)+z\sqrt{\varepsilon }$ in combination with the Taylor expansions

$$\begin{align*}{V}_{2}(y)& ={V}_{2}(q(t)+z\sqrt{\varepsilon })={V}_{2}(q(t))+\sqrt{\varepsilon }z\cdot \nabla {V}_{2}(q(t))\\ & \quad +\frac{\varepsilon }{2}\left\langle z,{\nabla }^{2}{V}_{2}(q(t))z\right\rangle +\mathcal{O}({\varepsilon }^{3/2}),\\ W(x,y)& =W(x,q(t)+z\sqrt{\varepsilon })=W(x,q(t))+\sqrt{\varepsilon }z\cdot {\nabla }_{y}W(x,q(t))+\mathcal{O}\left(\varepsilon \right),\end{align*}$$

we find:

$$\begin{align*}& \mathrm{i}\varepsilon {\partial }_{t}{\psi }_{\mathrm{app}}^{\varepsilon }+\frac{1}{2}{{\Delta}}_{x}{\psi }_{\mathrm{app}}^{\varepsilon }+\frac{{\varepsilon }^{2}}{2}{{\Delta}}_{y}{\psi }_{\mathrm{app}}^{\varepsilon }-V(x,y){\psi }_{\mathrm{app}}^{\varepsilon }\\ & \hspace{25.0pt}={\varepsilon }^{-d/4}\mathrm{exp}(\mathrm{i}p(t)\cdot z/\sqrt{\varepsilon }+\mathrm{i}S(t)/\varepsilon )\\ & \quad \hspace{25.0pt}\times \left(\left(p\cdot \dot{q}-\dot{S}-\frac{\vert p{\vert }^{2}}{2}-{V}_{2}(q)-{V}_{1}(x)-W(x,q)\right){u}_{\mathrm{app}}^{\varepsilon }\right.\\ & \quad \hspace{25.0pt}+\sqrt{\varepsilon }\left(-\mathrm{i}\dot{q}\cdot {\nabla }_{z}{u}_{\mathrm{app}}^{\varepsilon }-\dot{p}\cdot z\enspace {u}_{\mathrm{app}}^{\varepsilon }+\mathrm{i}p\cdot {\nabla }_{z}{u}_{\mathrm{app}}^{\varepsilon }-z\cdot \nabla {V}_{2}(q){u}_{\mathrm{app}}^{\varepsilon }\right.\\ & \hspace{25.0pt}\quad \left.-z\cdot {\nabla }_{y}W(x,q){u}_{\mathrm{app}}^{\varepsilon }\right)\\ & \quad \hspace{25.0pt}+\varepsilon \left(\mathrm{i}{\partial }_{t}{u}_{\mathrm{app}}^{\varepsilon }+\frac{1}{2\varepsilon }{{\Delta}}_{x}{u}_{\mathrm{app}}^{\varepsilon }+\frac{1}{2}{{\Delta}}_{z}{u}_{\mathrm{app}}^{\varepsilon }-\frac{1}{2}\left\langle z,{\nabla }^{2}{V}_{2}\left(q\right)z\right\rangle {u}_{\mathrm{app}}^{\varepsilon }\right)\\ & \quad \hspace{25.0pt}\left.+\mathcal{O}({\varepsilon }^{3/2})\right),\end{align*}$$

where the argument of ${u}_{\mathrm{app}}^{\varepsilon }$ and its derivatives are taken in $z=\frac{y-q(t)}{\sqrt{\varepsilon }}$. To cancel the first four terms in the $\sqrt{\varepsilon }$ line, it is natural to require

$$\begin{equation}\dot{q}=p,\hspace{25.0pt}q(0)={q}_{0},\hspace{25.0pt}\dot{p}=-\nabla {V}_{2}(q),\hspace{25.0pt}p(0)={p}_{0}.\end{equation} \tag{ 20 }$$

Now cancelling the first four terms in the first line of the right-hand side yields

$$\begin{equation}\dot{S}(t)=\frac{\vert p(t){\vert }^{2}}{2}-{V}_{2}(q(t)),\quad S(0)=0.\end{equation} \tag{ 21 }$$

In other words, (q(t), p(t)) is the classical trajectory in y, and S(t) is the associated classical action. At this stage, we note that the term $z\cdot {\nabla }_{y}W(x,q){u}_{\mathrm{app}}^{\varepsilon }$ is not compatible with decoupling the variables x and z (or equivalently, x and y). Using that ${\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}$ is assumed to be small, the above computation becomes

$$\begin{align*}& \mathrm{i}\varepsilon {\partial }_{t}{\psi }_{\mathrm{app}}^{\varepsilon }+\frac{1}{2}{{\Delta}}_{x}{\psi }_{\mathrm{app}}^{\varepsilon }+\frac{{\varepsilon }^{2}}{2}{{\Delta}}_{y}{\psi }_{\mathrm{app}}^{\varepsilon }-V(x,y){\psi }_{\mathrm{app}}^{\varepsilon }\\ & \hspace{25.0pt}={\varepsilon }^{-d/4}\mathrm{exp}(\mathrm{i}p(t)\cdot z/\sqrt{\varepsilon }+\mathrm{i}S(t)/\varepsilon )\\ & \hspace{25.0pt}\quad \times \left(\mathrm{i}\varepsilon {\partial }_{t}{u}_{\mathrm{app}}^{\varepsilon }+\frac{1}{2}{{\Delta}}_{x}{u}_{\mathrm{app}}^{\varepsilon }+\frac{\varepsilon }{2}{{\Delta}}_{z}{u}_{\mathrm{app}}^{\varepsilon }\right.\\ & \quad \hspace{25.0pt}\left.-\left(\frac{\varepsilon }{2}\left\langle z,{\nabla }^{2}{V}_{2}\left(q\right)z\right\rangle +{V}_{1}(x)+W(x,q)\right){u}_{\mathrm{app}}^{\varepsilon }\right.\\ & \quad \hspace{25.0pt}\left.+\mathcal{O}\left({\varepsilon }^{3/2}+\sqrt{\varepsilon }\enspace {\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}\right)\right).\end{align*}$$

In view of (19), we set

$$\begin{equation}\mathrm{i}\varepsilon {\partial }_{t}{\psi }_{1}^{\varepsilon }+\frac{1}{2}{{\Delta}}_{x}{\psi }_{1}^{\varepsilon }=\left({V}_{1}(x)+W(x,q)\right){\psi }_{1}^{\varepsilon };\quad {\psi }_{1\vert t=0}^{\varepsilon }={\varphi }_{0}^{x}\end{equation} \tag{ 22 }$$

$$\begin{equation}\mathrm{i}{\partial }_{t}{u}_{2}+\frac{1}{2}{{\Delta}}_{z}{u}_{2}=\frac{1}{2}\left\langle z,{\nabla }^{2}{V}_{2}\left(q\right)z\right\rangle {u}_{2};\quad {u}_{2\vert t=0}=a.\end{equation} \tag{ 23 }$$

Equation (23) is a Schrödinger equation with a time-dependent harmonic potential: it has a unique solution in L² as soon as $a\in {L}^{2}({\mathbb{R}}^{d})$. In addition, since a ∈ Σ^k for all $k\in \mathbb{N}$, ${u}_{2}\in C(\mathbb{R};{{\Sigma}}_{z}^{k})$ for all $k\in \mathbb{N}$. The validity of this approximation is stated in proposition 4 below. Note that if a is a Gaussian state, then u₂ too and its (time-dependent) parameters—width matrix and centre point—can be computed by solving ODEs (see e.g. [30, 33, 36, 37] and references therein).

Following e.g. [38, 39] or [37, section 2], we write

$$\begin{align*}{V}_{2}(y)& ={V}_{2}(q(t)+z\sqrt{\varepsilon })={\langle {V}_{2}\rangle }_{y}(t)+\sqrt{\varepsilon }z\cdot {\langle \nabla {V}_{2}\rangle }_{y}(t)+\frac{\varepsilon }{2}\enspace z\cdot {\langle {\nabla }^{2}{V}_{2}\rangle }_{y}(t)z+{\mathrm{v}}_{1},\\ W(x,y)& =W(x,q(t)+z\sqrt{\varepsilon })=\langle W{(x,\cdot )\rangle }_{y}(t)+{\mathrm{v}}_{2},\end{align*}$$

where the averages are with respect to $\vert {\psi }_{2}^{\varepsilon }(t,y){\vert }^{2}\mathrm{d}y$. For example,

$$\begin{align}{\langle {\nabla }^{2}{V}_{2}\rangle }_{y}(t)=\;\;& \frac{\int {\nabla }^{2}{V}_{2}(y)\vert {\psi }_{2}(t,y){\vert }^{2}\mathrm{d}y}{\int \vert {\psi }_{2}(t,y){\vert }^{2}\mathrm{d}y}\\ \hspace{25.0pt}=\;\;& \frac{1}{{\Vert}a{{\Vert}}_{{L}^{2}({\mathbb{R}}^{d})}^{2}}{\int }_{{\mathbb{R}}^{d}}{\nabla }^{2}{V}_{2}(q(t)+\sqrt{\varepsilon }z)\enspace {\left\vert {u}_{2}(t,z)\right\vert }^{2}\mathrm{d}z,\\ \langle W{(x,\cdot )\rangle }_{y}(t)=\;\;& \frac{1}{{\Vert}a{{\Vert}}_{{L}^{2}({\mathbb{R}}^{d})}^{2}}{\int }_{{\mathbb{R}}^{d}}W(x,q(t)+\sqrt{\varepsilon }z)\enspace {\left\vert {u}_{2}(t,z)\right\vert }^{2}\mathrm{d}z,\end{align} \tag{ 24 }$$

where we anticipate the fact that the ${L}_{y}^{2}$-norm of ${\psi }_{2}^{\varepsilon }(t)$ is independent of time. We almost literally repeat the previous argument and find that

$$\begin{align*}& \mathrm{i}\varepsilon {\partial }_{t}{\psi }_{\mathrm{app}}^{\varepsilon }+\frac{1}{2}{{\Delta}}_{x}{\psi }_{\mathrm{app}}^{\varepsilon }+\frac{{\varepsilon }^{2}}{2}{{\Delta}}_{y}{\psi }_{\mathrm{app}}^{\varepsilon }-V(x,y){\psi }_{\mathrm{app}}^{\varepsilon }\\ & \hspace{25.0pt}={\varepsilon }^{-d/4}\mathrm{exp}(\mathrm{i}p(t)\cdot z/\sqrt{\varepsilon }+\mathrm{i}S(t)/\varepsilon )\\ & \hspace{25.0pt}\quad \times \left(\mathrm{i}\varepsilon {\partial }_{t}{u}_{\mathrm{app}}^{\varepsilon }+\frac{1}{2}{{\Delta}}_{x}{u}_{\mathrm{app}}^{\varepsilon }+\frac{\varepsilon }{2}{{\Delta}}_{z}{u}_{\mathrm{app}}^{\varepsilon }-\left(\frac{\varepsilon }{2}\enspace z\cdot {\left\langle {\nabla }^{2}{V}_{2}\right\rangle }_{y}z\right.\right.\\ & \quad \hspace{25.0pt}\left.+{V}_{1}(x)+\langle W{(x,\cdot )\rangle }_{y}\right){u}_{\mathrm{app}}^{\varepsilon }\\ & \hspace{25.0pt}\quad \left.+\enspace {\tilde{r}}_{1}^{\varepsilon }+{\tilde{r}}_{2}^{\varepsilon }\right){\left.\right\vert }_{z=\frac{y-q(t)}{\sqrt{\varepsilon }}}\end{align*}$$

with ${\tilde{r}}_{j}^{\varepsilon }={\mathrm{v}}_{j}{u}_{\mathrm{app}}^{\varepsilon }$, j = 1, 2, and z is taken as $z=(y-q(t))/\sqrt{\varepsilon }$. The parameters satisfy the equations of motion

$$\begin{equation*}\quad \dot{q}=p,\quad q(0)={q}_{0},\quad \dot{p}=-{\langle \nabla {V}_{2}\rangle }_{y}(t),\quad p(0)={p}_{0},\quad \dot{S}(t)=\frac{\vert p(t){\vert }^{2}}{2}-{\langle {V}_{2}\rangle }_{y}(t),\enspace S(0)=0.\end{equation*}$$

We see that we can now define the approximate solution by:

$$\begin{equation}\mathrm{i}\varepsilon {\partial }_{t}{\psi }_{1}^{\varepsilon }+\frac{1}{2}{{\Delta}}_{x}{\psi }_{1}^{\varepsilon }=\left({V}_{1}(x)+{\left\langle W(x,\cdot )\right\rangle }_{y}(t)\right){\psi }_{1}^{\varepsilon };\quad {\psi }_{1\vert t=0}^{\varepsilon }={\varphi }_{0}^{x},\end{equation} \tag{ 25 }$$

$$\begin{equation*}\mathrm{i}{\partial }_{t}{u}_{2}+\frac{1}{2}{{\Delta}}_{z}{u}_{2}=\frac{1}{2}\enspace z\cdot {\left\langle {\nabla }^{2}{V}_{2}\right\rangle }_{y}(t)\enspace z\enspace {u}_{2};\quad {u}_{2\vert t=0}=a.\end{equation*}$$

Since the matrix ${\langle {\nabla }^{2}{V}_{2}\rangle }_{y}(t)$ is real-valued, we infer that the ${L}_{z}^{2}$-norm of u₂(t) is independent of time, hence ${\Vert}{\psi }_{2}(t){{\Vert}}_{{L}_{y}^{2}}={\Vert}{u}_{2}(t){{\Vert}}_{{L}_{z}^{2}}={\Vert}a{{\Vert}}_{{L}^{2}}$. The equation in u₂ is now nonlinear, and can be solved in Σ¹, since ∇V₂ is at most linear in its argument: ${u}_{2}\in C(\mathbb{R};{{\Sigma}}_{z}^{1})$, and higher Σ^k regularity is propagated. Here again, if a is a Gaussian, then so is u₂ and its width and centre can be computed by solving ODEs (see [30, 33]). Note also that, differently from the previous setting, u₂ is now ɛ-dependent via the quantity ${\langle {\nabla }^{2}{V}_{2}\rangle }_{y}(t)$ (see (24)). However, this dependence is very weak since a Taylor expansion in (24) shows that u₂ is close in any Σ^k  norm from the solution of the equation (23). For this reason, we do not keep memory of this ɛ-dependence and write u₂. By contrast, the ɛ-dependence of ${\psi }_{1}^{\varepsilon }$ is strong since it involves oscillatory features in time.

The main outcome of the approximations can be stated as follows, and is proved in appendix C.1:

Proposition 4. Let ψ^ɛ be the solution to (15) and (16), with g ^ɛ given by (17). Then with ${\psi }_{\mathrm{app}}^{\varepsilon }$ given either like in section 5.3 or like in section 5.4, there exist constants K₀, K₁ independent of ɛ such that for all t ⩾ 0,

$$\begin{equation*}{\Vert}{\psi }^{\varepsilon }(t)-{\psi }_{\mathrm{app}}^{\varepsilon }(t){{\Vert}}_{{L}^{2}}\leqslant {K}_{0}\left(\sqrt{\varepsilon }+\frac{{\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}}{\sqrt{\varepsilon }}\right){\mathrm{e}}^{{K}_{1}t}.\end{equation*}$$

Corollary 1. Assume $\eta {:=}{\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}\ll \sqrt{\varepsilon }$, then for all T > 0,

$$\begin{equation*}\underset{t\in [0,T]}{\mathrm{sup}}{\Vert}{\psi }^{\varepsilon }(t)-{\psi }_{\mathrm{app}}^{\varepsilon }(t){{\Vert}}_{{L}^{2}}=\mathcal{O}\left(\sqrt{\varepsilon }+\frac{\eta }{\sqrt{\varepsilon }}\right).\end{equation*}$$

Remark 15. Using the same techniques as in appendix B.2, one can prove estimates on higher regularity norms, using ɛ-derivatives in y and standard ones in x. For example, if a ∈ Σ⁴, then there exists K₀, K₁ independent of ɛ such that

$$\begin{align*}& {\Vert}\varepsilon {\nabla }_{y}{\psi }^{\varepsilon }(t)-\varepsilon {\nabla }_{y}{\psi }_{\mathrm{app}}^{\varepsilon }(t){{\Vert}}_{{L}^{2}}+{\Vert}y{\psi }^{\varepsilon }(t)-y{\psi }_{\mathrm{app}}^{\varepsilon }(t){{\Vert}}_{{L}^{2}}\\ & \hspace{25.0pt}\leqslant {K}_{0}\left(\sqrt{\varepsilon }{\int }_{0}^{t}{\mathrm{e}}^{{K}_{1}s}{\Vert}{u}_{2}(s){{\Vert}}_{{{\Sigma}}^{4}}\mathrm{d}s+\frac{{\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}}{\sqrt{\varepsilon }}{\int }_{0}^{t}{\mathrm{e}}^{{K}_{1}s}{\Vert}{u}_{2}(s){{\Vert}}_{{{\Sigma}}^{2}}\mathrm{d}s\right).\end{align*}$$

We refer to [40] (see also [30, chapter 12]) for more detailed computations.

Note that, in both approximations, the evolution of u₂ corresponds to the standard quadratic approximation. In particular, if a is Gaussian, then u₂ is Gaussian at all time, and solving the equation in u₂ amounts to solving ordinary differential equations. However, the equation (22) solved by ψ₁(t) is still quantum, such that a reduction of the total space dimension of the quantum system has been made from n + d to n.

Let us now discuss the approximation of observables that we choose as acting only in the variable y. Due to the presence of the small parameter ɛ, we choose semiclassical observables and associate with $b\in {C}_{c}^{\infty }({\mathbb{R}}^{2d})$ (b smooth and compactly supported) the operator op_ɛ (b) whose action on functions $f\in \mathcal{S}({\mathbb{R}}^{d})$ is given by

$$\begin{equation*}{\mathrm{o}\mathrm{p}}_{\varepsilon }(b)f(y)={(2\pi \varepsilon )}^{-d}{\int }_{{\mathbb{R}}^{2d}}b\left(\frac{y+{y}^{\prime }}{2},\xi \right)\mathrm{exp}(\mathrm{i}\xi \cdot (y-{y}^{\prime })/\varepsilon )f({y}^{\prime })\enspace \mathrm{d}\xi \enspace \mathrm{d}{y}^{\prime }.\end{equation*}$$

As before, the error estimate is better for quadratic observables than for the wave functions. More specifically, the following result, that is proved in appendix C.2, improves the error estimate from proposition 4 by a factor $\sqrt{\varepsilon }$.

Proposition 5. Let ψ^ɛ be the solution to (15) and (16), with g ^ɛ given by (17). Then with $b\in {C}_{c}^{\infty }({\mathbb{R}}^{2d})$ and ${\psi }_{\mathrm{app}}^{\varepsilon }$ given either like in section 5.3 or in section 5.4, there exist a constant K independent of ɛ such that for all t ⩾ 0,

$$\begin{equation*}\left\vert \langle {\psi }^{\varepsilon }(t),{\mathrm{o}\mathrm{p}}_{\varepsilon }(b){\psi }^{\varepsilon }(t)\rangle -\langle {\psi }_{\mathrm{app}}^{\varepsilon }(t),{\mathrm{o}\mathrm{p}}_{\varepsilon }(b){\psi }_{\mathrm{app}}^{\varepsilon }(t)\rangle \right\vert \leqslant K\enspace t\left(\varepsilon +{\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}\right).\end{equation*}$$

Remark 16. Of course, we could have considered a mixed setting consisting of pseudodifferential operators as in section 3.4 in the variable x, and semiclassical as above in the variable y. One would then obtain estimates mixing those of propositions 3 and 5.

Proof. To prove the estimates of example 3, recall that we have denoted

$$\begin{equation*}\eta ={\Vert}{\left\langle x\right\rangle }^{-{\sigma }_{x}}\langle {y\rangle }^{-{\sigma }_{y}}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}< \infty .\end{equation*}$$

The fundamental theorem of calculus yields

$$\begin{equation*}W(x,y)-W(x,{y}^{\prime })=(y-{y}^{\prime })\cdot {\int }_{0}^{1}{\nabla }_{y}W\left({y}^{\prime }+\theta (y-{y}^{\prime })\right)\mathrm{d}\theta ,\end{equation*}$$

so we have

$$\begin{equation*}\vert W(x,y)-W(x,{y}^{\prime })\vert \leqslant \vert y-{y}^{\prime }\vert {\left\langle x\right\rangle }^{{\sigma }_{x}}\eta {\int }_{0}^{1}{\left\langle {y}^{\prime }+\theta (y-{y}^{\prime })\right\rangle }^{{\sigma }_{y}}\mathrm{d}\theta \end{equation*}$$

$$\begin{equation*}\leqslant \vert y-{y}^{\prime }\vert {\left\langle x\right\rangle }^{{\sigma }_{x}}\enspace \mathrm{max}{\left(\left\langle y\right\rangle ,\left\langle {y}^{\prime }\right\rangle \right)}^{{\sigma }_{y}}\eta ,\end{equation*}$$

and we replace the pointwise estimate of δW with

$$\begin{equation*}\vert \delta W(x,{x}^{\prime },y,{y}^{\prime })\vert \leqslant \vert y-{y}^{\prime }\vert \left({\left\langle x\right\rangle }^{{\sigma }_{x}}+{\left\langle {x}^{\prime }\right\rangle }^{{\sigma }_{x}}\right)\mathrm{max}{\left(\left\langle y\right\rangle ,\left\langle {y}^{\prime }\right\rangle \right)}^{{\sigma }_{y}}\eta .\end{equation*}$$

The estimate on Σ_ϕ becomes

$$\begin{align*}{\Vert}{{\Sigma}}_{\phi }{{\Vert}}_{{L}^{2}}^{2}& \leqslant {\eta }^{2}{\int }_{{\mathbb{R}}^{n}}{\left({\int }_{{\mathbb{R}}^{n}}\left({\left\langle x\right\rangle }^{{\sigma }_{x}}+{\left\langle {x}^{\prime }\right\rangle }^{{\sigma }_{x}}\right)\vert {\phi }^{x}(t,{x}^{\prime }){\vert }^{2}\mathrm{d}{x}^{\prime }\right)}^{2}\vert {\phi }^{x}(t,x){\vert }^{2}\mathrm{d}x/{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}^{4}\\ & \quad \times \underset{=:{{\epsilon}}_{y}(t)}{\underbrace{{\int }_{{\mathbb{R}}^{d}}{\left({\int }_{{\mathbb{R}}^{d}}\vert y-{y}^{\prime }\vert \mathrm{max}\left({\left\langle y\right\rangle }^{{\sigma }_{y}},{\left\langle {y}^{\prime }\right\rangle }^{{\sigma }_{y}}\right)\vert {\phi }^{y}(t,{y}^{\prime }){\vert }^{2}\mathrm{d}{y}^{\prime }\right)}^{2}\vert {\phi }^{y}(t,y){\vert }^{2}\mathrm{d}y/{\Vert}{\varphi }_{0}^{y}{{\Vert}}_{{L}_{y}^{2}}^{4}}},\end{align*}$$

and we conclude by resuming the same estimates as above:

$$\begin{align*}& {\int }_{{\mathbb{R}}^{n}}{\left({\int }_{{\mathbb{R}}^{n}}\left({\left\langle x\right\rangle }^{{\sigma }_{x}}+{\left\langle {x}^{\prime }\right\rangle }^{{\sigma }_{x}}\right)\vert {\phi }^{x}(t,{x}^{\prime }){\vert }^{2}\mathrm{d}{x}^{\prime }\right)}^{2}\vert {\phi }^{x}(t,x){\vert }^{2}\frac{\mathrm{d}x}{{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}^{4}}\leqslant 4{\Vert}{\left\langle x\right\rangle }^{{\sigma }_{x}}{\phi }^{x}(t){{\Vert}}_{{L}_{x}^{2}}^{2},\end{align*}$$

and, in view of the inequality

$$\begin{align*}& \vert y-{y}^{\prime }\vert \mathrm{max}\left({\left\langle y\right\rangle }^{{\sigma }_{y}},{\left\langle {y}^{\prime }\right\rangle }^{{\sigma }_{y}}\right)\leqslant 2\enspace \mathrm{max}\left({\left\langle y\right\rangle }^{{\sigma }_{y}}\vert y\vert ,{\left\langle {y}^{\prime }\right\rangle }^{{\sigma }_{y}}\vert {y}^{\prime }\vert \right)\leqslant 2\left({\left\langle y\right\rangle }^{{\sigma }_{y}}\vert y\vert +{\left\langle {y}^{\prime }\right\rangle }^{{\sigma }_{y}}\vert {y}^{\prime }\vert \right),\end{align*}$$

we find

$$\begin{equation*}{{\epsilon}}_{y}(t)\leqslant 16{\Vert}{\left\langle y\right\rangle }^{{\sigma }_{y}}\vert y\vert {\phi }^{y}(t){{\Vert}}_{{L}_{y}^{2}}^{2}.\end{equation*}$$

□

To prove error estimates in ${H}^{1}({\mathbb{R}}^{n+d})$, we differentiate (12) in space, and two aspects must be considered: (i) in our framework, the operator ∇_x,y does not commute with H. (ii) We must estimate ∇_x,y Σ_ψ and ∇_x,y Σ_ϕ . Indeed, we compute

$$\begin{equation*}\mathrm{i}{\partial }_{t}{\nabla }_{x}{r}_{\psi }=H{\nabla }_{x}{r}_{\psi }+[{\nabla }_{x},H]{r}_{\psi }+{\nabla }_{x}{{\Sigma}}_{\psi },\end{equation*}$$

and

$$\begin{equation*}[{\nabla }_{x},H]={\nabla }_{x}H-H{\nabla }_{x}={\nabla }_{x}{V}_{1}+{\nabla }_{x}W.\end{equation*}$$

In the typical case where V₁ is harmonic, ∇_x V₁ is linear in x, and so xr_ψ appears as a source term. Note that in the general setting of assumption 1, $\vert {\nabla }_{x}{V}_{1}(x)\vert \lesssim \left\langle x\right\rangle$.

Remark 17. If ∇_x V₁ and ∇_x W are bounded, then lemma 2 yields

$$\begin{equation*}{\Vert}{\nabla }_{x}{r}_{\psi }(t){{\Vert}}_{{L}^{2}}\leqslant {\int }_{0}^{t}\left(C{\Vert}{r}_{\psi }(s){{\Vert}}_{{L}^{2}}+{\Vert}{\nabla }_{x}{{\Sigma}}_{\psi }(s){{\Vert}}_{{L}^{2}}\right)\mathrm{d}s.\end{equation*}$$

The term ${\Vert}{r}_{\psi }(s){{\Vert}}_{{L}^{2}}$ is estimated in proposition 1, and ${\Vert}{\nabla }_{x}{{\Sigma}}_{\psi }(s){{\Vert}}_{{L}^{2}}$ is estimated below.

Multiplying (12) by x, we find similarly

$$\begin{equation*}\mathrm{i}{\partial }_{t}(x{r}_{\psi })=H(x{r}_{\psi })+[x,H]{r}_{\psi }+x{{\Sigma}}_{\psi }=H(x{r}_{\psi })+{\nabla }_{x}{r}_{\psi }+x{{\Sigma}}_{\psi }.\end{equation*}$$

Energy estimates provided by lemma 2 applied to the equation for ∇_x r_ψ and xr_ψ then yield a closed system of estimates:

$$\begin{align*}& {\Vert}{\nabla }_{x}{r}_{\psi }(t){{\Vert}}_{{L}^{2}}+{\Vert}x{r}_{\psi }(t){{\Vert}}_{{L}^{2}}\\ & \hspace{25.0pt}\leqslant {\int }_{0}^{t}\left({\Vert}({\nabla }_{x}{V}_{1}+{\nabla }_{x}W){r}_{\psi }(s){{\Vert}}_{{L}^{2}}+{\Vert}{\nabla }_{x}{r}_{\psi }(s){{\Vert}}_{{L}^{2}}\right)\mathrm{d}s\\ & \hspace{25.0pt}\quad +{\int }_{0}^{t}\left({\Vert}{\nabla }_{x}{{\Sigma}}_{\psi }(s){{\Vert}}_{{L}^{2}}+{\Vert}x{{\Sigma}}_{\psi }(s){{\Vert}}_{{L}^{2}}\right)\mathrm{d}s\\ & \hspace{25.0pt}\leqslant C{\int }_{0}^{t}\left({\Vert}x{r}_{\psi }(s){{\Vert}}_{{L}^{2}}+{\Vert}{\nabla }_{x}{r}_{\psi }(s){{\Vert}}_{{L}^{2}}\right)\mathrm{d}s\\ & \hspace{25.0pt}\quad +{\int }_{0}^{t}\left({\Vert}{\nabla }_{x}{{\Sigma}}_{\psi }(s){{\Vert}}_{{L}^{2}}+{\Vert}x{{\Sigma}}_{\psi }(s){{\Vert}}_{{L}^{2}}\right)\mathrm{d}s,\end{align*}$$

where we have used the estimate |∇_x V₁ + ∇_x W| ⩽ C(1 + |x|), and the uncertainty principle (uncertainty in x, Cauchy–Schwarz in y),

$$\begin{equation*}{\Vert}f{{\Vert}}_{{L}^{2}}^{2}\leqslant \frac{2}{n}{\Vert}{\nabla }_{x}f{{\Vert}}_{{L}^{2}}{\Vert}xf{{\Vert}}_{{L}^{2}}.\end{equation*}$$

The Gronwall lemma then yields

$$\begin{equation*}{\Vert}{\nabla }_{x}{r}_{\psi }(t){{\Vert}}_{{L}^{2}}+{\Vert}x{r}_{\psi }(t){{\Vert}}_{{L}^{2}}\leqslant {\int }_{0}^{t}{\mathrm{e}}^{Cs}\left({\Vert}{\nabla }_{x}{{\Sigma}}_{\psi }(s){{\Vert}}_{{L}^{2}}+{\Vert}x{{\Sigma}}_{\psi }(s){{\Vert}}_{{L}^{2}}\right)\mathrm{d}s,\end{equation*}$$

for some C > 0. We compute

$$\begin{equation*}{\nabla }_{x}{{\Sigma}}_{\psi }=\left({\nabla }_{x}W(x,y)-{\nabla }_{x}W(x,0)\right){\psi }_{\mathrm{app}}+\delta W(x,0,y,0){\nabla }_{x}{\psi }_{\mathrm{app}}.\end{equation*}$$

The first term in controlled by $\vert y\vert {\Vert}{\nabla }_{x}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}\vert {\psi }_{\mathrm{app}}\vert$. The second term is controlled like in section 3.3, by replacing ψ_app with ∇_x ψ_app. We can of course resume the same approach when considering ∇_y r_ψ , and the analogue of the above first term is now controlled by $\vert x\vert {\Vert}{\nabla }_{x}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}\vert {\psi }_{\mathrm{app}}\vert$. Finally, in the case of r_ϕ , computations are similar (we do not keep track of the precise dependence of multiplicative constants here).

Proof. We use lemma 3 for the operators H and the approximate Hamiltonian H_bf,

$$\begin{equation}{H}_{\text{bf}}={H}_{x}+{H}_{y}+W(x,0)+W(0,y)-W(0,0),\end{equation} \tag{ B.1 }$$

to obtain

$$\begin{equation*}\vert {e}_{\psi }(t)\vert \leqslant {\int }_{0}^{t}\vert {\rho }_{\psi }(t,s)\vert \enspace \mathrm{d}s,\end{equation*}$$

where

$$\begin{equation*}{\rho }_{\psi }(t,s)=\left\langle {\psi }_{\mathrm{app}}(s),\left[B(t-s),H-{H}_{\text{bf}}\right]{\psi }_{\mathrm{app}}(s)\right\rangle ,\quad B(\sigma )={\text{e}}^{\text{i}\sigma H}B{\text{e}}^{-\text{i}\sigma H}.\end{equation*}$$

By Egorov theorem, see [41, theorem 11.1], the operator B(σ) is also a pseudodifferential operator, that is, B(σ) = op(b(σ)) for some function b(σ) that satisfies the growth condition (11). We have

$$\begin{equation*}\quad H-{H}_{\text{bf}}=W(x,y)-W(x,0)-W(0,y)+W(0,0)=\delta W(x,0,y,0){=:}\delta W(x,y),\end{equation*}$$

with the notations of appendix B.1. Then, by the direct estimate of lemma 4 below,

$$\begin{equation}{\Vert}[B(\sigma ),\delta W]{\psi }_{\mathrm{app}}(s){{\Vert}}_{{L}^{2}}\leqslant {C}_{b(\sigma )}\left({\Vert}\nabla (\delta W){\psi }_{\mathrm{app}}(s){{\Vert}}_{{H}^{1}}+{C}_{2}(\delta W){\Vert}{\psi }_{\mathrm{app}}(s){{\Vert}}_{{L}^{2}}\right),\end{equation} \tag{ B.2 }$$

where C_b(σ) > 0 depends on derivative bounds for the function b(σ) and

$$\begin{equation*}{C}_{2}(\delta W)=\sum\limits _{2\leqslant \vert \alpha \vert \leqslant {N}_{n+d}}{\Vert}{\partial }^{\alpha }\delta W{{\Vert}}_{{L}^{\infty }}.\end{equation*}$$

We therefore obtain

$$\begin{equation*}\vert {\rho }_{\psi }(t,s)\vert \leqslant \enspace {C}_{b(t-s)}\enspace \left({\Vert}\nabla (\delta W){\psi }_{\mathrm{app}}(s){{\Vert}}_{{H}^{1}}+{C}_{2}(\delta W)\enspace {\Vert}{\psi }_{0}{{\Vert}}_{{L}^{2}}\right)\enspace {\Vert}{\psi }_{0}{{\Vert}}_{{L}^{2}}.\end{equation*}$$

Using the rectangular n × d matrix M(x, y) introduced in remark 4, the gradient of δW(x, y) can be written as

$$\begin{equation*}\nabla (\delta W)(x,y)=\left(\begin{matrix}\hfill {\nabla }_{x}W(x,y)-{\nabla }_{x}W(x,0)\hfill \\ \hfill {\nabla }_{y}W(x,y)-{\nabla }_{y}W(0,y)\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {\int }_{0}^{1}M(x,\eta y)y\enspace \mathrm{d}\eta \hfill \\ \hfill {\int }_{0}^{1}\enspace {}^{t}M(\theta x,y)x\enspace \mathrm{d}\theta \hfill \end{matrix}\right)\end{equation*}$$

We estimate the Sobolev norm by

$$\begin{align*}{\Vert}\nabla (\delta W){\psi }_{\mathrm{app}}(s){{\Vert}}_{{H}^{1}}& \leqslant {\Vert}\nabla M{{\Vert}}_{{L}^{\infty }}\left({\Vert}x{\psi }_{\mathrm{app}}(s){{\Vert}}_{{L}^{2}}+{\Vert}y{\psi }_{\mathrm{app}}(s){{\Vert}}_{{L}^{2}}\right)\\ & \quad +{\Vert}M{{\Vert}}_{{L}^{\infty }}\left({\Vert}\nabla (x{\psi }_{\mathrm{app}}(s)){{\Vert}}_{{L}^{2}}+{\Vert}\nabla (y{\psi }_{\mathrm{app}}(s)){{\Vert}}_{{L}^{2}}\right),\end{align*}$$

so that integration in time provides

$$\begin{align*}\vert {e}_{\psi }(t)\vert & \leqslant {C}_{\mathrm{b}}{\Vert}\nabla M{{\Vert}}_{{L}^{\infty }}{\Vert}{\psi }_{0}{{\Vert}}_{{L}^{2}}{\int }_{0}^{t}\left({\Vert}x{\psi }_{\mathrm{app}}(s){{\Vert}}_{{L}^{2}}+{\Vert}y{\psi }_{\mathrm{app}}(s){{\Vert}}_{{L}^{2}}\right)\mathrm{d}s\\ & \quad +{C}_{\mathrm{b}}{\Vert}M{{\Vert}}_{{L}^{\infty }}{\Vert}{\psi }_{0}{{\Vert}}_{{L}^{2}}{\int }_{0}^{t}\left({\Vert}\nabla (x{\psi }_{\mathrm{app}}(s)){{\Vert}}_{{L}^{2}}+{\Vert}\nabla (y{\psi }_{\mathrm{app}}(s)){{\Vert}}_{{L}^{2}}\right)\mathrm{d}s\\ & \quad +{C}_{\mathrm{b}}\enspace {C}_{2}(\delta W)\enspace t\enspace {\Vert}{\psi }_{0}{{\Vert}}_{{L}^{2}}^{2},\end{align*}$$

where the constant C_b = max_σ∈[0,t] C_b(σ) depends on derivatives of b. In the mean-field case, the approximate Hamiltonian is time-dependent,

$$\begin{equation}{H}_{\mathrm{mf}}(t)={H}_{x}+{H}_{y}+{\left\langle W\right\rangle }_{y}(t)+{\left\langle W\right\rangle }_{x}(t)-\left\langle W\right\rangle (t).\end{equation} \tag{ B.3 }$$

The difference of the Hamiltonians is also a function, which is now time-dependent, H − H_mf(t) = W + ⟨W⟩(t) − ⟨W⟩_x (t) − ⟨W⟩_y (t). However, it is easy to check that a similar estimate can be performed, leading to an analogous conclusion. □

We now explain the commutator estimate used in the previous subsection:

Lemma 4. Let N ⩾ 1 and b = b(z, ζ) be a smooth function on ${\mathbb{R}}^{2N}$ satisfying the Hörmander growth condition (11). Let δW be a smooth function on ${\mathbb{R}}^{N}$ with bounded derivatives. Then, there exist constants C_b > 0 and M_N > 0 such that

$$\begin{equation*}{\Vert}[\mathrm{o}\mathrm{p}(b),\delta W]\psi {{\Vert}}_{{L}^{2}}\leqslant {C}_{\mathrm{b}}\enspace \left({\Vert}\nabla (\delta W)\psi {{\Vert}}_{{H}^{1}}+\sum\limits _{2\leqslant \vert \alpha \vert \leqslant {M}_{N}}{\Vert}{\partial }^{\alpha }(\delta W){{\Vert}}_{\infty }{\Vert}\psi {{\Vert}}_{{L}^{2}}\right)\end{equation*}$$

for all $\psi \in {H}^{1}({\mathbb{R}}^{N})$.

Proof. We explicitly write the commutator as

$$\begin{align*}& [\mathrm{o}\mathrm{p}(b),\delta W]\psi (z)={(2\pi )}^{-N}{\int }_{{\mathbb{R}}^{2N}}b\left(\frac{z+{z}^{\prime }}{2},\zeta \right)\mathrm{exp}(\mathrm{i}\zeta \cdot (z-{z}^{\prime }))\left(\delta W({z}^{\prime })\right.\\ & \hspace{25.0pt}\qquad \left.-\delta W(z)\right)\psi ({z}^{\prime })\enspace \mathrm{d}\zeta \enspace \mathrm{d}{z}^{\prime }.\end{align*}$$

We Taylor expand the function δW(z) around the point z', so that

$$\begin{equation*}\delta W(z)-\delta W({z}^{\prime })=\nabla (\delta W)({z}^{\prime })\cdot (z-{z}^{\prime })+(z-{z}^{\prime })\cdot \delta {R}_{2}(z,{z}^{\prime })(z-{z}^{\prime })\end{equation*}$$

with

$$\begin{equation*}\delta {R}_{2}(z,{z}^{\prime })={\int }_{0}^{1}(1-\vartheta ){\nabla }^{2}(\delta W)({z}^{\prime }+\vartheta (z-{z}^{\prime }))\enspace \mathrm{d}\vartheta .\end{equation*}$$

Corresponding to the above decomposition, we write [op(b), δW]ψ(z) = f₁(z) + f₂(z) and estimate the two summands separately. We observe that (z − z')exp(iζ ⋅ (z − z')) = −i∇_ζ  exp(iζ ⋅ (z − z')) and perform an integration by parts to obtain

$$\begin{align*}& {\int }_{{\mathbb{R}}^{2N}}b\left(\frac{z+{z}^{\prime }}{2},\zeta \right)\mathrm{exp}\left(\right.\mathrm{i}\zeta \cdot (z-{z}^{\prime })\enspace \nabla (\delta W)({z}^{\prime })\cdot (z-{z}^{\prime })\enspace \psi ({z}^{\prime })\enspace \mathrm{d}\zeta \mathrm{d}{z}^{\prime }\\ & \hspace{25.0pt}=\mathrm{i}{\int }_{{\mathbb{R}}^{2N}}\nabla (\delta W)({z}^{\prime })\cdot {\nabla }_{\zeta }b\left(\frac{z+{z}^{\prime }}{2},\zeta \right)\mathrm{exp}(\mathrm{i}\zeta \cdot (z-{z}^{\prime }))\psi ({z}^{\prime })\enspace \mathrm{d}\zeta \mathrm{d}{z}^{\prime }\end{align*}$$

Therefore,

$$\begin{equation*}{\Vert}{f }_{1}{{\Vert}}_{{L}^{2}}\leqslant {C}_{\mathrm{b}}\enspace {\Vert}\nabla (\delta W)\psi {{\Vert}}_{{H}^{1}},\end{equation*}$$

where the constant C_b > 0 depends on derivative bounds of the function b. For the remainder term of the above Taylor approximation we write

$$\begin{align*}& {\int }_{{\mathbb{R}}^{2N}}b\left(\frac{z+{z}^{\prime }}{2},\zeta \right)\mathrm{exp}(\mathrm{i}\zeta \cdot (z-{z}^{\prime }))\enspace (z-{z}^{\prime })\cdot \delta {R}_{2}(z,{z}^{\prime })(z-{z}^{\prime })\enspace \psi ({z}^{\prime })\enspace \mathrm{d}\zeta \enspace \mathrm{d}{z}^{\prime }\\ & \hspace{25.0pt}={\int }_{{\mathbb{R}}^{2N}}\mathrm{t}\mathrm{r}\left(\delta {R}_{2}(z,{z}^{\prime }){\nabla }_{\zeta }^{2}b\left(\frac{z+{z}^{\prime }}{2},\zeta \right)\right)\mathrm{exp}(\mathrm{i}\zeta \cdot (z-{z}^{\prime }))\psi ({z}^{\prime })\enspace \mathrm{d}\zeta \enspace \mathrm{d}{z}^{\prime },\end{align*}$$

and obtain that

$$\begin{equation*}{\Vert}{f }_{2}{{\Vert}}_{{L}^{2}}\leqslant {C}_{\mathrm{b}}^{\prime }\sum\limits _{2\leqslant \vert \alpha \vert \leqslant {M}_{N}}{\Vert}{\partial }^{\alpha }(\delta W){{\Vert}}_{\infty }{\Vert}\psi {{\Vert}}_{{L}^{2}},\end{equation*}$$

where C_b' > 0 depends on derivative bounds of b, and M_N > 0 depends on the dimension N. □

Here we provide an elementary ad hoc proof for energy conservation of the time-dependent Hartree approximation, lemma 1.

Proof. A first observation is that

$$\begin{equation*}\left\langle {\phi }_{\mathrm{app}}(t),{H}_{\mathrm{mf}}(t){\phi }_{\mathrm{app}}(t)\right\rangle =\left\langle {\psi }_{0},{H}_{\mathrm{mf}}(0){\psi }_{0}\right\rangle \quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace \mathrm{a}\mathrm{l}\mathrm{l}\enspace t\geqslant 0.\end{equation*}$$

Indeed,

$$\begin{equation*}\frac{d}{\mathrm{d}t}\left\langle {\phi }_{\mathrm{app}}(t),{H}_{\mathrm{mf}}(t){\phi }_{\mathrm{app}}(t)\right\rangle =\left\langle {\phi }_{\mathrm{app}}(t),{\partial }_{t}{H}_{\mathrm{mf}}(t){\phi }_{\mathrm{app}}(t)\right\rangle \end{equation*}$$

$$\begin{equation*}=\left\langle {\phi }_{\mathrm{app}}(t),{\partial }_{t}{W}_{\mathrm{app}}(t){\phi }_{\mathrm{app}}(t)\right\rangle \end{equation*}$$

with ${W}_{\mathrm{app}}(t)={\langle W\rangle }_{y}(t)+{\langle W\rangle }_{x}(t)-\left\langle W\right\rangle (t).$ We deduce

$$\begin{align*}& \frac{\mathrm{d}}{\mathrm{d}t}\left\langle {\phi }_{\mathrm{app}}(t),{H}_{\mathrm{mf}}(t){\phi }_{\mathrm{app}}(t)\right\rangle \\ & \hspace{25.0pt}=\int W(x,y)\left({\partial }_{t}\vert {\phi }^{x}(t,x){\vert }^{2}\enspace \vert {\phi }^{y}(t,y){\vert }^{2}+\vert {\phi }^{x}(t,x){\vert }^{2}\enspace {\partial }_{t}\vert {\phi }^{y}(t,y){\vert }^{2}\right)\mathrm{d}x\mathrm{d}y\\ & \hspace{25.0pt}\quad -\int W(x,y){\partial }_{t}\left(\vert {\phi }^{x}(t,x){\vert }^{2}\enspace \vert {\phi }^{y}(t,y){\vert }^{2}\right)\mathrm{d}x\mathrm{d}y=0,\end{align*}$$

where we have used the self-adjointness of H_mf(t) and norm-conservation in the multiplicative components. Secondly, since

$$\begin{equation*}\left\langle {\psi }_{0},{W}_{\mathrm{app}}(0){\psi }_{0}\right\rangle =\left\langle {\psi }_{0},\left({\left\langle W\right\rangle }_{x}(0)+{\left\langle W\right\rangle }_{y}(0)-\left\langle W\right\rangle (0)\right){\psi }_{0}\right\rangle \end{equation*}$$

$$\begin{equation*}=2\left\langle W\right\rangle (0)-\left\langle W\right\rangle (0)=\left\langle W\right\rangle (0),\end{equation*}$$

the approximate energy coincides with the actual energy, and we obtain the aimed for result. □

In this section, we prove proposition 4, and comment on the constants K₀, K₁, which may be analyzed more explicitly in some cases.

C.1.1. Approximation by partial Taylor expansion

Proof. Section 5.3 defines an approximate solution of the from

$$\begin{equation*}\quad {\psi }_{\mathrm{app}}^{\varepsilon }(t,x,y)={\varepsilon }^{-d/4}\mathrm{exp}(\mathrm{i}p(t)\cdot (y-q(t))/\varepsilon +\mathrm{i}S(t)/\varepsilon ){u}_{2}\left(t,\frac{y-q(t)}{\sqrt{\varepsilon }}\right){\psi }_{1}^{\varepsilon }(t,x)\end{equation*}$$

with

$$\begin{equation*}\mathrm{i}\varepsilon {\partial }_{t}{\psi }_{1}^{\varepsilon }+\frac{1}{2}{{\Delta}}_{x}{\psi }_{1}^{\varepsilon }=\left({V}_{1}(x)+W(x,q)\right){\psi }_{1}^{\varepsilon };\quad {\psi }_{1\vert t=0}^{\varepsilon }={\varphi }_{0}^{x}\end{equation*}$$

and

$$\begin{equation*}\begin{gathered}\quad \mathrm{i}{\partial }_{t}{u}_{2}+\frac{1}{2}{{\Delta}}_{z}{u}_{2}=\frac{1}{2}\left\langle z,{\nabla }^{2}{V}_{2}\left(q\right)z\right\rangle {u}_{2};\quad {u}_{2\vert t=0}=a,\\ \quad \dot{q}=p,\enspace q(0)={q}_{0},\quad \dot{p}=-\nabla {V}_{2}(q),\enspace p(0)={p}_{0},\enspace \enspace S(t)={\int }_{0}^{t}\left(\frac{\vert p(s){\vert }^{2}}{2}-{V}_{2}(q(s))\right)\mathrm{d}s.\end{gathered}\end{equation*}$$

It solves the equation

$$\begin{align*}& \mathrm{i}\varepsilon {\partial }_{t}{\psi }_{\mathrm{app}}^{\varepsilon }+\frac{1}{2}{{\Delta}}_{x}{\psi }_{\mathrm{app}}^{\varepsilon }+\frac{{\varepsilon }^{2}}{2}{{\Delta}}_{y}{\psi }_{\mathrm{app}}^{\varepsilon }-V{\psi }_{\mathrm{app}}^{\varepsilon }={\varepsilon }^{-d/4}\mathrm{exp}(\mathrm{i}p\cdot z/\sqrt{\varepsilon }+\mathrm{i}S/\varepsilon )\left({r}_{1}^{\varepsilon }+{r}_{2}^{\varepsilon }\right),\end{align*}$$

where the remainder ${r}_{1}^{\varepsilon }$ is due to the Taylor expansion in V₂, and satisfies the pointwise estimate

$$\begin{equation*}\vert {r}_{1}^{\varepsilon }(t,x,z)\vert \leqslant \frac{1}{6}\times {\varepsilon }^{3/2}{\Vert}{\nabla }^{3}{V}_{2}{{\Vert}}_{{L}_{y}^{\infty }}\vert {\psi }_{1}^{\varepsilon }(t,x)\vert \times \vert z{\vert }^{3}\vert {u}_{2}(t,z)\vert \end{equation*}$$

while the remainder ${r}_{2}^{\varepsilon }$ is due to the Taylor expansion in W, and satisfies the pointwise estimate

$$\begin{equation*}\vert {r}_{2}^{\varepsilon }(t,x,z)\vert \leqslant \sqrt{\varepsilon }\enspace {\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}\vert {\psi }_{1}^{\varepsilon }(t,x)\vert \times \vert z\vert \vert {u}_{2}(t,z)\vert .\end{equation*}$$

This implies for the L²-norm,

$$\begin{equation*}{\Vert}{r}_{1}^{\varepsilon }(t){{\Vert}}_{{L}^{2}}\leqslant \frac{{\varepsilon }^{3/2}}{6}{\Vert}{\nabla }^{3}{V}_{2}{{\Vert}}_{{L}_{y}^{\infty }}{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\Vert}\vert z{\vert }^{3}{u}_{2}(t){{\Vert}}_{{L}_{z}^{2}},\end{equation*}$$

$$\begin{equation*}{\Vert}{r}_{2}^{\varepsilon }(t){{\Vert}}_{{L}^{2}}\leqslant \sqrt{\varepsilon }\enspace {\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\Vert}z{u}_{2}(t){{\Vert}}_{{L}_{z}^{2}}.\end{equation*}$$

Lemma 2 then yields, with now h = ɛ,

$$\begin{equation*}{\Vert}{\psi }^{\varepsilon }(t)-{\psi }_{\mathrm{app}}^{\varepsilon }(t){{\Vert}}_{{L}^{2}}\leqslant \frac{\sqrt{\varepsilon }}{6}{\Vert}{\nabla }^{3}{V}_{2}{{\Vert}}_{{L}_{y}^{\infty }}{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\int }_{0}^{t}{\Vert}\vert z{\vert }^{3}{u}_{2}(s){{\Vert}}_{{L}_{z}^{2}}\mathrm{d}s\end{equation*}$$

$$\begin{equation*}\quad +\frac{{\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}}{\sqrt{\varepsilon }}{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\int }_{0}^{t}{\Vert}\vert z\vert {u}_{2}(s){{\Vert}}_{{L}_{z}^{2}}\mathrm{d}s.\end{equation*}$$

According to the signature of ∇² V₂(q(t)), the quantities ${\Vert}\vert z{\vert }^{3}{u}_{2}(s){{\Vert}}_{{L}_{z}^{2}}$ and ${\Vert}\vert z\vert {u}_{2}(s){{\Vert}}_{{L}_{z}^{2}}$ may be bounded uniformly in s ⩾ 0 or not. For instance, they are bounded if ∇² V₂ is uniformly positive definite, or at least uniformly positive definite along the trajectory q. On the other hand, we always have an exponential bound, even if it may not be sharp,

$$\begin{equation*}{\Vert}\vert z{\vert }^{3}{u}_{2}(s){{\Vert}}_{{L}_{z}^{2}}+{\Vert}\vert z\vert {u}_{2}(s){{\Vert}}_{{L}_{z}^{2}}\leqslant {C}_{0}{\mathrm{e}}^{{C}_{1}s},\end{equation*}$$

for some constants C₀, C₁ > 0. This control is sharp in the case where ∇² V₂ is uniformly negative definite. See e.g. [30, lemma 10.4] for a proof of the exponential control, and [30, section 10.5] for a discussion on its optimality. In particular, for bounded time intervals, the (relative) error is small if ${\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}\ll \sqrt{\varepsilon }\ll 1$. □

Remark 18. If ∇_y W is not bounded, e.g. ${\nabla }_{y}W(x,y)=\eta {\left\langle x\right\rangle }^{\gamma }$, then we can replace the previous error estimate with

$$\begin{equation*}{\Vert}{\psi }^{\varepsilon }(t)-{\psi }_{\mathrm{app}}^{\varepsilon }(t){{\Vert}}_{{L}^{2}}\leqslant \frac{\sqrt{\varepsilon }}{6}{\Vert}{\nabla }^{3}{V}_{2}{{\Vert}}_{{L}_{y}^{\infty }}{\Vert}{\varphi }_{0}^{x}{{\Vert}}_{{L}_{x}^{2}}{\int }_{0}^{t}{\Vert}\vert z{\vert }^{3}{u}_{2}(s){{\Vert}}_{{L}_{z}^{2}}\mathrm{d}s\end{equation*}$$

$$\begin{equation*}\quad +\frac{\eta }{\sqrt{\varepsilon }}{\int }_{0}^{t}{\Vert}{\left\langle x\right\rangle }^{\gamma }{\psi }_{1}^{\varepsilon }(s){{\Vert}}_{{L}_{x}^{2}}{\Vert}\vert z\vert {u}_{2}(s){{\Vert}}_{{L}_{z}^{2}}\mathrm{d}s.\end{equation*}$$

In other words, the cause for the unboundedness of ∇_y W is transferred to a weight for ${\psi }_{1}^{\varepsilon }$. Similarly, if ∇_y W is unbounded in y, we may change the weight in the terms ${\Vert}\vert z{\vert }^{k}{u}_{2}{{\Vert}}_{{L}_{z}^{2}}$, after substituting y with $q+z\sqrt{\varepsilon }$.

C.1.2. Approximation by partial averaging

Proof. The semiclassical approximation obtained by partial averaging reads:

$$\begin{equation*}{\psi }_{\mathrm{app}}^{\varepsilon }(t,x,y)={\varepsilon }^{-d/4}\mathrm{exp}(\mathrm{i}p(t)\cdot (y-q(t))/\varepsilon +\mathrm{i}S(t)/\varepsilon ){u}_{2}\left(t,\frac{y-q(t)}{\sqrt{\varepsilon }}\right){\psi }_{1}^{\varepsilon }(t,x)\end{equation*}$$

with

$$\begin{equation*}\mathrm{i}\varepsilon {\partial }_{t}{\psi }_{1}^{\varepsilon }+\frac{1}{2}{{\Delta}}_{x}{\psi }_{1}^{\varepsilon }=\left({V}_{1}(x)+{\left\langle W(x,\cdot )\right\rangle }_{y}(t)\right){\psi }_{1}^{\varepsilon };\quad {\psi }_{1\vert t=0}^{\varepsilon }={\varphi }_{0}^{x}\end{equation*}$$

and

$$\begin{equation*}\begin{gathered}\quad \mathrm{i}{\partial }_{t}{u}_{2}+\frac{1}{2}{{\Delta}}_{z}{u}_{2}=\frac{1}{2}\enspace z\cdot {\left\langle {\nabla }^{2}{V}_{2}\right\rangle }_{y}(t)z\enspace {u}_{2};\quad {u}_{2\vert t=0}=a,\\ \quad \dot{q}=p,\quad q(0)={q}_{0},\quad \dot{p}=-{\langle \nabla {V}_{2}\rangle }_{y}(t),\quad p(0)={p}_{0},\enspace \enspace \dot{S}(t)=\frac{\vert p(t){\vert }^{2}}{2}-{\langle {V}_{2}\rangle }_{y}(t).\end{gathered}\end{equation*}$$

To estimate the size of ${\tilde{r}}_{1}$ and ${\tilde{r}}_{2}$ introduced in section 5.4, we might argue again via Taylor expansion. Indeed,

$$\begin{equation*}{\Vert}a{{\Vert}}_{{L}^{2}}^{2}{\left\langle {V}_{2}\right\rangle }_{y}=\int {V}_{2}(q(t)+\sqrt{\varepsilon }z)\enspace \vert {u}_{2}(t,z){\vert }^{2}\mathrm{d}z\end{equation*}$$

$$\begin{equation*}={V}_{2}(q(t))+\sqrt{\varepsilon }\enspace \nabla {V}_{2}(q(t))\cdot \int z\vert {u}_{2}(t,z){\vert }^{2}\mathrm{d}z+{r}_{3}^{\varepsilon }(t),\end{equation*}$$

where

$$\begin{equation*}\vert {r}_{3}^{\varepsilon }(t)\vert \leqslant \frac{\varepsilon }{2}{\Vert}{\nabla }^{2}{V}_{2}{{\Vert}}_{{L}^{\infty }}{{\Vert}\vert z\vert {u}_{2}(t,z){\Vert}}_{{L}_{z}^{2}}^{2}.\end{equation*}$$

Hence, we have for all averages f = V₂, ∇V₂, ∇² V₂, W(x, ⋅) that

$$\begin{equation*}{\left\langle f\right\rangle }_{y}(t)=f\left(q(t)\right)+\mathcal{O}(\sqrt{\varepsilon }),\end{equation*}$$

where the error constant depends on moments of |u₂|². In particular, if u₂(0) is Gaussian, the odd moments of |u₂(t, z)|² vanish, and the above estimate improves to $\mathcal{O}(\varepsilon )$. Hence, the L²-norm of ${\tilde{r}}_{1}$ is $\mathcal{O}(\sqrt{\varepsilon })$ close to the L²-norm of r₁, and the L²-norm of ${\tilde{r}}_{2}$ is $\mathcal{O}(\eta \sqrt{\varepsilon })$, $\eta ={\Vert}{\nabla }_{y}W{{\Vert}}_{{L}^{\infty }}$, close to the L²-norm of r₂ (with each time an extra $\sqrt{\varepsilon }$ gain in the above mentioned Gaussian case). In particular, the order of magnitude for the difference between exact and approximate solution is the same as in the previous subsection, only multiplicative constants are affected. We emphasize that the constants C₀ and C₁ from the previous subsection are in general delicate to assess. On the other hand, in specific cases (typically when u₂ is Gaussian and ∇² V₂ is known), they can be computed rather explicitly. □

The proof of proposition 5 is discussed in the next two sections.

C.2.1. Approximation by Taylor expansion

Proof. Taylor expansion yields a time-dependent Hamiltonian ${H}_{\mathrm{app}}^{\varepsilon }={H}_{\text{te}}^{\varepsilon }$ with

$$\begin{equation*}{H}_{\text{te}}^{\varepsilon }{:=}-\frac{1}{2}{{\Delta}}_{x}-\frac{{\varepsilon }^{2}}{2}{{\Delta}}_{y}+{V}_{1}(x)+W(x,q)+{V}_{2}(q)+(y-q)\cdot \nabla {V}_{2}(q)\end{equation*}$$

$$\begin{equation*}\quad +\frac{1}{2}\left\langle y-q,{\nabla }^{2}{V}_{2}(q)(y-q)\right\rangle ,\end{equation*}$$

where q = q(t). In particular, the difference ${H}^{\varepsilon }-{H}_{\text{te}}^{\varepsilon }$ is a function,

$$\begin{equation*}{H}^{\varepsilon }-{H}_{\text{te}}^{\varepsilon }=W(x,y)-W(x,q)+{V}_{2}(y)-{V}_{2}(q)-(y-q)\cdot \nabla {V}_{2}(q)\end{equation*}$$

$$\begin{equation*}\quad -\frac{1}{2}\left\langle y-q,{\nabla }^{2}{V}_{2}(q)(y-q)\right\rangle {=:}\delta W(t,x,y).\end{equation*}$$

In view of lemma 3, if B = op_ɛ (b) with $b\in {C}_{c}^{\infty }({\mathbb{R}}^{2d})$, it yields (a posteriori estimate)

$$\begin{equation*}\left\vert \left\langle {\psi }^{\varepsilon }(t),B{\psi }^{\varepsilon }(t)\right\rangle -\left\langle {\psi }_{\mathrm{app}}^{\varepsilon }(t),B{\psi }_{\mathrm{app}}^{\varepsilon }(t)\right\rangle \right\vert \leqslant \frac{1}{\varepsilon }{\int }_{0}^{t}\vert {\rho }^{\varepsilon }(t,x)\vert \mathrm{d}s,\end{equation*}$$

where

$$\begin{equation*}\begin{aligned}{\rho }^{\varepsilon }(t,s)& =\left\langle {\psi }_{\mathrm{app}}^{\varepsilon }(s),[B(t-s),\delta W(s)]{\psi }_{\mathrm{app}}^{\varepsilon }(s)\right\rangle ,\\ \quad B(\sigma )& =\mathrm{exp}(\mathrm{i}\sigma H/\varepsilon )B\mathrm{exp}(-\mathrm{i}\sigma H/\varepsilon ).\end{aligned}\end{equation*}$$

By Egorov theorem [41, theorem 11.1], B(σ) = ɛ  op_ɛ (b(σ)) for a function $b(\sigma )\in {C}_{c}^{\infty }({\mathbb{R}}^{d})$. Therefore, by semiclassical calculus,

$$\begin{equation*}\frac{1}{\mathrm{i}\varepsilon }[B(t-s),\delta W(s)]={\mathrm{o}\mathrm{p}}_{\varepsilon }\left(\left\{b(t-s),\delta W(s)\right\}\right)+{\varepsilon }^{2}{\mathrm{o}\mathrm{p}}_{\varepsilon }\left({r}^{\varepsilon }(s,t)\right),\end{equation*}$$

where ${\Vert}{\mathrm{o}\mathrm{p}}_{\varepsilon }({r}^{\varepsilon }(s,t)){{\Vert}}_{\mathcal{L}({L}^{2})}$ is bounded uniformly in ɛ, whence the estimate of proposition 5. □

C.2.2. Approximation by partial averaging

Proof. The time-dependent Hamiltonian is ${H}_{\mathrm{app}}^{\varepsilon }={H}_{\text{pa}}^{\varepsilon }$ with

$$\begin{align*}{H}_{\text{pa}}^{\varepsilon }& =-\frac{1}{2}{{\Delta}}_{x}-\frac{{\varepsilon }^{2}}{2}{{\Delta}}_{y}+{V}_{1}(x)+\langle W{(x,\cdot )\rangle }_{y}+{\langle {V}_{2}\rangle }_{y}(t)+(y-q)\cdot {\langle \nabla {V}_{2}\rangle }_{y}(t)\\ & \quad +\frac{1}{2}(y-q)\cdot {\langle {\nabla }^{2}{V}_{2}\rangle }_{y}(t)(y-q),\end{align*}$$

where q = q(t). In particular, as in the preceding case, the difference ${H}^{\varepsilon }-{H}_{\text{pa}}^{\varepsilon }$ is a time-dependent function

$$\begin{align*}{H}^{\varepsilon }-{H}_{\text{pa}}^{\varepsilon }& =W(x,y)-\langle W{(x,\cdot )\rangle }_{y}+{V}_{2}(y)-{\langle {V}_{2}\rangle }_{y}(t)-(y-q)\cdot {\langle \nabla {V}_{2}\rangle }_{y}(t)\\ & \quad -\frac{1}{2}(y-q)\cdot {\langle {\nabla }^{2}{V}_{2}\rangle }_{y}(t)(y-q)\\ & {=:}\tilde{\delta W}(t,x,y),\end{align*}$$

and the arguments developed above also apply. □

The evolution equations for the quantum part of the system, equations (22) and (25), can be written as an adiabatic problem:

$$\begin{equation*}\mathrm{i}\varepsilon {\partial }_{t}{\psi }_{1}^{\varepsilon }(t)=\mathfrak{h}(t){\psi }_{1}^{\varepsilon }(t),\enspace \enspace {\psi }_{1}^{\varepsilon }(0)={\varphi }_{0}^{x},\end{equation*}$$

where $\mathfrak{h}(t)$ is one of the time-dependent self-adjoint operators on ${L}^{2}({\mathbb{R}}^{n})$

$$\begin{equation*}\quad {\mathfrak{h}}_{\text{te}}(t)=-\frac{1}{2}{\Delta}+{V}_{1}(x)+W(x,q(t));\text{and}\quad {\mathfrak{h}}_{\text{pa}}(t)=-\frac{1}{2}{\Delta}+{V}_{1}(x)+{\langle W(x,.)\rangle }_{y}(t).\end{equation*}$$

We assume here that $\mathfrak{h}(t)$ has a compact resolvent and thus, that its spectrum consists in a sequence of time-dependent eigenvalues

$$\begin{equation*}{{\Lambda}}_{1}(t)\leqslant {{\Lambda}}_{2}(t)\leqslant \cdots \leqslant {{\Lambda}}_{k}(t)\underset{k\to +\infty }{\to }\infty .\end{equation*}$$

We also assume that some eigenvalue Λ_j (t) is separated from the remainder of the spectrum for all $t\in \mathbb{R}$ and that the initial datum ${\varphi }_{0}^{x}$ is in the eigenspace of $\mathfrak{h}(0)$ for the eigenvalue Λ_j (0):

$$\begin{equation}\mathfrak{h}(0){\varphi }_{0}^{x}={{\Lambda}}_{j}(0){\varphi }_{0}^{x}.\end{equation} \tag{ C.1 }$$

Then adiabatic theory as developed by Kato [42] states that ψ₁(t) stays in the eigenspace of Λ_j (t) on finite time, up to a phase.

Proposition 6 (Kato [42]). Assume we have (C.1) and that Λ_j (0) is a simple eigenvalue of $\mathfrak{h}(0)$ such that there exists δ₀ > 0 for which

$$\begin{equation*}d\left(\left\{{{\Lambda}}_{j}(t)\right\},\mathrm{S}\mathrm{p}(\mathfrak{h}(t)){\backslash}\left\{{{\Lambda}}_{j}(t)\right\}\right)\geqslant {\delta }_{0}.\end{equation*}$$

Denote by ${{\Phi}}_{j}^{x}(t)$ a family of normalized eigenvectors of $\mathfrak{h}(t)$ such that

$$\begin{equation*}{{\Phi}}_{j}^{x}(0)={\varphi }_{0}^{x},\quad \left\langle {{\Phi}}_{j}^{x}(t),{\partial }_{t}{{\Phi}}_{j}^{x}(t)\right\rangle =0.\end{equation*}$$

Then, for all T > 0, there exists a constant C_T > 0 such that

$$\begin{equation*}{{\Vert}{\psi }_{1}^{\varepsilon }(t)-{\mathrm{e}}^{-\frac{\mathrm{i}}{\varepsilon }{\int }_{0}^{t}{{\Lambda}}_{j}(s)\mathrm{d}s}{{\Phi}}_{j}^{x}(t,x){\Vert}}_{{L}_{x}^{2}}\leqslant {C}_{T}\varepsilon .\end{equation*}$$

In contrast to the Born–Oppenheimer point of view recalled in remark 14, we obtain the following time-adiabatic extension for our wave-packet approximation:

Corollary 2. In the setting of propositions 4 and 6, we obtain the following approximate solution

$$\begin{align*}& {\psi }_{\mathrm{app}}^{\varepsilon }(t,x,y)=\mathrm{exp}\left(-\frac{\mathrm{i}}{\varepsilon }{\int }_{0}^{t}{{\Lambda}}_{j}(s)\mathrm{d}s+\frac{\mathrm{i}}{\varepsilon }S(t)+\frac{\mathrm{i}}{\varepsilon }p(t)\cdot (y-q(t))\right){{\Phi}}_{j}^{x}(t,x){u}_{2}\left(t,\frac{y-q(t)}{\sqrt{\varepsilon }}\right).\end{align*}$$