Variational Gaussians for the magnetic Schrödinger equation]
Variational Gaussian approximation for the magnetic Schrödinger equation

S. Burkhard]Selina Burkhard
Institute for Applied and Numerical Mathematics,
Karlsruhe Institute of Technology, 76149 Karlsruhe, Germany

{selina.burkhard, benjamin.doerich,

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID 258734477 – SFB 1173.

B. Dörich]Benjamin Dörich
M. Hochbruck]Marlis Hochbruck
C. Lasser]Caroline Lasser

Department of Mathematics,
Technical University of Munich,
85748 Garching bei München, Germany

In the present paper we consider the semiclassical magnetic

equation, which describes the dynamics of particles under the influence of a magnetic field. The solution of the Schrödinger equation is approximated by a
single Gaussian wave packet via the time-dependent Dirac–Frenkel variational principle.
For the approximation we derive ordinary differential equations of motion for the parameters of the variational solution. Moreover, we prove $L^2$-error bounds and observable error bounds for the approximating Gaussian wave packet.

§ MAIN RESULTS

In the remaining paper we consider (<ref>) and
(<ref>) for  initial data satisfying

\begin{equation}\label{eq:gauss_initialdata}
\sol_0 = \vsol_0\in \Mf\quad \text{ and } \quad \normLtwo{\vsol_0} =1.
\end{equation}

Our first step is to derive equations of motions for the
parameters defining the variational solution $\vsol$. Then we show that in the limit $\scp \to 0$, these equations tend to classical equations of motions.
Moreover, we study geometric properties of the solution and the variational approximation. Finally, we
state error bounds for the solution in the $L^2$-norm and for averages of observables.
Our work generalizes the results in <cit.>
in the sense that we
treat time-dependent, magnetic Hamiltonians. We also generalize the results of <cit.> from the position and momentum operator to sublinear observables in the sense of <Ref>.

For the sake of readability, we postpone the proofs to <Ref>.

§.§ Variational equations of motion

In order to write equations of motion for the parameters of a $u \in \Mf$ we use the short notation

= ,         = ,

v    = (v_j)_j=1^,
= (_j)_j=1^,

J_   = (∂_j _k)_j,k = 1^,
(D^2_,v)_k,l    = ∑_j=1^∂_l∂_k _j v_j.

We start by deriving two equivalent sets of equations for
$0 < \scp \ll 1$. In the following section, we discuss the limit $\scp \to 0$
and show that the two sets lead to the classical
equations of motion for charged particles in a magnetic field
given by  the Hamiltonian function

\begin{align}
\cham(t,\varone,\vartwo)
%&= \frac{1}{2}|\vartwo|^2-\mgPot(t,\varone)\cdot \vartwo +\frac{1}{2}|\mgPot(t,\varone)|^2 + \mltPot(t,\varone)\\
&= \frac{1}{2}|\vartwo|^2-\mgPot(t,\varone)\cdot \vartwo + \mPot(t,\varone),\quad
\end{align}

cf. <cit.>.
The first set of equations of motion reads:

Let $\vsol_0$ satisfy (<ref>) and  be given by its parameters $\pos_0, \mom_0, \wm_0, \pha_0$ defined in (<ref>).

Then, the parameters of the solution $u\in \Mf$ of (<ref>) satisfy

= - ⟨⟩_,

=  /2ł⟨∇ł(J_A^T )̊⟩̊_+ ⟨J_⟩_^T - ⟨∇⟩_,

= - ^2  + ⟨D^2_,⟩_+ ⟨J_⟩_^T + ⟨J_⟩_- ⟨∇^2 ⟩_.

+ /2 ⟨∇^2 ł(J_^T )̊⟩_,

= 1/2 ||^2 + /2 ⟨ł(J_^T )̊⟩_+ /2 ()

- /4((/2 ⟨∇^2 ł(J_^T )̊⟩_+
ł⟨J_⟩̊_^T +⟨J_⟩_+

- ⟨⟩_+ /4(⟨∇^2 ⟩_u),

with initial data $(\pos(0),\mom(0), \wm(0), \pha(0)) = (\pos_0, \mom_0, \wm_0, \pha_0)$.

The proof of <Ref> is given in <Ref>. We observe that in terms
of the classical Hamiltonian function $h$ defined in (<ref>), the equations of motion (<ref>) can be rewritten as

\beginalign+

= -⟨∇_⟩_,

=-⟨∇_ ⟩_- ⟨∇_ ⟩_- ⟨∇_ ⟩_-⟨∇_ ⟩_,

= -⟨⟩_+/4  (B  ^-1) + ^T ⟨∇_ ⟩_.
\endalign+

with the matrix $B\in\C^{d\times d}$ given by

\[
B = \begin{pmatrix}\Id,\wm \end{pmatrix}\langle \nabla^2 \cham\rangle_\vsol
\begin{pmatrix}\Id\\ \wm\end{pmatrix}.
\]

Later on, in <Ref> we extend these findings to the variational dynamics
induced by a general a subquadratic Hamiltonian.

In order to solve (<ref>) numerically, one might adapt the
Boris algorithm originally proposed in <cit.> and recently analyzed
in <cit.>. This algorithm is constructed for
the classical equations of motion for charged particle systems.

Details or an efficient numerical algorithm are ongoing work which will be presented elsewhere.

An alternative approach presented in <cit.> makes use of
a factorization of the width matrix $\wm$ due to Hagedorn.

For the magnetic equation, it leads to differential equations for the factors of $\wm$ instead of (<ref>).

By <cit.>, we can write

\begin{equation} \label{eq:Hag_fac}
\wm = \fcP \fcQ^{-1} \quad \mathrm{and}\quad \im \wm = (\fcQ\fcQ^*)^{-1},
\end{equation}

with complex, invertible, and symplectic matrices $\fcP$ and $\fcQ$. The latter means that for

\begin{equation}\label{eq:PQSym}
Y \coloneqq
\begin{pmatrix}
\re Q\quad&  \im Q\\
\re P \quad&\im P
\end{pmatrix}
\qquad
\text{and}
\qquad
J \coloneqq
\begin{pmatrix}
0& -\Id\\
\Id&0\\
\end{pmatrix}
\in \R^{2d \times 2d} %\quad \mathrm{symplektisch}
\end{equation}

it holds $Y^T JY =J$, or equivalently

\beginalign+
^T- ^T   = 0,

^*- ^*   = 2  Id.
\endalign+

In fact, if $\fcQ$ and $\fcP$ are complex matrices satisfying (<ref>), then $\fcQ$ and $\fcP$ are invertible and the matrix $\wm = \fcP \fcQ^{-1}$ is symmetric with positive definite imaginary part $(\fcQ \fcQ^*)^{-1}$.

This allows us to write the (<ref>) as

\begin{equation}\label{eq:GWP_factorization}
\vsol(\cdot,x) = \exp \Bigl(\frac{\ii}{\scp}\Bigl(\frac{1}{2}(x-\pos)^T \fcP\fcQ^{-1}(x-\pos)+ \mom^T (x-\pos)+ \pha\Bigr)\Bigr)
\end{equation}

and to derive equations of motion for the parameters $(\pos, \mom, Q,P,\pha)$.

Let $\vsol_0$ satisfy (<ref>) and be
given by the parameters $\pos_0, \mom_0$, $\wm_0, \pha_0$. Then the (<ref>) with parameters $(\pos,\mom,\fcQ,\fcP,\pha)$ solving

\beginalign+

= - ł⟨J_⟩̊_,

= ⟨J_⟩_^T + /2 ⟨∇^2 ł(J_)̊⟩_+
⟨∇^2 ⟩_,
\endalign+

and (<ref>), (<ref>), and (<ref>) is the variational solution (<ref>) with initial data

\begin{equation}
(\pos(0),\mom(0), \wm(0), \pha(0)) = (\pos_0, \mom_0, \wm_0, \pha_0).
\end{equation}

If the initial matrices $\fcQ_0$ and $\fcP_0$ are symplectic, then $\fcQ(t)$ and $\fcP(t)$ are symplectic for all times $t \in\R$.

The proof of <Ref> is given in <Ref>.

§.§ Equations of motion in the limit $\scp\to 0$

The classical Hamiltonian function (<ref>) induces the non-autonomous classical Hamiltonian system

\begin{equation}\label{eq:cldiff}
\begin{aligned}
\begin{pmatrix}
\dot{\varone}(t)\\
\dot{\vartwo}(t)
\end{pmatrix}
&= J^{-1}\nabla \cham(t,\varone(t), \vartwo(t))
\\
\begin{pmatrix}
\vartwo(t) -\mgPot(t,\varone(t))\\
-\nabla \mPot(t,\varone(t))
\end{pmatrix}%\\
% \begin{pmatrix}
% \pos(s)\\
% \mom(s)
% \end{pmatrix}
% &=
% \begin{pmatrix}
% \pos_s\\
% \mom_s
% \end{pmatrix}
\end{aligned}
%\dot{q}(t)&= \nabla_p \cham(q,p,t)= p -\mgPot(q,t),		&		 q(s) &= q_s,\\
%\dot{p}(t)&= -\nabla_q \cham(q,p,t)= J_A^T(q,t)p-\frac{1}{2}\nabla |\mgPot(q,t)|^2-\nabla \mltPot(q,t),	&		 p(s) &=p_s,
\end{equation}

with initial data $\varone(s)=\varone_s,\vartwo(s)= \vartwo_s$ and with $J$ defined in (<ref>).

Since $\mgPot(t,\varone)$ and $\mPot(t,\varone)$ are sublinear and subquadratic with respect to $\varone$, the classical Hamiltonian function $\cham(t,\varone,\vartwo)$ is subquadratic with respect to $(\varone,\vartwo)$. This provides a globally Lipschitz continuous
right-hand side for the ordinary differential equation (<ref>), and the Picard–Lindelöf theorem guarantees a unique global solution.

Since $\mgPot(t,\varone)$ and $\mPot(t,\varone)$ are sublinear and subquadratic with respect to $\varone$,
the right-hand side for the ordinary differential equation (<ref>) is locally Lipschitz continuous. There is no , since

\begin{align}
\frac12\partial_t \big(\abs{\varone}^2 + \abs{\vartwo}^2\big) &= \varone^T (\vartwo - \mgPot(\varone)) + \vartwo^T (J_\mgPot^T(\varone)\vartwo -\nabla \mPot(\varone))\\
&\le C\big(1+\abs{\varone}^2 + \abs{\vartwo}^2\big),
\end{align}

where the constant $C>0$ depends on bounds of the potentials.
By Gronwall's lemma, there is no finite time .

This provides the existence of a unique global solution.

The bound in <cit.> states that $\langle \cdot\rangle_\vsol$ tend to point evaluations at $\pos$ as $\scp \to 0$, i.e.,

$\langle \mgPot \rangle_\vsol \to \mgPot(\pos)$.

Hence, we observe that the magnetic equations of motion (<ref>) and (<ref>) tend to classical equations (<ref>) as ${\scp \to 0}$ and (<ref>) to

\begin{align}
\dot{\pha}= \frac{1}{2}|\mom|^2  -  \mPot(\cdot,\pos).
\end{align}

In order to link the set of equations (<ref>) to classical mechanics, we consider the linearization of (<ref>) along the position and momentum parameters $(\pos, \mom)$, i.e.,

\begin{equation} \label{eq:eqmo_lin_classic}
\begin{aligned}
\begin{pmatrix}
\dot{\fcQ}\\
\dot{\fcP}
\end{pmatrix}
&= J^{-1}\nabla^2 \cham(\cdot,\varone,\vartwo)
\begin{pmatrix}
\fcQ\\
\fcP
\end{pmatrix} \\
\fcP - J_\mgPot(\cdot,\varone)\fcQ\\
\left(D^2_{\mgPot(\cdot,\varone), \vartwo}  - \nabla^2 \mPot(\cdot,\varone)\right)\fcQ + J_\mgPot(\cdot,\varone)^T \fcP
\end{pmatrix}.
%&\fcP - J_\mgPot(\pos, \cdot)\fcQ,\\
%&\left(D^2_{\mgPot(\pos, \cdot), \mom}  - \nabla^2 \mPot(\pos,\cdot)\right)\fcQ + J_\mgPot(\pos,\cdot)^T \fcP.
\end{aligned}
\end{equation}

By the same reasoning, we observe that the
equations (<ref>) tend to the linearized equations classical equations (<ref>) as ${\scp \to 0}$.

§.§ Averages

A further remarkable property of is the conservation of several physical quantities. In the following, we recall the definitions of the linear and angular momentum for quantum dynamical systems.

Let $x = (x_1, \ldots, x_N)$, where $x_k \in \R^3,\, k = 1,\ldots, N$
and $d = 3N$, be position variables. We recall the follwoing definition given in <cit.>.

* The quantum mechanical total linear momentum operator is given by

\begin{equation}%\label{eq:lin_Mom}
\linMom := -\ii \scp\sum_{k=1}^N{ \nabla_{x_k}}.
\end{equation}

* The quantum mechanical total angular momentum operator is given by

\begin{equation}%\label{eq:ang_Mom}
\angMom := \sum_{k=1}^N{x_k \times \left( -\ii \scp \nabla_{x_k}\right)} = -\ii \scp \sum_{k=1}^N\begin{pmatrix}
x_{k_3}\partial_{k_1} - x_{k_1}\partial_{k_3}\\
x_{k_1}\partial_{k_2} -x_{k_2}\partial_{k_1}
\end{pmatrix}. %\,\mathrm{Hertischer\, Teil\, von\,}
\end{equation}

Next, we state sufficient conditions on the potentials $\mgPot$ and $\mltPot$, which lead to the conservation of
averages of the observables from <Ref>.

We call a potential $W = (W_j)_{j=1, \ldots, \dim}: (\R^3)^N \to \R^\dim$

* translation invariant, if

\begin{equation}
W_j(x_1, \ldots ,x_N) = W_j(x_1 + r, \ldots , x_N + r),
\end{equation}

for all $r\in \R^3$ and $j= 1, \ldots, \dim$,

* rotation invariant if for all orthogonal matrices $R \in \R^{3\times 3}$ with $\det R =1$ it holds

\begin{equation}
W_j(x_1, \ldots , x_N) = W_j(R x_1 , \ldots , R x_N ),
\end{equation}

where $j=1, \ldots, \dim$.

In the next lemma we provide a representation for the energy and state conservation properties of the momenta.

The following assertions hold.

* We have $\normLtwo{\sol(t)} = \normLtwo{\vsol(t)}  = \normLtwo{\vsol_0}$ for all $t\in\R$.

* If the potentials $\mgPot$ and $\mltPot$ are both , then

\begin{equation}
\langle \Ham \rangle_{\sol(t)} = \langle \Ham \rangle_{\sol_0}\quad\text{ and }\quad\langle \Ham \rangle_{\vsol(t)} = \langle \Ham 	\rangle_{\vsol_0}.
\end{equation}

* For $\varphi = \sol, \vsol$ the energy $\langle \Ham \rangle_\varphi$ is given by

\begin{align}
\langle \Ham(t)\rangle_{\varphi(t)} &= \langle  \Ham(0)\rangle_{\varphi(0)} + \int_0^t{ \bigl\langle  \ii\scp \partial_s \mgPot(s)\cdot \nabla \bigr\rangle_{\varphi(s)} + \bigl\langle\partial_s\mPot(s)\bigr\rangle_{\varphi(s)} }\,\dd s .
\end{align}

* For $\linMom$ and $\angMom$ from <Ref> we have:

(i) If $\mltPot$ and $\mgPot = (\mgPot_j)_{j=1}^\dim$ given in  <Ref> are invariant under translations

\begin{equation}
\langle \linMom \rangle_{\sol(t)} = \langle \linMom \rangle_{\sol_0}
\quad\text{ and }\quad
\langle \linMom \rangle_{\vsol(t)} = \langle \linMom \rangle_{\vsol_0}.
\end{equation}

(ii) If $\mPot$ defined in (<ref>) is invariant under rotations and $\mgPot(\cdot,x)= \alpha(\cdot) x$ for some $\alpha(\cdot) \in \R$, then

\begin{equation}
\langle \angMom \rangle_{\sol(t)} = \langle \angMom \rangle_{\sol_0}
\quad\text{ and }\quad
\langle \angMom \rangle_{\vsol(t)} = \langle \angMom \rangle_{\vsol_0}.
\end{equation}

The proof of <Ref> is given in <Ref>.

§.§ L2-error bound

In this section, we present the approximation property of the Gaussian wave packet the $L^2$-norm.

Since our error bounds depend on parameters characterizing the

Gaussian wave packet in (<ref>), we first consider

the boundedness of these parameters up to a fixed but arbitrary finite time $\etime>0$ specified by -theory.

For all times $\etime >0$, the set of equations (<ref>) is on $[0,\etime]$ independently of $\scp$. Furthermore, the solution parameters are bounded independently of $\scp$, i.e.

\begin{equation}
|\param|\leq c_{\param_0}, \qquad \text{ for all }~ \param \in \{\pos, \mom, \wm, \pha\},
\end{equation}

uniformly on $[0,\etime]$, where $c_{\param_0}$ depends on the parameters of the initial Gaussian $\vsol_0$, on the potentials $\mltPot, \mgPot$, and on $\etime$.

We note that by <Ref> the matrix $\IC$ is real symmetric, positive definite for all times $t$. To formulate the following results, we denote by $\rho >0$ a lower bound on the smallest eigenvalue of $\IC$ on the finite time horizon $[0,\etime]$. For a discussion of relevant time scales on which $\rho$ is sufficiently large compared to $\eps$, called the Ehrenfest time, we refer to <cit.>.

With this,

we can state our approximation result.

Let $\sol, \vsol$ be the solution of (<ref>) and (<ref>), respectively,
and let $\vsol_0$ satisfy (<ref>).

Then the error bound

\begin{equation}
\normLtwo{\sol(t) -\vsol(t)} \leq t c \sqrt{\varepsilon}, \qquad t  \in [0,T],
\end{equation}

holds with a constant $c$ which depends on $\rho$, the bounds on the parameters from
and on the potentials, but is independent of $\varepsilon$ and $t$.

We provide the

details and the proof
of the theorem

in <Ref>.

§.§ Observable error bound

In classical mechanics physical states are described by the position and momentum parameters $\varone, \vartwo \in \R^{\dim}$. Observables are functions depending smoothly on $(\varone, \vartwo ) \in \R^{\dim}\times\R^{\dim}$, see, for example, <cit.>. Classical mechanics can be linked to quantum mechanics via , which asigns a classical observable to a quantum mechanical one using Fourier transformation, cf. <cit.> or <cit.>. Formally, for $ \varphi \in \SchF(\R^\dim)$ and an observable $\clobs$, we define

\begin{align}\label{eq:Weyl_def}
\op(\clobs)\varphi(x) &\coloneqq \frac{1}{(2\pi\scp)^{d}}\int_{\R^{2\dim}}{\clobs\Bigl(\frac{x+\varone}{2},\vartwo\Bigr)\e^{\ii \vartwo \cdot(x-\varone)/\scp}\varphi(\varone)}\,\dd(\varone,\vartwo).
%\text{ for all } \varphi &\in \SchF(\R^\dim)
%&= \frac{1}{(2\pi\eps)^\frac{d}{2}}\int{K_b(x,\varone)\varphi(\varone)}\,\dd \varone,\nonumber
\end{align}

The of the projections to the first or second component of the classical variables are

\begin{align}
\op(\vartwo)\varphi = -\ii\varepsilon \nabla \varphi\quad\text{ and } \quad\op(\varone)\varphi = x\varphi.
\end{align}

Further examples of physically relevant observables stemming from classical symbols are

\begin{align}
\op(|\vartwo|^2)\sol(x) &= -\scp^2\Delta \sol(x)
\end{align}

and, due to $\divergence \mgPot = 0$,

\begin{align}
\op(\mgPot(\varone)\cdot \vartwo)\sol(x)
&= \tfrac12\left(\mgPot(x) \cdot(-\ii\scp \nabla) + (-\ii\scp\nabla)\cdot \mgPot(x)\right) \sol(x)\\
&= \left(\mgPot(x) \cdot(-\ii\scp \nabla)\right) \sol(x),
\end{align}

and, of course,

\[
\op(\cham(t))\sol(x) = \Ham(t)\sol(x)
\]

for the Hamiltonian function (<ref>) and the magnetic Schrödinger operator (<ref>).
An observable $\obs= \op(\clobs)$ defines for an $L^2$-normalised function
$\varphi\in\SchF(\R^\dim)$ an expectation value,

\[
\bigl\langle \varphi| \obs \varphi \bigr\rangle = \int_{\R^\dim} \overline{\varphi(x)} (\obs\varphi)(x) \dd x,
\]

and we investigate how expectation values issued by the variational approximation $\vsol(t)$ differ from the ones of the true solution $\sol(t)$. For an error estimate relying on $L^2$ bounds, we have to restrict ourselves to sublinear classical observables.

The class of sublinear classical symbols is defined as smooth functions $\clobs : \R^{2\dim} \to \R$ such that for $\alpha \in \N^{2\dim}_0$ with $\abs{\alpha} \ge 1$ there exists $C_\alpha>0$

\begin{equation}
\abs{\partial^\alpha \clobs(\varone,\vartwo)} \le C_\alpha
\end{equation}

for all $(\varone,\vartwo)\in\R^\dim\times\R^\dim$.

For the expectation values of classical sublinear observables, we obtain the following error
estimate that generalizes and improves the findings of N. King and T. Ohsawa <cit.>, where asymptotic accuracy of the order $\scp^{3/2}$ has been observed and proved for the variational position and momentum expectation value.

Let $\sol, \vsol$ be the solution of (<ref>) and (<ref>), respectively,
and let $\vsol_0$ satisfy (<ref>).

Moreover, let $\obs= \op(\clobs)$ be an observable stemming from a classical sublinear observable $\clobs$ in the sense of <Ref>

such that $\egor{t}{s}$ is sublinear.

Then we have the error bound

\begin{equation}
\bigl|\bigl\langle \sol(t)| \obs \sol(t) \bigr\rangle - \bigl\langle \vsol(t)| \obs \vsol(t) \bigr\rangle \bigr|\leq
t\,c\,\scp^2, %{\vaiindex{m}}{n+j},
\end{equation}

for all $t\in[0,T]$. The error constant $c$ depends on the parameter bounds of <Ref> for the time-interval $[0,T]$, in particular on
the bounds for the width matrix $\wm$, on the potentials, and on $\clobs$, but is independent of $\scp$ and $t$.

Note that the convergence in the observables is of order $\eps^2$,

while the convergence in the $L^2$-norm presented in <Ref> is of order $\sqrt{\eps}$. This is
an improvement of the results obtained in <cit.>, where $\mathcal O(\sqrt\scp)$ norm accuracy and an $\mathcal O(\scp)$ estimate for the non-magnetic observable error were proved.

The rest of the paper is devoted to the proofs of the equations of motion and the error estimates presented in this section.

§ EQUATIONS OF MOTIONS: PROOF OF <REF>

In this section we derive equations of motion for the parameters $(\pos, \mom, \wm, \pha)$ as well as for the factorization matrices $\fcQ$ and $\fcP$.

To do so, we compute both sides of (<ref>) and compare the coefficients.

In order to use the formula for the orthogonal projection derived in <cit.> for (<ref>), we observe that derivatives $x$ of a turn into
scalar functions of

times $\vsol$. For notational simplicity, we omit the and in the potentials $\mgPot$ and $\mltPot$ we omit the space variable $x$. In particular, we have

\beginalign+
·∇   = -·((x-)+) ,

-^2/2Δ   = (1/2(x-)^T^2(x-) + ^T(x-)+1/2||^2 - /2()  ) ,
\endalign+

and for the time derivative it holds that

\begin{equation}\label{eq:formula-dt-vsol}
\ii\scp \pt \vsol(\cdot,x) =
\Bigl(-\frac{1}{2}(x-\pos)^T\dot{\wm}(x-\pos) + \dot{\pos}^T\wm(x-\pos) - \dot{\mom}^T(x-\pos) + \mom^T \dot{\pos} -\dot{\zeta}\Bigr) \vsol.
\end{equation}

Motivated by the classical magnetic Hamiltonian system (<ref>),

we eliminate one degree of freedom by setting $\dot{\pos} = \mom - \langle \mgPot\rangle_\vsol$, see <cit.>.

Incorporating the above formulas, we compare the coefficients in $x$ on both sides of (<ref>) and arrive at equations of motions of the form

\begin{align}
\dot{\pos} &= \mom - \langle \mgPot\rangle_\vsol,
\\
\dot{\mom} &= \bigl\langle J_\mgPot^T \RC (x-\pos)\bigr\rangle_\vsol + \bigl\langle J_\mgPot \bigr\rangle_\vsol^T \mom  - \bigl\langle \nabla \mPot \bigr\rangle_\vsol,\\
\dot{\wm} &= - \wm^2 + \bigl\langle D^2_{\mgPot,\RC(x-\pos)} \bigr\rangle_\vsol + \bigl\langle  D^2_{\mgPot,\mom}\bigr\rangle_\vsol + \bigl\langle J_\mgPot\bigr\rangle_\vsol^T \wm + \wm \bigl\langle J_\mgPot\bigr\rangle_\vsol  - \bigl\langle \nabla^2 \mPot \bigr\rangle_\vsol,\\
\dot{\pha} &= \frac12 |\mom|^2
\l\langle \mgPot^T \RC(x-\pos) \r\rangle_\vsol
\frac{\ii \eps}{2} \tr (\wm)  \\
&\quad - \frac{\eps}{4}\tr
\bigl(\ICinv
\bigl( \bigl\langle D^2_{\mgPot,\RC(x-\pos)} \bigr\rangle_\vsol
\l\langle J_\mgPot\r\rangle_\vsol^T \RC
\RC \langle J_\mgPot\rangle_\vsol
\langle  D^2_{\mgPot,\mom}\rangle_\vsol\bigr)
\bigr)\\
&\quad  - \langle \mPot\rangle_\vsol + \frac{\eps}{4}\tr\bigl(\ICinv \bigl\langle \nabla^2 \mPot\bigr\rangle_u \bigr).
%- \frac12\l\langle |\mgPot|^2\r\rangle_\vsol + \frac{\eps}{8}\tr\l(\ICinv \l\langle \nabla^2 |\mgPot|^2\r\rangle_\vsol\r)
\end{align}

It remains to extract the additional power of $\scp$ from the terms that contain the difference $x-\pos$.

\begin{equation}
|\vsol(x)|^2 = \exp\Bigl(-\frac{1}{\scp} (x-\pos)^T \IC (x-\pos) - \frac{2}{\scp} \im \zeta\Bigr)
\end{equation}

we obtain the derivative

\begin{align}
%|\vsol(x)|^2 &= \exp\left(-\frac{1}{\scp} (x-\pos)^T \IC (x-\pos) - \frac{2}{\scp} \im \zeta\right),\\ %\label{eq:Betrag u}
\nabla |\vsol(x)|^2 &= -\frac{2}{\scp}\IC(x-\pos)|\vsol(x)|^2,\label{eq:nabla_abs_u} %\\
%	\nabla^2 |\vsol(x)|^2 &=  \frac{1}{\scp}\left(\frac{4}{\scp} \IC (x-\pos) (x-\pos)^T\IC -2 \IC \right)|\vsol(x)|^2,\label{eq:nabla2_abs_u}
\end{align}

and apply integration by parts to obtain

\begin{align}
\bigl\langle \mgPot^T \RC (x-\pos)\bigr\rangle_\vsol &= \bigl\langle \mgPot^T \RC \ICinv \IC (x-\pos)\bigr\rangle _\vsol
\\
&= \int_{\R^\dim} \mgPot^T \RC \ICinv \IC (x-\pos) |\vsol(x)|^2 \dd x
\\
%	&= \frac{\eps}{2}\Bigl\langle \sum_{j,k,l =1}^\dim \partial_j \mgPot_k \RCkl{k}{l} \ICinvkl{lj} \Bigr\rangle\\%_{{\mkern-10mu}\vsol} \\
&= \frac{\eps}{2} \bigl\langle \tr\l(J^T_A \RC \ICinv\r) \bigr\rangle_\vsol .
\end{align}

Similarly, we gain an order of $\scp$ for

\begin{align}
\bigl( \bigl\langle J_\mgPot^T \RC (x-\pos)\bigr\rangle_\vsol \bigr)_{i}
= \bigl(\bigl\langle J_\mgPot^T \RC \ICinv \IC (x-\pos)\bigr\rangle_{\vsol} \bigr)_i
=\frac{\eps}{2} \bigl\langle\partial_i \tr\l(J^T_A \RC \ICinv\r) \bigr\rangle_\vsol,
%	&= \int J_\mgPot^T \RC \ICinv \IC (x-\pos) |\vsol(x)|^2 \dd x\\
%	&= \frac{\eps}{2}\l\langle \sum_{j,k,l =1}^\dim \partial_i \partial_j \mgPot_k \RCkl{k}{l} \ICinvkl{lj} \r\rangle\\%_{{\mkern-10mu}\vsol} \\
%	&= \frac{\eps}{2}\partial_i \tr\l(J_A \RC \ICinv\r).
\end{align}

as well as for

\begin{align}
\bigl(\bigl\langle D^2_{A,\RC(x-q)} \bigr\rangle_\vsol\bigr)_{ij}
%	&=  \bigl(\bigl\langle D^2_{A,\RC \IC \ICinv(x-q)} \bigr\rangle_\vsol\bigr)_{ij} \\
%	&= \Bigl\langle \sum_{k,l,m =1}^\dim \partial_i \partial_j \mgPot_k \RCkl{k l}\ICinvkl{lm} \sum_{n=1}^d \ICkl{m n}(x_n-q_n) \Bigr\rangle_\vsol\\
&= \frac{\eps}{2} \Bigl\langle \sum_{k,l,m =1}^d \partial_m\partial_i \partial_j \mgPot_k \RCkl{k l}\ICinvkl{lm}  \Bigr\rangle_\vsol.
\end{align}

By the identity

\begin{align}
\partial_{ij} \tr\l(J_\mgPot^T \RC \ICinv\r) & = \sum_{k,m,l=1}^\dim \partial_{ij} \partial_m \mgPot_k\RCkl{k l} \ICinvkl{lm},
\end{align}

we conclude the equations of motion stated in (<ref>).

We now turn to the equations of motion for the Hagedorn factorization (<ref>).
The idea is to show that the product $\fcP \fcQ^{-1}$ solves the same differential equation as $\wm$ and conclude with the uniqueness of the variational solution $\vsol$.

We employ the differential identity

\begin{equation}
\pt(\fcQ^{-1} ) = -\fcQ^{-1}\pt{\fcQ}\fcQ^{-1} ,
%		 \quad \text{ where } \quad
%		 \pt \fcQ = \dot{\fcQ},
\end{equation}

and the product rule to find that $\wm = \fcP\fcQ^{-1}$ satisfies the differential equation

\begin{equation}
\dot{\wm} = -\fcP\fcQ^{-1}\dot{\fcQ}\fcQ^{-1} +\dot{\fcP}\fcQ^{-1}
\end{equation}

with $\pt \fcQ = \dot{\fcQ}$. Then, using (<ref>), we see that this is the differential equation for $\wm$ in

Concerning the symplectic relation in (<ref>), we have

\begin{equation}
\pt (\fcQ^T \fcP - \fcP^T \fcQ) = \dot{\fcQ}^T \fcP + \fcQ^T \dot{\fcP} - \dot{\fcP}^T \fcQ - \fcP^T \dot{\fcQ},
\end{equation}

and by inserting the differential equations of $\fcP, \fcQ$ given in (<ref>), we see that $\fcQ^T \fcP - \fcP^T \fcQ$ is constant.

The same calculation holds for $\pt(\fcQ^* \fcP - \fcP^* \fcQ)$ with $^*$ replaced by $^T$. Since $p,q,A$ and $V$ are real valued, we conclude

\begin{equation}
\pt(\fcQ^* \fcP - \fcP^* \fcQ) = \dot{\fcQ}^* \fcP + \fcQ^* \dot{\fcP} - \dot{\fcP}^* \fcQ - \fcP^* \dot{\fcQ} =0,
\end{equation}

which means that (<ref>) holds true for all times.

§.§ Equations of motion for a general Hamiltonian

The findings of <Ref> for the magnetic  Schrödinger operator $\Ham(t)$ extend to the dynamics for general Hamiltonian operators that are the Weyl quantization of a smooth function
$\cham:\R\times\R^{2\dim}\to\R$ of subquadratic growth, that is, for all
$\alpha\in\N_0^{2\dim}$ with $|\alpha|\ge 2$ there exists $C_\alpha>0$ such that

\begin{equation}\label{eq:general-ham-subquadr}
|\partial^\alpha \cham(t,\varone,\vartwo)| \le C_\alpha
\end{equation}

for all $t\in\R$ and $(\varone,\vartwo)\in\R^{2\dim}$.
Note that the classical magnetic Hamiltonian function (<ref>) is not subquadratic, but our analysis works for both cases.

A first step for the generalization is the construction of a suitable orthonormal basis of the tangent space of a Gaussian wave packet, which is done in <cit.> for the non-magnetic case, where only the modulus squared of the wave packet matters. For convenience, we state the representation formulas of the basis functions that we use. Consider a Gaussian wave packet $\vsol\in\Mf$ of unit norm, $\|\vsol\| =1$.  The family $\{\varphi_n\}_{|n|\le 2}$ with

\beginalign+
φ_0    = ,

φ_e_j    = √(2/) (Q^-1(x-))_j u,

φ_e_j+e_k    = 1/√(δ_kj+1) (2/ (Q^-1(x-))_j  (Q^-1(x-))_k  - (Q^*Q^-T)_j,k ) ,

\endalign+

is an orthonormal basis of the tangent space
$\mathcal T_u\mathcal M$ of $\mathcal M$ at $u$.

For calculating the orthogonal projection to the tangent space, we make use of another representation via the raising and lowering operators $\mathcal A_j^\dagger$ and $\mathcal A_j$. These are the $j$th component of the vector-valued operators

\begin{align}
\mathcal A^\dagger &= \frac{\ii}{\sqrt{2\scp}}\left( P^*\op(\varone-\pscone) -Q^*\op(\vartwo-\psctwo)\right),\\
\mathcal A &= -\frac{\ii}{\sqrt{2\scp}}\left( P^T\op(\varone-\pscone) -Q^T\op(\vartwo-\psctwo)\right),
\end{align}

Using the complete family of Hagedorn functions constructed by the infinite ladder process, we obtain that $\{\varphi_n\}_{|n|\le 2}$  with

\begin{align}\label{eq:onb-ladder-all}
\varphi_0 &= \vsol, \quad
\varphi_{e_j} = \mathcal A_j^\dagger \vsol, \quad
\varphi_{e_k + e_j} = \frac{1}{\sqrt{\delta_{kj} + 1}}\ \mathcal A_j^\dagger \ \mathcal A_k^\dagger \vsol,
\end{align}

see also <cit.> or <cit.>.

Equipped with the orthonormal basis (<ref>) and (<ref>), we can give an explicit formula for the quadratic polynomial generated by the orthogonal projection when acting on a general Hamiltonian operator.

Let $\cham:\R^{2\dim}\to\R$ be smooth and of growth. Let $\vsol\in\Mf$ be a Gaussian wave packet
of unit norm, $\|\vsol\| = 1$, with phase space center $z_0=(\pscone,\psctwo)\in\R^{\dim}\times\R^\dim$. Then,

\begin{equation}
P_\vsol (\op(\cham)\vsol) = p_2 \vsol,
\end{equation}

where $p_2$ is the quadratic polynomial

\begin{align}
& p_2:\R^{\dim}\to\C,\\
& p_2(x) = \beta + b^T(x-\pscone) + \frac12 (x-\pscone)^T B(x-\pscone)
\end{align}

given by the complex coefficients

\begin{align}
\beta &= \langle \cham\rangle_\vsol -\frac{\scp}{4}\, \tr(B \ \IC^{-1}),\\
b &= \begin{pmatrix}\Id & \wm\end{pmatrix} \langle \nabla \cham\rangle_\vsol\in\C^d,\\
B &= \begin{pmatrix}\Id & \wm \end{pmatrix}\langle \nabla^2 \cham\rangle_\vsol
\begin{pmatrix}\Id\\ \wm\end{pmatrix}\in\C^{\dim\times \dim}.
\end{align}

The notation $\langle \clobs\rangle_\vsol = \langle \vsol\mid\op(\clobs)\vsol\rangle$ refers to the expectation value of a quantized smooth observable $\clobs:\R^{2\dim}\to\R^L$ with respect to the Gaussian state $\vsol$.

We use the Hagedorn wave packets $\{\varphi_n\}_{|n|\le 2}$ associated with the Gaussian wave packet $u$ as an orthonormal basis of the tangent space $\mathcal T_u\mathcal M$, see (<ref>) and (<ref>), and write the orthogonal projection as

\[
P_u (\op(h)u) = \sum_{|n|\le 2} \langle\varphi_n\mid\op(h)u\rangle\ \varphi_n.
\]

Starting with the contribution for $n=0$, we have

\begin{align}
\langle\varphi_0\mid\op(h)u\rangle &= \langle u\mid\op(h)u\rangle =\langle h\rangle_u.
\end{align}

For the following, it will be useful to introduce the slim rectangular matrix
$Z = (Q;P)\in\C^{2d\times d}$ with column vectors $Z_1,\ldots,Z_d\in\C^{2d}$ and to write the ladder operators more compactly as

\begin{align}
\mathcal A^\dagger &= \frac{\ii}{\sqrt{2\scp}}\ Z^* J\op(\pstup-\psc) ,\qquad
\mathcal A = -\frac{\ii}{\sqrt{2\scp}}\ Z^TJ\op(\pstup-\psc).
\end{align}

For $n=e_j$ we have by (<ref>), <cit.> that

\begin{align}
\langle\varphi_{e_j}\mid\op(h)u\rangle &=
\langle u\mid \mathcal A_j\op(h)u\rangle
=\langle u\mid [\mathcal A_j,\op(h)]u\rangle.
\end{align}

Since the symbol of $\mathcal A_j$ is linear, we can use pseudodifferential calculus without remainders and obtain that the commutator satisfies

\begin{align}
[\mathcal A_j,\op(h)] &= -\frac{\ii}{\sqrt{2\scp}} \  \left[\op(Z_j^TJ(\pstup-\psc)),\op(\cham)\right]
\\
= -\frac{\ii}{\sqrt{2\scp}} \frac{\scp}{\ii}  \op(\{Z_j^TJ(\pstup-\psc),\cham\})\\
&= \sqrt{\frac{\scp}{2}} \ \op(Z_j^T\nabla h), \label{eq:com}
\end{align}

where we have calculated the Poisson bracket according to

\begin{align}
\{Z_j^TJ z,h\} &= \nabla (Z_j^TJz) \cdot J\nabla h = -Z_j^T\nabla h.
\end{align}

\[
\langle\varphi_{e_j}\mid\op(h)u\rangle = \sqrt{\frac{\scp}{2}} Z_j^T\langle\nabla h\rangle_u.
\]

After summation, we therefore obtain that

\begin{align}
\sum_{j=1}^d \langle\varphi_{e_j}\mid\op(h)u\rangle \,\varphi_{e_j} &=
\sum_{j=1}^d \langle\nabla h\rangle_u^T Z e_j e_j^T Q^{-1}(x-\pscone)u\\
&=\langle\nabla h\rangle_u^T Z Q^{-1}(x-\pscone)u\\
&=\langle\nabla h\rangle_u^T \begin{pmatrix}\Id\\ \wm\end{pmatrix}(x-\pscone)u,
\end{align}

which concludes the computation of the first order contributions.
For the second order wave packets, we analogously compute the projection coefficient as

\begin{align}
\langle \varphi_{e_j+e_k}\mid\op(h)u\rangle
%	&= \frac{1}{\sqrt{\delta_{kj}+1}} \langle \mathcal A_k^\dagger \mathcal A_j^\dagger u\mid\op(h)u\rangle\\
%	&= \frac{1}{\sqrt{\delta_{kj}+1}} \langle \mathcal A_j^\dagger u\mid[\mathcal A_k,\op(h)]u\rangle\\
&= \frac{1}{\sqrt{\delta_{kj}+1}} \langle u\mid[\mathcal A_j,[\mathcal A_k,\op(h)]]u\rangle.
\end{align}

Using (<ref>) twice, we obtain that the double commutator satisfies

\begin{align}
\left[\mathcal A_j,[\mathcal A_k,\op(h)]\right]
&= \sqrt{\frac{\scp}{2}}\left[\mathcal A_j,\op(Z_k^T\nabla h)\right]
%	\\
%	&= \frac{\scp}{2}\,\op(Z_j^T\nabla(Z_k^T\nabla h))\\
%	&
= \frac{\scp}{2}\,\op(Z_j^T\nabla^2 h Z_k).
\end{align}

This implies for the coefficient that

\begin{equation}
\langle \varphi_{e_j+e_k}\mid\op(h)u\rangle  =  \frac{\scp}{2\sqrt{\delta_{kj}+1}}\,Z_j^T\langle\nabla^2 h\rangle_u Z_k.
\end{equation}

We now calculate the sum of all the second order contributions. We have

\begin{align}
%	&
\sum_{|n|=2} \langle \varphi_n\mid\op(h)u\rangle \varphi_n
%	\\
&= \sum_{j=1}^d \sum_{k=1}^j \langle \varphi_{e_j+e_k}\mid\op(h)u\rangle \varphi_{e_j+e_k}\\
&= \sum_{j,k=1}^d \frac{\scp}{2\sqrt{2}}\,Z_j^T\langle\nabla^2 h\rangle_u Z_k \ \frac{1}{\sqrt2} \mathcal A_j^\dagger \mathcal A_k^\dagger u,
\end{align}

where the complete summation over the full square of indices is compensated by a change in
normalisation of the contributions for $j\neq k$. For the part of the sum that generates a constant prefactor for the Gaussian, we have

\begin{align}
-\frac{\eps}{4}\sum_{j,k=1}^d \,Z_j^T\langle\nabla^2 h\rangle_u Z_k \ (Q^*Q^{-T})_{j,k}
&=-\frac{\eps}{4}\tr( Q^*Q^{-T} Z^T \langle\nabla^2 h\rangle_u Z)\\
&=-\frac{\eps}{4}\tr(\begin{pmatrix}\Id & \mathcal C\end{pmatrix} \langle\nabla^2 h\rangle_u \begin{pmatrix}\Id\\ \mathcal C\end{pmatrix} QQ^*).
\end{align}

For the quadratic prefactor, we similarly obtain

\begin{align}
\frac{1}{2}\sum_{j,k=1}^d \,Z_j^T\langle\nabla^2 h\rangle_u Z_k \
\left(Q^{-1}(x-\pscone)\right)_j \,\left(Q^{-1}(x-\pscone)\right)_k
\\
%	&=\frac12 (x-\pscone)^T Q^{-T} Z^T \langle\nabla^2 h\rangle_u Z Q^{-1} (x-\pscone)\\
&=\frac12 (x-\pscone)^T \begin{pmatrix}\Id & \mathcal C\end{pmatrix} \langle\nabla^2 h\rangle_u \begin{pmatrix}\Id\\ \mathcal C\end{pmatrix}
\end{align}

Let $\cham:\R\times\R^{2\dim}\to\R$ be continuous with respect to time $t\in\R$, and smooth, and of subquadratic growth in the sense of (<ref>). Denote $H(t) = \op(\cham(t))$. Then, the time-dependent Schrödinger equation

\[
i\scp\partial_t\sol(t) = H(t)\sol(t),\quad \psi(0) = \psi_0
\]

has a unique solution $\sol(t) = U(t,0)\psi_0$ for all times $t\in\R$ for all square integrable initial data $\psi_0\in L^2(\R^\dim)$, see <cit.> or <cit.>. The corresponding variational Gaussian wave packet obeys the following equations of motion.

Let $\vsol_0\in\Mf$ satisfy (<ref>) and  be given by its parameters $\pos_0, \mom_0, \wm_0, \pha_0$ defined in (<ref>). Then, the parameters of the variational approximation

\begin{equation}
\ii\scp \partial_t \vsol(t) = P_{\vsol(t)} (H(t)\vsol(t)),\quad u(0) = u_0
\end{equation}

satisfy the following set of ordinary differential equations (<ref>)

subject to initial data $(\pos(0),\mom(0), \wm(0), \pha(0)) = (\pos_0, \mom_0, \wm_0, \pha_0)$, where $\cham$ is now the given general subquadratic classical Hamiltonian function.
The Hagedorn parameter matrices of the variational wave packet

\begin{align}
\dot{\fcP} &= -\langle \nabla_{\pos\pos} \cham\rangle_\vsol \fcQ -  \langle\nabla_{\pos\mom} \cham\rangle_\vsol\fcP, \qquad
%			\label{eq:hag-eqmo-general-P}\\
\dot{\fcQ} =   \langle\nabla_{\mom\pos} \cham \rangle_\vsol \fcQ + \langle\nabla_{\mom\mom} \cham\rangle_\vsol \fcP. \label{eq:hag-eqmo-general-Q}
\end{align}

Moreover, the matrix factors $\fcQ$, $\fcP$ are symplectic, provided that the initial matrices $\fcQ_0$, $\fcP_0$ of the factorization $\wm_0 = \fcP_0 \fcQ_0^{-1}$ are symplectic.

We again use (<ref>) and <Ref> and compare the coefficients with respect to the spatial variable $x$. We have one degree of freedom and set, inspired by  (<ref>),

\begin{equation}
\dot{\pos} =\langle \nabla_\mom \cham\rangle_\vsol.
\end{equation}

Now, the claim follows by a direct calculation.

The equations of motion given in <Ref> are indeed a generalization of the magnetic ones derived in <Ref> as we verify next.

In the special space of the magnetic Hamiltonian given in (<ref>) we rediscover the equations of motion (<ref>).

Moreover, if $\scp \to 0$ and averages tend to point evaluations at the center point $\pos$, then the equations (<ref>) and (<ref>) tend to classical equations of motion for a general classical Hamiltonian function $\cham$.

We have that

\begin{align}
\langle \nabla_\mom \cham\rangle_\vsol &= \mom -\langle \mgPot\rangle_\vsol, \quad\text{ and }\quad
-\langle \nabla_\pos \cham\rangle_\vsol
%	&
= -\ii\scp\langle J_A^T\nabla \rangle_\vsol - \langle \nabla \mPot \rangle_\vsol.
\end{align}

Furthermore, it is

\begin{equation}
\nabla^2 \cham(\cdot,\pos, \mom) =
\begin{pmatrix}
\nabla^2 \mPot(\cdot,\pos) - D^2_{\mgPot(\cdot,\pos),\mom}  & -J_A^T\\
-J_A & \Id
\end{pmatrix}
\end{equation}

such that the trace part appearing in (<ref>) contains the terms

-⟨∇_ ⟩_    =        ⟨D^2_,-∇⟩_-⟨∇^2 ⟩_,
- ⟨∇_ ⟩_    =        ⟨J_A^T⟩_,

- ⟨∇_ ⟩_    =         ⟨J_A⟩_,

-⟨∇_ ⟩_    =        -^2.

For the scalar contribution of the projection <Ref> we observe by (<ref>),

\begin{align}
\langle \cham \rangle_\vsol &= -\frac{\scp^2}{2} \langle \Delta \rangle_\vsol +\ii\scp\langle \mgPot\cdot  \nabla \rangle_\vsol +  \langle \mPot \rangle_\vsol\\
%&= \frac12\abs{\mom}^2 +\frac{\eps}{4}\tr\Big((\RC^2 + \IC^2)\ICinv\Big)
%-\langle \mgPot^T \big((\RC+ \ii\IC)((x-\pos)+\mom)\big) \rangle_\vsol+  \langle \mPot \rangle_\vsol\\
&=\frac12\abs{\mom}^2 +\frac{\eps}{4}\tr\Big((\RC^2 + \IC^2)\ICinv\Big)
-\langle \mgPot^T \big(\RC(x-\pos)+\mom\big) \rangle_\vsol+  \langle \mPot \rangle_\vsol.
\end{align}

Finally, we calculate the following trace, appearing in (<ref>), as

\begin{align}
-\tr\Big((\RC^2 + \IC^2)\ICinv\Big) + \tr\Big(\wm^2 \ICinv\Big)
%&= \tr\Big(\big(\RC^2 - \IC^2 + \ii(\IC \RC + \RC\IC)\big)\ICinv\Big)\\
&= -2\,\tr(\IC) + 2\ii \,\tr(\RC)
= 2\ii\, \tr (\wm),
\end{align}

such that, together with

\begin{align}
\mom^T \langle\nabla_{\mom} \cham\rangle_\vsol &= \abs{\mom}^2 - \mom^T \langle \mgPot \rangle_\vsol,
\end{align}

we obtain the differential equation (<ref>).

§ L2-ERROR BOUND

This section is devoted to the of the equations of motion (<ref>)
and the approximation quality of the variational solution in the $L^2$-norm.

We first state

the following lemma which will be used frequently to obtain error bound $\scp$.

We recall that the lower bound on the eigenvalues of $\IC$ was denoted by $\rho >0$.

For any $m \geq 0$ there exists a constant $c_{m}$ such that for all $\scp > 0$ it holds

\begin{align}
&(\pi\scp)^{-\frac{\dim}{4}}\det(\IC)^\frac{1}{4}\Bigl(\int{|x|^{2m} \exp\bigr(-\frac{1}{\scp}x^T \IC x\bigr)}\,\dx\Bigr)^{\frac{1}{2}}\leq c_{m}\Bigl(\frac{\scp}{\rho}\Bigr)^\frac{m}{2} ,
\end{align}

where $c_m$ is independent of $\scp$ and $\rho$.

We now prove the result for (<ref>) and show the boundedness of the parameters solving (<ref>).

We show that the right-hand side of (<ref>) satisfies a local Lipschitz condition with Lipschitz constant independent of $\scp$. To this end it is sufficient if the derivatives parameters $\pos, \mom, \RC, \IC,\pha$ are bounded on a bounded domain. Then, we obtain a local solution and, as in <Ref>, we can show that there is no .

The potentials in the averages of the equations of motion in (<ref>) do not depend on $\scp$.
However, we need to carefully treat

the absolute values of the , since they contain $\scp$ in the denominator.
By the chain rule, it is sufficient to first calculate the derivatives of averages of some arbitrary potential $\refPot$, which is independent of the parameters. Then, the average has the form

\begin{align}
\langle \refPot\rangle_\vsol
%	&= \int \refPot(x)|\vsol(x)|^2\,\dd x\\
&= \frac{\sqrt{\det(\IC)}}{(\pi \scp)^\frac{\dim}{2}} \int \refPot(x) \exp\l(-\frac{1}{\scp}(x-\pos)^T \IC (x-\pos)\r) \,\dd x,
\end{align}

from which we see that, in this case, the average only depends on $\pos$ and $\IC$. Let $\vsol$ be a with $\normLtwo{\vsol} =1$. By (<ref>) we obtain

\begin{align}
& \frac{\sqrt{\det(\IC)}}{(\pi \scp)^\frac{\dim}{2}} \ \partial_{\pos} \exp\l(-\frac{1}{\scp}(x-\pos)^T \IC (x-\pos)\r) \\
&=  \frac{\sqrt{\det(\IC)}}{(\pi \scp)^\frac{\dim}{2}} \ \frac{2}{\scp} \IC(x-\pos) \exp\l(-\frac{1}{\scp}(x-\pos)^T \IC (x-\pos)\r)\\
&=- \nabla |\vsol(x)|^2,
\end{align}

thus, using integration by parts, the derivative of the average to $\pos$ is given by

\begin{align}%\label{Lemma_Abl_nach_q}
\partial_\pos \langle \refPot(x) \rangle_\vsol &= \langle \nabla \refPot(x) \rangle_\vsol.
\end{align}

We continue with derivatives $\IC$. For a differentiable matrix function $F:\R^{\dim\times\dim}\rightarrow\R$ and a general invertible, symmetric matrix $\matrix = (m_{ij})_{i,j=1,\ldots,\dim}$ we define the componentwise derivation matrix

\begin{equation}
\partial_\matrix F(\matrix) \coloneqq (\partial_{m_{ij}} F(\matrix))_{i,j=1,\ldots,\dim} \in \R^{\dim\times\dim}.
\end{equation}

By <cit.> we have

\begin{align}
\partial_{\matrix} a^T \matrix b = a b^T \quad \text{ and }\quad\partial_{\matrix} \det(\matrix)  ={\det(\matrix)} \matrix^{-1}
\end{align}

and consequently,

\begin{align}
\partial_{\matrix} \sqrt{\det(\matrix)} = \frac{1}{2} \sqrt{\det(\matrix)} M^{-1}.
\end{align}

Hence, it follows that

\begin{equation}
\partial_\IC\exp\left(-\frac{1}{\scp}(x-\pos)^T\IC(x-\pos)\right) = -\frac{1}{\scp}(x-\pos)(x-\pos)^T \exp\left(-\frac{1}{\scp}(x-\pos)^T\IC(x-\pos)\right)
\end{equation}

\begin{align}
\partial_{\IC} \langle \refPot(x) \rangle_\vsol &= -\frac{1}{\scp} \bigl\langle (x-\pos)(x-\pos)^T \refPot(x) \bigr\rangle_\vsol
\frac{1}{2}I^{-1}\bigl\langle  \refPot(x) \bigr\rangle_\vsol.
\end{align}

By <Ref> we have $\bigl|\bigl\langle (x-\pos)(x-\pos)^T \refPot(x) \bigr\rangle_\vsol\bigr| \leq C\scp$ for parameters on a bounded domain.

For potentials depending on the parameters, we use again that we are on a bounded domain and that the dependence on $\scp$ of the potentials in (<ref>) is such that $\scp$ does not enter the denominator.

We now turn to the $L^2$-error bound
and adapt the
proof of <cit.> to the magnetic case
and note that the multiplication potential $\mPot$ is already covered. In order to demonstrate the dependence of the constant in the error bound, we carry out the proof for the advection term.

From the proof of <cit.> we know that

\begin{equation}
\normLtwobig{\sol(t)-\vsol(t)} \leq \int_0^t \frac{1}{\scp}\,\normLtwobig{ Hu - P_{u}(Hu)}\,\dd s.
\end{equation}

We write the action of the magnetic Schrödinger operator $H$ on a Gaussian $u$ with width $\wm$ and phase space center $(\pos,\mom)$ as

\[
H\vsol = -\frac{\scp^2}{2}\Delta \vsol + Y \vsol + \mPot \vsol
\]

\begin{equation}\label{eq:advPot}
\advPot \coloneqq -\mgPot \cdot (\wm (x-\pos) + \mom).
\end{equation}

We perform a second order Taylor expansion of the potentials $\advPot$ and $\mPot$ around
the point $\pos$ and denote by $\advPotTaylor$ and $\PotTaylor$ the respective remainders. Then,

\begin{align}
%	H\vsol-P_\vsol(H\vsol) &=
(\Id - P_\vsol)(H\vsol) = (\Id-P_\vsol)(\advPotTaylor\vsol + \PotTaylor\vsol)
\end{align}

\begin{equation}
\normLtwobig{\sol(t)-\vsol(t)} \leq \int_0^t{\frac{1}{\scp}\,\normLtwobig{ \advPotTaylor \vsol + \PotTaylor \vsol}}\,\dd s.
\end{equation}

\begin{equation}\label{eq:Wq_taylor_remainder}
\advPotTaylor = \frac{1}{2}\sum_{|\alpha|=3}{(x-\pos)^\alpha\int_0^1{(1-\theta)^2 \partial^\alpha \advPot\bigl(\pos+\theta (x-\pos)\bigr)}\dd \theta},
\end{equation}

we bound $\normLtwobig{\advPotTaylor u}$ by finding a bound on ${\partial^\alpha \advPot \bigl(q+\theta (x-q)\bigr)}$, which then leads us to

\begin{align}
|\advPotTaylor(x)|^2 \leq C|x-q|^6.
\end{align}

By norm conservation and <Ref> the claim
that $\normLtwobig{\advPotTaylor\vsol} = \mathcal O(\scp^{3/2})$ follows.

For the third derivative of $\partial_{lmn} \advPot$ where $l,m,n = 1,\ldots,d$, we have

\begin{equation}\label{eq:thirdderivY}
\begin{aligned}
\partial_{lmn} \advPot  &= (\partial_{lmn} \mgPot)^T \wm({x}-\pos)+ (\partial_{lmn} \mgPot)^T \mom \\
&\quad + \big((\partial_{lm} \mgPot)^T \wm\big)_n +\big((\partial_{ln} \mgPot)^T \wm\big)_m + \big((\partial_{mn} \mgPot)^T \wm\big)_l,
%	&\quad + \sum_{j=1}^d{\partial_{lm} \mgPot_j \wm_{jn}} + \sum_{j=1}^d{\partial_{ln} \mgPot_j \wm_{jm}} + \sum_{j=1}^d{\partial_{mn} \mgPot_j \wm_{jl}}.
\end{aligned}
\end{equation}

where $\partial_{lmn} \mgPot$ is meant component wise. The term $x-\pos$ in (<ref>) evaluated at $x = \pos+\theta(x-\pos)$ has the form

\begin{equation}
\theta (\partial_{lmn} \mgPot)^T \wm({x}-\pos).
%	\sum_{j,k=1}^d{\partial_{lmn} \mgPot_j \wm_{jk}\theta(x_k-\pos_k)}.
\end{equation}

By <Ref> we gain additional orders of $\scp$, and we thus neglect the first summand in (<ref>). The remaining terms are bounded again using <Ref>.

(b) In the general subquadratic case, we use that the action of a semiclassical pseudodifferential operator on a Gaussian wave packet can be approximated by a polynomial prefactor, see
<cit.>. For any $\ell \in \N$ there exists a polynomial $\mathcal{Q}_\ell$ of degree $\ell$, such that

\begin{equation}
Hu  = \mathcal{Q}_\ell \, \vsol + \O(\scp^{(\ell+1)/2}),
\end{equation}

i.e., we have

\begin{equation}\label{eq:remainder-pot-sq}
Hu - P_{u}(Hu) = \remPot_{\vsol,\ell}\, \vsol + \O(\scp^{(\ell+1)/2}),
\end{equation}

with a remainder potential $\remPot_{\vsol, \ell}$.

We now fix $\ell=2$ and denote the corresponding cubic remainder potential $\remPot_\vsol =\remPot_{\vsol, 2}$. The proof then works along the lines of the magnetic case.

§ EXPECTATION VALUES: PROOF OF <REF>

In this section we adapt the proofs of <cit.> on conservation properties to the , magnetic case. Due to , the energy will not be a conserved quantity.

Let $\sol$ be the exact solution of (<ref>) and $\vsol$ the variational solution (<ref>) such that (<ref>) holds.

The proof of norm conservation and the energy formula can be done in the same way as in <cit.>.

We only show the conservation of total linear and angular momentum.

By <cit.> or <cit.> it is sufficient to show that $\Ham(t)$ commutes with $P$ and $L$, respectively, for each $t \in [0,\etime]$.
By <cit.> it follows that $\linMom \mgPot_{k_j} = 0$ for all $k \in \{1,\ldots,N\}$ and $j\in \{1,2,3\}$. We further calculate

\begin{align}
\linMom(\mgPot\cdot \nabla)\psi	&= \sum_{k=1}^N{\sum_{j=1}^3{(\linMom \mgPot_{k_j}) \partial_{k_j}\psi +  \mgPot_{k_j} \linMom \partial_{k_j}\psi}}= (\mgPot\cdot \nabla)\linMom\psi.
\end{align}

Furthermore, a tedious calculation shows that $ (\mgPot\cdot \nabla)\angMom\psi = \angMom(\mgPot\cdot \nabla)\psi$ if and only if

\begin{align} %\label{A_und_L}
\sum_{l=1}^N{
\begin{pmatrix}
\mgPot_{{l_2}}\partial_{l_3} -\mgPot_{{l_3}}\partial_{l_2} \\
\mgPot_{{l_3}}\partial_{l_1} - \mgPot_{{l_1}}\partial_{l_3}\\
\mgPot_{{l_1}}\partial_{l_2} - \mgPot_{{l_2}}\partial_{l_1}
\end{pmatrix}} \psi =
\sum_{k=1}^N{\sum_{j=1}^3{\sum_{l=1}^N{\begin{pmatrix}
x_{l_2}(\partial_{l_3}\mgPot_{k_j})\partial_{k_j} -x_{l_3}(\partial_{l_2}\mgPot_{k_j})\partial_{k_j} \\
x_{l_3}(\partial_{l_1}\mgPot_{k_j})\partial_{k_j}  - x_{l_1}(\partial_{l_3}\mgPot_{k_j})\partial_{k_j} \\
x_{l_1}(\partial_{l_2}\mgPot_{k_j})\partial_{k_j} -x_{l_2}(\partial_{l_1}\mgPot_{k_j})\partial_{k_j}
\end{pmatrix}}}}\psi
\end{align}

holds true. This condition

is fulfilled if

\begin{equation}
\partial_{l_m}\mgPot_{k_j} = \alpha \delta_{{l_m},{k_j}}\quad\mathrm{und} \quad \mgPot_{{l_n}} = \alpha x_{l_n}
\end{equation}

holds true for some $\alpha \in \R$, $j,n,m \in \{1,2,3\}$, and $k,l\in \{1,\ldots,N\}$ and thus, if $\mgPot(\cdot,x)= \alpha(\cdot) x$ holds.

§ ERROR BOUND FOR AVERAGES OF OBSERVABLES: PROOF OF <REF>

In this section we give the proof of <Ref>. We proceed in three steps:
First, we follow <cit.> and establish an integral representation for the error
that involves a commutator with
the time-evolved observable. Second, we prove Egorov's theorem for the time-evolution
of observables in the general context of magnetic Schrödinger operators. Third, we derive
a expansion of averages with respect to Gaussian wave packets. The combination
of these steps then allows us to prove <Ref>. We  note that the expansion
of the averages is crucial for improving the observable estimate in <cit.>.
This section applies for both the magnetic and the general subquadratic hamiltonian case. For better readability, some arguments will be provided for the magnetic case only, but with natural slight modifications they also apply for the general subquadratic case.

§.§ Error representation

We start with a useful a posteriori representation for the observable error. To this end,
let $\evof(t,s)$ be the evolution family given by <Ref> and $\obs$ an  observable. We introduce the notation

\begin{equation}\label{eq:evof_obs}
\hobs(t,s) \coloneqq \evof(s,t)\obs\evof(t,s), \qquad t, s \in\R.
\end{equation}

Let $\sol$ be the solution of (<ref>) and $\vsol$ the solution of (<ref>). If the initial value $\sol_0 = \vsol_0\in \Mf$ is a with $\normLtwo{\vsol_0}=1$, then the error of the observables takes the form

\begin{align}\label{eq:obs_bound_aux}
&\bigl\langle \sol(t)|\obs\sol(t)\bigr\rangle
\bigl\langle \vsol(t)|\obs \vsol(t)\bigr\rangle
\\
\int_0^t{\frac{1}{\ii\eps}\bigl\langle u(s) \Big|
\left( \overline\remPot_{\vsol(s)}\hobs(t,s)-\hobs(t,s)\remPot_{\vsol(s)}\right)\vsol(s)\bigr\rangle}\,\dd s,
\end{align}

where the remainder potential $\remPot_\vsol:\R^\dim\to\C$ depends on the Gaussian wave
packet $\vsol$. In the general subquadratic case, it has been previously defined in (<ref>).
In the magnetic Schrödinger case, it satisfies

\begin{equation}\label{eq:remPot}
\begin{aligned}
\remPot_\vsol &= \auxPot(\pos) -\langle \auxPot\rangle_\vsol
+ \frac{\scp}{4} \tr(\IC^{-1}\langle \nabla^2 \auxPot\rangle_\vsol)
+ \left(\nabla \auxPot(\pos) - \langle \nabla \auxPot\rangle_\vsol\right)^T(x-\pos)\\
&\quad + \frac{1}{2}(x-\pos)^T\left(\nabla^2\auxPot(\pos)-\langle \nabla^2 \auxPot\rangle_\vsol\right)(x-\pos)
%&\quad +\frac{1}{3!}\sum_{\abs{k}=3} (x-\pos)^k \big(\partial^k \auxPot(\pos) - \langle\partial^k \auxPot \rangle_\vsol\big)
+ R(\auxPot) ,
\end{aligned}
\end{equation}

with $\auxPot = \advPot + \mPot$ defined in (<ref>) and (<ref>), respectively, and $R(\auxPot)$ being the remainder potential of the
quadratic Taylor expansion of $\auxPot$ around the point $\pos$.
For the non-magnetic Schrödinger case $\mgPot = 0$, we have $\advPot = 0$ and

Let $\evof(t,s)$ be the evolution family, such that the exact solution of (<ref>) is given by (<ref>). Using $\sol_0 = \vsol_0$ and $\evof(t,t) = \Id$ we calculate

\begin{align}
\langle u(t)|\obs u(t)\rangle &- \langle \sol(t)|\obs \psi(t)\rangle\\*[1ex]
&= \langle u(t)|\evof(t,t)\obs\evof(t,t) u(t)\rangle - \langle \evof(t,0)u(0)|\obs \evof(t,0)u(0)\rangle
\\*[1ex]
&= \langle u(t)|U(t,t)\obs \evof(t,t) u(t)\rangle - \langle u(0)|U(0,t)\obs \evof(t,0)u(0)\rangle \\*[1ex]
&= \int_0^t{\frac{\partial}{\partial s}\langle u(s)|\underbrace{\evof(s,t)\obs\evof(t,s)}_{=\hobs(t,s)}u(s)\rangle}\,\dd s.
\end{align}

Employing the differential properties of the evolution family, that is,
$\ii\scp\partial_t \evof(t,s) = H(t)\evof(t,s)$ and $-\ii\scp\partial_t \evof(s,t) = \evof(s,t)H(t)$,
we obtain

\begin{align}
\frac{\partial}{\partial s} \hobs(t,s)&= \frac{1}{\ii\eps}\bigl(\Ham(s)\evof(s,t)\obs\evof(t,s)-\evof(s,t)\obs\evof(t,s)\Ham(s)\bigr)\\
&= \frac{1}{\ii\eps}\bigl(\Ham(s)\hobs(t,s)-\hobs(t,s)\Ham(s)\bigr).
\label{eq:Heisenberg-A-tilde}
\end{align}

Since the variational evolution satisfies $\ii\scp\partial_t \vsol(t) = P_{\vsol(t)}\Ham(t)\vsol(t)$, we then have

\begin{align}
\frac{\partial}{\partial s}\langle u(s)|\hobs(t,s)u(s)\rangle &=
\frac{1}{\ii\scp} \left( \langle  (\Id-P_{\vsol(s)})\Ham(s)\vsol(s)| \hobs(t,s) \vsol(s)\rangle \right)\\
&\left.- \langle \vsol(s)| \hobs(t,s) (\Id - P_{\vsol(s)})\Ham(s) \vsol(s)\rangle\right).
\end{align}

We arrive at (<ref>), using that

\begin{align}
(\Id-P_{\vsol(s)})\Ham(s)\vsol(s) &=
X_{\vsol(s)}\vsol(s) - P_{\vsol(s)}(X_{\vsol(s)}\vsol(s))  = W_{u(s)}u(s).
%	  \\
%	 & =:
%	 W_{u(s)}u(s).
\end{align}

The claimed form of the remainder potential $W_{\vsol(s)}:\R^\dim\to\C$ follows
from <cit.>,
since the proof of the projection formula there also applies for
the potential function $X_{\vsol(s)}$ even though it is complex-valued.

§.§ Egorov's theorem

Further, to prove <Ref> we have to establish a variant of Egorov's theorem, which connects the time-evolved quantum observable $\hobs(t,s)$, in case it originates from
a $\obs = \op(\clobs)$, with the evolution map
of the classical Hamiltonian system.

Since the classical Hamiltonian function $\cham(t)$ is subquadratic, the ordinary differential equation (<ref>) has a globally Lipschitz continuous
right-hand side, and the Picard–Lindelöf theorem provides a unique global solution.

Recall that since $\mgPot$ and $\mPot$ are sublinear and subquadratic, respectively, we obtain a unique global solution to the ordinary differential equation (<ref>).
We denote by

\begin{align}
\flow^{t,s} : \R^{2\dim} \to \R^{2\dim},\quad (\varone_s,\vartwo_s) \mapsto (\varone_s(t), \vartwo_s(t))
\end{align}

the classical propagator, which maps initial values at time $s$ to the solution of (<ref>) at time $t$.
For any $\pstup=(\varone,\vartwo)$, it satisfies the evolution equation

\begin{align}\label{eq:reform_cham}
\pt\Phi^{t,s}(\pstup) &= -J(\nabla_\pstup \cham)(t,\Phi^{t,s}(\pstup)),\\
\Phi^{s,s}(\pstup) &= \pstup.
\end{align}

Both in the magnetic and the general subquadratic case, the classical propagator  $\Phi^{t,\tau}$ is a diffeomorphism with inverse $(\Phi^{t,\tau})^{-1} = \Phi^{\tau,t}$.
For time-independent, subquadratic Hamiltonians it is well-established that

\[
\hobs(t,0) = \op\l(\egor{t}{0}\r) + \mathcal O(\scp^2).
\]

However, to the best of our knowledge, in the literature a proof of the Egorov approximation for the non-autonomous case is not available, and the proofs presented for example in
<cit.>, <cit.>, or <cit.> assume time-independent or compactly supported Hamiltonians and thus do not cover our more general situation.
The main difficulties are the of the Hamiltonian operator $\Ham(t)$, which prevents energy conservation,

and the allowed sublinear growth of the observables.

Let $\obs = \op(\clobs)$ be a quantum observable stemming from a smooth, sublinear classical observable $\clobs$ in the sense of <Ref>. Further, let $\auxobs:\R\times\R\times\R^{2d}\to\R$, $(t,s,\pstup)\mapsto \auxobs(t,s,\pstup)$ be defined by

\begin{equation}\label{eq:egorov-obs}
\auxobs(t,s,\pstup) = \egor{t}{s}(\pstup).
\end{equation}

We consider two cases.

* The Hamiltonian operator stems from a classical, subquadratic  function $\cham$. Then, the observable given in (<ref>) is sublinear and for all $\varphi\in L^2(\R^\dim)$ we have

\begin{align}\label{eq:egorov-estimate-h-sq}
\normLtwobig{\bigl(\hobs(t,s) - \op\l(\auxobs(t,s)\r)\bigr)\varphi} \leq
C\, \scp^2\, \e^{C|t-s|}\, \normLtwo{\varphi}
\end{align}

for all $s,t\in\R$.

* The Hamiltonian operator is a magnetic Schrödinger operator. We assume that the observable given in (<ref>) is of time-exponential growth in the following sense. There exists a smooth nonnegative function $\Gamma(t,s)\ge0$ such that for any $\alpha\in\N^{2\dim}$ there exists $C_\alpha>0$ with

\begin{equation}\label{eq:sublinear_time}
|\partial^{\alpha}_\pstup\auxobs(t,s,\pstup)| \le C_{\clobs,\alpha}\ \exp(|\alpha|\,\Gamma(t,s))
\end{equation}

for all $\pstup\in\R^{2\dim}$ and all $t,s\in\R$. Then, for any $\varphi\in L^2(\R^\dim)$ such that $\op(\pstup)\varphi\in L^2(\R^\dim)$, we then have

\begin{align}
\normLtwoBig{\Big(\hobs(t,s)- \op\l(\auxobs(t,s)\r)\varphi\Big)} \le
C\, \scp^2\, \e^{C|t-s|}\, \normLtwo{\op(\pstup)\varphi},
\end{align}

for all $s,t\in \R$.

for all $s,t\in\R$.
The constant $C>0$ depends on derivative bounds of the potentials $\mgPot$, $\mltPot$ and the observable $\clobs$,
but not on $\scp,t,s$. In particular, $C=0$ for $\mgPot$ linear and $\mltPot$ quadratic.

(1) We start by discussing the growth of the function $\auxobs(t,s,\pstup)$ for case (a).
For first order derivatives $(\varone, \vartwo)$ of the classical propagator we have

\begin{align}
D\flowts{t}{s} &= \Id + J^{-1} \int_s^t \nabla^2 \cham(\tau, \flowts{\tau}{s})\, D \flowts{\tau}{s}\,\dd \tau,
\end{align}

and thus

\begin{align}
\norm{\infty}{D\flowts{t}{s}} &\le 1 + \int_s^t \sup_{\pstup\in\R^{2\dim}}\|\nabla_\pstup^2 \cham(\tau,\pstup)\| \norm{\infty}{D \flowts{\tau}{s}}\,\dd \tau.
\end{align}

Since the Hamiltonian function $\cham(t,\cdot)$ is subquadratic, we have

\begin{equation}
\Gamma(t,s) \coloneqq \int_s^t \sup_{\pstup\in\R^{2\dim}}\|\nabla_\pstup^2 \cham(\tau,\pstup)\| \dd \tau<\infty,
\end{equation}

and by Gronwall's lemma

\[
\norm{\infty}{D\Phi^{t,s}}\le  \exp(\Gamma(t,s)).
\]

Moreover, for any $\alpha\in\N^{2\dim}$ with $|\alpha|\ge 1$ there exists a constant $C_\alpha>0$ such that

\[
|\partial^{\alpha}_\pstup\Phi^{t,s}(\pstup)| \le C_\alpha \exp(|\alpha|\,\Gamma(t,s))
\]

for all $t,s\in\R$ and all $\pstup\in\R^{2\dim}$, see <cit.> for
a proof that literally applies to the non-autonomous case. Then, the same argument as for
<cit.> yields that for every $\alpha\in\N^{2\dim}$ with $|\alpha|\ge 1$ there exists a constant $C_{\clobs,\alpha}>0$ such that

\[
|\partial^{\alpha}_\pstup\auxobs(t,s,\pstup)| \le C_{\clobs,\alpha}\ \exp(|\alpha|\,\Gamma(t,s))
\]

for all $t,s\in\R$ and all $\pstup\in\R^{2\dim}$. In particular, $\auxobs(t,s,\cdot)$ is sublinear.

(1) We start by discussing the growth of the function $\auxobs(t,s,\pstup)$. Since the Hamiltonian function $\cham(t,\cdot)$ is subquadratic, we have

\[
\Gamma(t,s) \coloneqq \int_s^t \sup_{z\in\R^{2\dim}}\|\nabla_\pstup^2 \cham(\tau,\pstup)\| \d\tau<\infty,
\]

and there exists for any $\alpha\in\N^{2\dim}$ with $|\alpha|\ge 1$ a constant $C_\alpha>0$ such that

\[
|\partial^{\alpha}_\pstup\Phi^{t,s}(\pstup)| \le C_\alpha \exp(|\alpha|\,\Gamma(t,s))
\]

for all $t,s\in\R$ and all $\pstup\in\R^{2\dim}$, see <cit.> for
a proof that literally applies to the non-autonomous case. Then, the same argument as for
<cit.> yields that for every $\alpha\in\N^{2\dim}$ with $|\alpha|\ge 1$ there exists a constant $C_{\clobs,\alpha}>0$ such that

\[
|\partial^{\alpha}_\pstup\auxobs(t,s,\pstup)| \le C_{\clobs,\alpha}\ \exp(|\alpha|\,\Gamma(t,s))
\]

for all $t,s\in\R$ and all $\pstup\in\R^{2\dim}$. In particular, $\auxobs(t,s,\cdot)$ is sublinear.

(2) Next we compare the operators $\op(\auxobs(t,s))$ and $\hobs(t,s)= \evof(s,t)\obs\evof(t,s)$. Since on the diagonal
$\auxobs(t,t,\cdot) = \clobs$ and $\evof(s,s) = \Id$, we obtain similarly as for (<ref>)

\begin{align}
&  \hobs(t,s) - \op\l(\auxegor{t}{s}\r)\\
=&\int_s^t \evof(s,\tau)
\Big(\frac{\ii}{\scp} \big[\Ham(\tau), \op\l(\auxegor{t}{\tau}\r)\big] +\op\l(\partial_\tau\auxegor{t}{\tau}\r) \Big)
\evof(\tau,s)
\,\dd \tau\\
=&\int_s^t \evof(s,\tau)
\Big(\op(\{\cham(\tau), \auxegor{t}{\tau}\})  +\op\l(\partial_\tau\auxegor{t}{\tau}\r) \Big)
\evof(\tau,s)
\,\dd \tau
\remEgorov(t,s) ,
% \O(\scp^2),
\end{align}

where the last equation relies on the product rule of Weyl quantization <cit.>. Here,

\begin{align}
\{\cham(\tau), \auxegor{t}{\tau}\} &=
% \nabla_\vartwo\cham(\tau)\cdot\nabla_\varone\auxegor{t}{\tau}-
\nabla_\pstup \cham(\tau)\cdot J\nabla_\pstup \auxegor{t}{\tau}
\end{align}

denotes the Poisson bracket of $\cham(\tau)$ and $\auxegor{t}{\tau}$. It remains to show that the integral vanishes and that the remainder $\remEgorov(t,s)$ is of order $\scp^2$.

(3)  For the estimation of the remainder, we use that

\[
\remEgorov(t,s) = \scp^2 \int_s^t \evof(s,\tau) \op(\bm{r}(t,\tau)) \evof(\tau,s)\,\dd \tau,
\]

where $\bm{r}(t,\tau,\cdot)$ is a smooth function depending on the derivatives of the order $\ge 3$ of
the function $\cham(\tau,\cdot)$ and of the sublinear $\auxobs(t,\tau,\cdot)$.
In order to estimate

\begin{align}
\normLtwobig{\rho(t,s)\varphi} &\le  \scp^2 \int_s^t
\normLtwobig{\op(\bm{r}(t,\tau)) \evof(\tau,s)\varphi}\,\dd \tau,
\end{align}

we investigate the above integrand.

If $\cham$ is subquadratic, then, due to the estimates given in (a), for all $\alpha\in\N_0^{2d}$ there exist
$c_{1,\alpha},c_{2,\alpha}>0$ such that

\[
|\partial_\pstup^\alpha \bm{r}(t,\tau,\pstup)| \le c_{1,\alpha}\exp(c_{2,\alpha}|t-\tau|)
\]

for all $t,\tau\in\R$ and $\pstup\in\R^{2\dim}$, and the Calderón–Vaillancourt Theorem, see e.g. <cit.>, provides the claimed constant $C>0$ for part (a).

In the magnetic case, we rewrite the remainder function

\begin{equation}
\bm{r}(\varone,\vartwo)= \bm{r}(\cdot, \cdot, \varone,\vartwo) = \bm{b}_0(\varone,\vartwo) + \bm{b}(\varone,\vartwo)^T\vartwo,
\end{equation}

where $\bm{b}_0 : \R \times \R \times \R^{2\dim} \to \R$ and $\bm{b} : \R \times \R \times \R^{2\dim} \to \R^\dim$ are bounded with all their derivatives. For the first summand, we proceed as in the subquadratic case, using Calderón–Vaillancourt. For the second summand containing an unbounded linearity in $\vartwo$,

we use the product rule and obtain that

\begin{align}
\op\big(\bm{r}_2(\varone, \vartwo)\big) = -\op\big(\bm{b}(\varone,\vartwo)\big)\cdot \ii \scp \nabla + \O(\scp).
\end{align}

Then, the boundedness of $\bm{b}$ provides $C_{\bm{b}}>0$ such that

\begin{align}
\normLtwobig{\op(\bm{b}(t,\tau)) \scp\nabla \big(\evof(\tau,s)\varphi\big)} \le C_{\bm{b}}\,\normLtwobig{ \scp\nabla \big(\evof(\tau,s)\varphi\big)}.
\end{align}

In the next step, we analyse $\normLtwobig{ \scp\nabla \big(\evof(\tau,s)\varphi\big)}$.

In particular, due to the estimates given in (a), for all $\alpha\in\N_0^{2d}$ there exist
$c_{1,\alpha},c_{2,\alpha}>0$ such that

\[
|\partial_\pstup^\alpha \bm{r}(t,\tau,\pstup)| \le c_{1,\alpha}\exp(c_{2,\alpha}|t-\tau|)
\]

for all $t,\tau\in\R$ and $\pstup\in\R^{2\dim}$. Therefore, <Ref> in combination with the Calderón–Vaillancourt Theorem, see e.g. <cit.>, provides the claimed constant $C>0$ such that

\begin{align}
\normLtwobig{\rho(t,s)\varphi} &\le  \scp^2 \int_s^t
\normLtwobig{\op(\bm{r}(t,\tau)) \evof(\tau,s)\varphi}\,\dd \tau\\
&\le C\, \scp^2\, \e^{C|t-s|}\, \normLtwo{\varphi}.
\end{align}

(4) Let $t\ge s$ and set $f(t) =  \op(\pstup) \evof(t,s)\varphi$. We argue as in the proof for
<cit.> and observe that $f(t)$ solves the perturbed magnetic Schrödinger equation

\begin{align}
\ii\scp \pt f(t) &=  \op(\pstup)\Ham(t)\evof(t,s)\varphi =  \Ham(t)f(t) + \delta(t)
\end{align}

with source term

\begin{align}
\delta(t) &= [\op(\pstup), \Ham(t)]\evof(t,s)\varphi = \frac{\scp}{\ii}\begin{pmatrix}
\op(\vartwo - \mgPot(\varone))\\
\op\big(\nabla \mPot(\varone)-\nabla \mgPot(\varone) \cdot\vartwo \big)
\end{pmatrix}
\evof(t,s)\varphi,
\end{align}

where we used the product rule for the second equation. In the same spirit as in step (3), we estimate

\begin{align}
\normLtwobig{\delta(t)} \le C \scp \normLtwobig{\op(\pstup)\evof(t,s)\varphi} = C\scp\normLtwobig{f(t)},
\end{align}

where we exploited the sublinearity of $\mgPot$ and that $\mPot$ is subquadratic.

By the variation of constants formula

followed by Gronwall's lemma, we obtain that

\begin{equation}
\normLtwobig{f(t)} \le \e^{C(t-s)}\normLtwobig{f(s)} = \e^{C(t-s)}\normLtwobig{\op(\pstup)\varphi}.
\end{equation}

(5) In the following step we show that 	$\auxobs$

satisfies the transport equation

\beginalign+
∂_τtτ    = -{(τ),tτ},

tt    = \endalign+

for $\tau\in[s,t]$. Then the integrand in question indeed vanishes, and we obtain

\begin{equation}
\hobs(t,s) = \op(\auxegor{t}{s}) + \mathcal O(\eps^2),
\end{equation}

as claimed. We rewrite the transport equation (<ref>) as

\begin{align}\label{eq:reform_transport}
\partial_\tau \auxegor{t}{\tau} &= J\nabla_\pstup \cham(\tau) \cdot\nabla_\pstup\auxegor{t}{\tau},\\
\auxegor{t}{t} &= \clobs.
\end{align}

The following argument crucially uses that $\Phi^{t,\tau}$ is a diffeomorphism with inverse $(\Phi^{t,\tau})^{-1} = \Phi^{\tau,t}$. We observe that
$\auxobs(t,\tau,\flow^{\tau,t}(\pstup)) = \clobs(\pstup)$ for all $\pstup\in\R^{2\dim}$ and calculate

\begin{align}
0 &= \partial_\tau \clobs(\pstup)\\
&= \partial_\tau \auxobs(t,\tau,\flow^{\tau,t}(\pstup))\\
%&=(\partial_\tau \auxobs)(t,\tau,\flow^{\tau,t}(\pstup))
%+ \partial_\tau\flow^{\tau,s}(\pstup)\cdot(\nabla_\pstup \auxobs)(t,\tau,\flow^{\tau,t}(\pstup))\\
&=(\partial_\tau \auxobs)(t,\tau,\flow^{\tau,t}(\pstup))
-J (\nabla_\pstup \cham)(\tau,\flow^{\tau,t}(\pstup))\cdot(\nabla_\pstup \auxobs)(t,\tau,\flow^{\tau,t}(\pstup)),
%&=\left(\partial_\tau \tilde a -J\nabla_\pstup h\right)(\tau,s,\Phi^{\tau,s}(\pstup)).
\end{align}

where, we used in the last step, the chain rule and

Since $\Phi^{\tau,t}$ is a diffeomorphism, this proves that
$\auxobs(t,\tau)$ indeed solves the transport equation (<ref>).

§.§ Averages with respect to Gaussian wave packets

The a posteriori error representation of <Ref> involves an average with
respect to the variational solution. By Egorov's theorem, <Ref>, the time-evolved quantum observable can be approximated by the classical observable evolved along the classical flow. We therefore derive an asymptotic expansion of averages of operators with respect to Gaussian wave packets. For obtaining this expansion, the following phase space moments will be useful.

We consider a Gaussian $u\in\mathcal M$ of unit norm, $\|\vsol\|=1$, with phase space center $\psc=(\pscone,\psctwo)\in\R^{2\dim}$ and width matrix
$\mathcal C\in\C^{\dim\times\dim}$. We denote by

\begin{equation}\label{eq:rho-ell-wm}
%	\covm = \begin{pmatrix} \IC + \RC \ICinv \RC & -\RC \ICinv\\ \ICinv\RC & \ICinv
%\quad  \text{ and }\quad
\rho_\ell(\wm) =\pi^{-d}
\int_{\R^{2\dim}} \pstup^\ell \,\exp(-\pstup\cdot \covm \pstup) \, \dd \pstup
\quad  \text{ with }\quad
\covm = \begin{pmatrix} \IC + \RC \ICinv \RC & -\RC \ICinv\\ \ICinv\RC & \ICinv
\end{pmatrix},
\end{equation}

where $G \in \R^{2\dim\times 2\dim}$ is symmetric, positive definite and symplectic. Then, for any multi-index $\ell=(\ell_1,\ldots,\ell_{2\dim})\in\N_0^{2\dim}$, we have

\begin{equation}
\langle (\pstup-\psc)^\ell\rangle_\vsol  = \scp^{|\ell|/2} \rho_\ell(\wm).
\end{equation}

If the length $|\ell|$ of the multi-index is odd, then we have $\langle (\pstup-\psc)^\ell\rangle_\vsol = 0$.

The claimed representation becomes evident, when using the Wigner function of the
Gaussian wave packet $\vsol$. The Wigner function of a Gaussian wave packet centered in $\psc$ satisfies

\[
\mathcal W_\vsol(\pstup) = (\pi\scp)^{-\dim} \exp(-\tfrac1\scp (\pstup-\psc)\cdot \covm(\pstup-\psc)),
\]

where the matrix $\covm$

is symplectic, symmetric, positive definite,
see <cit.>. The average of any observable can be written as the phase space integral of the symbol versus the Wigner function,
see for example <cit.>. In particular,

\begin{align}
\langle (\pstup-\psc)^\ell\rangle_\vsol
&= \int_{\R^{2\dim}} (\pstup-\psc)^\ell  \,\mathcal W_u(\pstup) \, \dd \pstup\\\label{eq:gauss_mom}
&= (\pi\scp)^{-\dim} \int_{\R^{2\dim}} (\pstup-\psc)^\ell  \exp(-\tfrac1\scp (\pstup-\psc)\cdot \covm(\pstup-\psc)) \, \dd \pstup\\
&=  \pi^{-\dim} \scp^{|\ell|/2} \int_{\R^{2\dim}} \pstup^\ell  \exp(-\pstup\cdot \covm\pstup) \, \dd \pstup,
\end{align}

where we have used that symplecticity implies $\det(\covm) = 1$.
We observe, that if the length $|\ell|$ of the multi-index is odd, then the above integral vanishes,  and consequently $\langle (\pstup-\psc)^\ell\rangle_\vsol = 0$ as well.

We now use these
moments for expanding Gaussian averages with respect to general observables.

We consider a Gaussian $u\in\mathcal M$ of unit norm, $\|\vsol\|=1$, with phase space center $\psc=(\pscone,\psctwo)\in\R^{2\dim}$ and complex width matrix
$\mathcal C\in\C^{\dim\times\dim}$. Then, for any smooth function
$\clobs:\R^{2\dim}\to\R$ with bounded sixth order derivatives,

\[
\langle \clobs\rangle_\vsol = \clobs(\psc) + \scp f_2(\clobs,\wm) + \scp^2 f_4(\clobs,\wm) + \rho^\scp(\clobs,\wm),
\]

\begin{equation}
f_k(\clobs,\wm) = \sum_{|\ell|=k} \frac{1}{\ell !}\, \partial^\ell\clobs(\psc) \,\rho_\ell(\wm),
\qquad k=2,4.
\end{equation}

The second order contribution satisfies

\begin{equation}
f_2(\clobs,\wm) = \tfrac14 \tr(\nabla^2\clobs(\psc)_\wm \,\IC^{-1})
\end{equation}

\begin{equation}\label{eq:hessian-symbol-wm}
\nabla^2\clobs(\psc)_\wm =
\begin{pmatrix}\Id & \wm^*\end{pmatrix} \nabla^2\clobs(\psc)
\begin{pmatrix}\Id\\ \wm\end{pmatrix}\in\C^{d\times d}.
\end{equation}

The remainder satisfies $|\rho^\scp(\clobs,\wm)|\le C \scp^3$ with a constant $C>0$ that only depends on sixth order derivatives of $\clobs$ as well as on the width matrix $\wm$.

We start by Taylor expanding the symbol around $\psc$ with sixth order remainder,

\begin{equation}
\clobs(\pstup) = \sum_{|k|\le 5} \frac{1}{k!} \partial^k \clobs(\psc) (\pstup-\psc)^k + \bm r_6(\pstup;\psc),
\end{equation}

\begin{equation}
\bm r_6(\pstup;\psc) = \sum_{|k|=6} r_k(\pstup;\psc) (\pstup-\psc)^k, \qquad  r_k(\pstup;\psc) = \frac{6}{k!}\int_0^1 (1-\vartheta)^5 \partial^k \clobs(\psc + \vartheta(\pstup-\psc)) \ \dd \vartheta.
\end{equation}

We have

\begin{equation}
\langle \bm r_6(\pstup;\psc)\rangle_u = \sum_{|k|=6} \int_{\R^{2\dim}} r_k(\pstup;\psc) \,(\pstup-\psc)^k\, \mathcal W_\vsol(\pstup) \,\dd \pstup.
\end{equation}

Therefore, using <Ref>,

\begin{align}
\langle \clobs\rangle_u
%	&= \clobs(\psc) + \sum_{|k|=1}^5 \frac{1}{k!} \partial^k \clobs(\psc) \langle(\pstup-\psc)^k\rangle_u + \langle \bm r_6(\pstup;\psc)\rangle_u\\
%	&= \clobs(\psc) + \sum_{|k|=2,4} \frac{1}{k!} \partial^k \clobs(\psc) \langle(\pstup-\psc)^k\rangle_u + \langle \bm r_6(\pstup;\psc)\rangle_u\\*[1ex]
&=\clobs(\psc) + \scp f_2(\clobs,\wm) + \scp^2 f_4(\clobs,\wm) + \langle \bm r_6(\pstup;\psc)\rangle_u,
\end{align}

and, with the same substitution as in the proof of <Ref>, we bound

\begin{align}
\left|\langle \bm r_6(\pstup;\psc)\rangle_u\right| \le C(\clobs,\wm) \, \scp^3 \quad\text{ with }\quad  C(\clobs,\wm) =
\sum_{|k|=6} \|r_k(\cdot;\psc)\|_\infty\
\abs{\rho_k(\wm)}.
%		\pi^{-\dim}
%		 \int_{\R^{2\dim}} |\pstup|^{|k|} \ \exp(-\pstup\cdot \covm \pstup) \ \dd \pstup,
\end{align}

The constant $C(\clobs,\wm)>0$  depends on fourth order derivatives of $\clobs$ and on the width matrix $\wm$. It remains to rewrite the second order contribution as

\begin{align}
f_2(\clobs,\wm) &=
\pi^{-\dim} \sum_{|\ell|=2} \frac{1}{\ell !} \partial^\ell \clobs(\psc)  \int_{\R^{2d}} (G^{-1/2}\pstup)^\ell
\exp(-|\pstup|^2) \, \dd \pstup \\
&= \frac{1}{2}\ \pi^{-\dim} \int_{\R^{2\dim}} \pstup \cdot G^{-1/2}\nabla^2\clobs(\psc)G^{-1/2}\pstup\
\exp(-|\pstup|^2) \, \dd \pstup \\*[1ex]
&= \frac{1}{4}\, \tr(\nabla^2 \clobs(\psc) G^{-1}).
\end{align}

Since $G$ is symplectic and symmetric, its inverse satisfies

\begin{align}
G^{-1} &= JGJ^{-1}
%	\\
%	&
= \begin{pmatrix}  \ICinv & \ICinv\RC \\ \RC\ICinv &  \IC + \RC \ICinv \RC
\end{pmatrix}.
\end{align}

We decompose the Hessian $\nabla^2 \clobs(\psc)$ in block form as

\begin{equation}\label{eq:hess-block-form}
\nabla^2\clobs(\psc) = \begin{pmatrix}A & B\\ B^T & D\end{pmatrix}.
\end{equation}

Using the cyclicity of the trace, we calculate that

\begin{align}
%	&
\tr(\nabla^2 \clobs(\psc) G^{-1})
%	 \\
%	&= \tr(A\ICinv) + \tr(B\RC\ICinv) + \tr(B^T\ICinv\RC) + \tr(D(\IC + \RC\ICinv\RC))\\
%	&
&= \tr((A+B\RC + \RC B^T + \IC D\IC + \RC D\RC)\ICinv)\\
&= \tr(\nabla^2\clobs(\psc) _{\wm} \, \ICinv),
\end{align}

where $\nabla^2\clobs(\psc) _{\wm}$ was defined in (<ref>) and has the form

\begin{align}
\nabla^2\clobs(\psc) _{\wm}
%&\coloneqq \begin{pmatrix}\Id & \wm^*\end{pmatrix}\nabla^2 \clobs(z_0)
%\begin{pmatrix}\Id\\ \wm\end{pmatrix}\\
&= A+B\wm+ \wm^*B^T + \wm^*D\wm,
\end{align}

which gives the claim.

For our analysis of the observable error, we will use <Ref> also for observables
that are products of two functions. One of the factors will have a controlled expansion, when evaluated in the position center of the variational solution.

In the situation of <Ref> applied to a sublinear classical observable
$\clobs:\R^{2\dim}\to\R$, we additionally consider a smooth and subquadratic function $b^\scp:\R^\dim\to\R$, $x\mapsto b^\scp(x)$. Then,

\begin{equation}
f_2(b^\scp,\wm)  = \tfrac14 \tr(\nabla^{2} b^\scp(\pos)\IC^{-1}).
\end{equation}

* If the function satisfies

\[
b^\scp(\pos), \nabla b^\scp(\pos) = \mathcal O(\scp),
\]

\[
\langle \clobs b^\scp\rangle_u = \clobs(\psc) \left( b^\scp(\pos) +
\scp f_2(b^\scp,\wm)\right) +
\mathcal O(\scp^2).
\]

* If the function satisfies

\[
b^\scp(\pos) = \mathcal O(\scp^2),\quad
\nabla b^\scp(\pos),\nabla^2 b^\scp(\pos) = \mathcal O(\scp),
\]

\begin{align}
\langle \clobs b^\scp\rangle_\vsol &=
\clobs(\psc)\left( b^\scp(\pos) + \scp f_2(b^\scp,\wm) + \scp^2 f_4(b^\scp,\wm)\right)\\*[1ex]
&+\scp F_{1,1}(\clobs,b^\scp,\wm) + \scp^2 F_{1,3}(\clobs,b^\scp,\wm) + \mathcal O(\scp^3)
\end{align}

\begin{align}
\quad
F_{1,n}(\clobs,b^\scp,\wm) =
\sum_{|\ell|=n+1} \sum_{\beta\le\ell, |\beta|=1} \frac{1}{(\ell-\beta)!}\ \partial^\beta\clobs(\psc)\ \partial^{\ell-\beta} b^\scp(\pos)\ \rho_\ell(\wm), \quad n=1,3.
\end{align}

For the trace formula, it is enough to observe that the matrices $B$ and $D$ in the block matrix (<ref>) vanish, since $b^\scp$ only depends on $x$.

For proving the expansions of the averages, we crucially use the Leibniz formula for the $\ell$th derivative of the product, that is,

\[
\partial^\ell(\clobs b^\scp)(\psc) =
\sum_{\beta\le\ell} \binom{\ell}{\beta}\ \partial^\beta\clobs(\psc)\ \partial^{\ell-\beta} b^\scp(\pos)
\]

for any multi-index $\ell\in\N_0^{2\dim}$.

* In the situation of statement (a), we only consider $|\ell| =2$ and obtain

\[
\scp \partial^\ell(\clobs b^\scp)(\psc) = \scp \clobs(\psc) \partial^\ell b^\scp(\pos) + \mathcal O(\scp^2).
\]

Then, <Ref> implies

\begin{align}
\langle \clobs b^\scp\rangle_u &=  \clobs(\psc) b^\scp(\pos) + \scp f_2(\clobs b^\eps,\wm)
+ \mathcal O(\scp^2)\\
&=\clobs(\psc) \left( b^\scp(\pos) + \scp f_2(b^\scp,\wm)\right) +
\mathcal O(\scp^2).
\end{align}

* In the situation of statement (b), we aim for a higher order expansion and need to consider
second and fourth derivatives.
In the same spirit as the proof of part (a), <Ref> implies the claimed expansion of the average
$\langle \clobs b^\scp\rangle_\vsol$.

The estimates of <Ref> also apply to functions $b^\eps(x)$ of the form
$b^\scp(x) = B^\scp(x)\cdot (x-q)$, where $B^\scp(x)\in\R^\dim$ is sublinear with uniform bounds in $\scp$.
A derivative $\partial^\alpha b^\scp(x)$
additively decomposes into a bounded function and the function $\partial^\alpha B^\scp(x)\cdot (x-q)$,
which can be controlled by the arguments used in the proof of <Ref> (a).

§.§ Proof of <Ref>

We now have everything at hand to estimate the error of observables
and to conclude our final main result. In the following proof, we use
<Ref> on the potentials and the representation (<ref>)
of the remainder potential $\remPot_\vsol$
only to the extent that the arguments literally also apply to the dynamics induced by general subquadratic hamiltonians.

Thus, the proof improves known observable error estimates in full generality.

By <Ref> we only have to bound the commutator in
the representation formula (<ref>).

(a) We start by recalling, that in the proof of <Ref>, we have estimated

\begin{align}\label{eq:Wu_L2}
\normLtwo{\remPot_{\vsol(s)}\vsol(s)} = \normLtwo{(\Id-P_{\vsol(s)})\Ham(s)\vsol(s)}
\le C\scp^{3/2}.
\end{align}

(b) We denote $\auxobs(t,s) = \egor{t}{s}$ and expand

\begin{align}
&\frac{1}{\ii\eps}\big\langle  \overline\remPot_{\vsol(s)}\hobs(t,s)-\hobs(t,s)\remPot_{\vsol(s)}\big\rangle_{\vsol(s)}\\
=~&	\frac{1}{\ii\eps}\big\langle  \overline\remPot_{\vsol(s)}\op\l(\auxobs(t,s)\r)-\op\l(\auxobs(t,s)\r)\remPot_{\vsol(s)}\big\rangle_{\vsol(s)} +r_1(s,t).
\end{align}

Using first Cauchy-Schwarz and then (<ref>) together with <Ref> and norm conservation, we bound the remainder by

\begin{equation}%\label{eq_Rechnung_r_1}
\begin{aligned}
|r_1(s,t)| &\leq \frac{2}{\scp}	\normLtwo{\remPot_{\vsol(s)}\vsol(s)} \normLtwobig{\bigl(\hobs(t,s) -\op\l(\auxobs(t,s)\r)\bigr)\vsol(s)}%\\
%&\qquad + \frac{1}{\scp}\normLtwobig{\bigl(\hobs(t,s) -\op\l(\egor{t}{s}\r)\bigr)\remPot_{\vsol(s)} \vsol(s)}\\
%&\leq c\,\normLtwobig{\bigl(\hobs(t,s) -\op\l(\egor{t}{s}\r)\bigr)\vsol(s)} \\
%&\qquad+ c\scp \normweiSob{\remPot_{\vsol(s)} \vsol(s)}{\vaiindex{m}}{j}\\
%&\leq c \scp^2, %\int_0^t\|\psi_{t-s}(s)\|_{L^2, m+d+1, l+j, \eps} \,\dd s.
%		&
\leq c\, \scp^{5/2}.
\end{aligned}
\end{equation}

As in the proof of <Ref>, we use the product rule of Weyl calculus and expand the commutator. For notational simplicity, we suppress the dependence on $t$ and $s$. We have by symmetry of the real part the $L^2$ scalar product and anti-symmetry of the Poisson bracket

\begin{align}
%	&=\big\langle\overline\remPot_{\vsol}\auxobs-\auxobs \remPot_{\vsol}\big\rangle_{\vsol} +
%	\frac{\scp}{2\ii} \big\langle\{\overline\remPot_{\vsol},\auxobs\}-\{\auxobs, \remPot_{\vsol}\}\big\rangle_{\vsol}\\
%	 &\quad+\frac{\scp^2}{8} \big\langle\nabla^2\overline\remPot_{\vsol}J\nabla^2\auxobs J
%	 -\nabla^2\auxobs J \nabla^2\remPot_{\vsol}J\big\rangle_{\vsol}
%	 +	\mathcal O(\scp^3)\\
~=~&\frac{2}{\ii} \,\big\langle\im\remPot_{\vsol}\auxobs\big\rangle_{\vsol}
+ \frac{\scp}{\ii} \big\langle\{\re\remPot_{\vsol},\auxobs\}\big\rangle_{\vsol}
+\frac{\scp^2}{4\ii} \big\langle\nabla^2\im\remPot_{\vsol}J\nabla^2\auxobs J\big\rangle_{\vsol}
+	\mathcal O(\scp^3),
\end{align}

where the constant in $\O(\scp^3)$ depends on phase space derivatives of the remainder potential $\remPot_{\vsol(s)}$ and of $\auxobs(t,s)$ of the order $\ge 3$. 	Since

is sublinear and $\remPot_{\vsol(s)}$ consists of subquadratic summands and a non-subquadratic summand which can be handled by <Ref>,

the Calderón–Vaillancourt Theorem applies for the remainder term.

We will prove below that

\begin{align}
&\big\langle\op\l(\im\remPot_{\vsol}\auxobs\r)\big\rangle_{\vsol} = \O(\scp^3), \label{eq:obs-rem-imaginary-clobs-prod}\\
&\big\langle\op\l(\{\re\remPot_{\vsol},\auxobs\}\r)\big\rangle_{\vsol} = \mathcal O(\scp^2), \label{eq:obs-rem-real-poiss}\\
&\big\langle\op\l(\nabla^2\im\remPot_{\vsol}J\nabla^2\auxobs J\r)\big\rangle_{\vsol} = \mathcal O(\scp),\label{eq:obs-rem-imaginary-hess}
\end{align}

which allows us to conclude that

\begin{align}
\frac{1}{\ii\eps}\big\langle \overline\remPot_{\vsol}\hobs-\hobs\remOp_{\vsol}\big\rangle_{\vsol}
= \O(\scp^2).
\end{align}

In order to do so, we first aim at the application of <Ref> statement (a) for $b^\scp = \partial_j\re\remPot_\vsol$ and statement (b) for $b^\eps = \im\remPot_\vsol$;
see also <Ref> for the non-subquadratic terms in $\remPot_\vsol$.

(d) From now on, we notationally focus on the magnetic Schrödinger case, but the analysis works the same for the general case. We denote the phase space center of the variational Gaussian $\vsol$ by $\psc=(\pscone,\psctwo)$. The width matrix of $\vsol$ is $\wm$ and has imaginary part $\IC$.

We recall that the cubic remainder $R(\auxPot)$ in (<ref>) vanishes together with
its first and second derivatives when evaluated in $\pos$.

We first apply the analysis to the Poisson bracket that involves the real part of the remainder

For any $j=1,\ldots,\dim$, we use (<ref>) and  <Ref> and obtain $	\partial_j\re\remPot_\vsol(\pos) = \mathcal O(\scp)$. Furthermore, by <cit.> we have

\begin{align}
&\nabla \partial_j\re\remPot_\vsol(\pos) = \nabla\partial_j\re\auxPot(q)-\langle\nabla\partial_j\re\auxPot\rangle_u = \mathcal O(\scp),\\
&\nabla^2\partial_j\re\remPot_\vsol(\pos) = \nabla^2\partial_j \re R(\auxPot)(q) =
\nabla^2\partial_j \re \auxPot(q).
\end{align}

Hence, the function $b^\scp = \partial_j\re\remPot_\vsol$ fulfills the assumptions of statement (a) in <Ref>, and together with <Ref>, we obtain

\begin{align}
&\langle\partial_j\re\remPot_\vsol \partial_{p_j}\auxobs\rangle_\vsol
= \mathcal O(\scp^2).
\end{align}

After summation over $j$, we have proven (<ref>).

Similarly, the first and second derivatives of $\im \remPot_\vsol$ satisfy

\begin{equation}\label{eq:deriv-im-remPot}
\nabla \im\remPot_\vsol(\pos),~ \nabla^2\im\remPot_\vsol(\pos) 	= \mathcal O(\scp).
\end{equation}

Moreover, <Ref> implies for the point evaluation of the imaginary part of the remainder potential that

\begin{align}
\im\remPot_\vsol(\pos)
%	&=
%	-\langle\im\auxPot\rangle_\vsol + \frac\scp4 \tr(\langle\nabla^2\im\auxPot\rangle_\vsol\,\IC^{-1}) \\
&= \frac\scp4\, \tr\!\left(\left(\langle\nabla^2\im\auxPot\rangle_\vsol - \nabla^2\im\auxPot(\pos)\right)\IC^{-1}\right) + \mathcal O(\scp^2)
\\ %*[1ex]
= \mathcal O(\scp^2).\label{eq:im-remPot}
\end{align}

At this point, a simple application of <Ref> and (<ref>) yields (<ref>).

(f) The expansions in (<ref>) and (<ref>)

show that $b^\eps = \im\remPot_\vsol$  satisfies the assumptions of statement (b) in <Ref>.

In order to prove (<ref>), we analyse the expansion obtained from <Ref> in two steps, aiming at

\begin{align}\label{eq:a3}
&\im\remPot_\vsol(\pos) + \scp f_2(\im\remPot_\vsol,\wm) + \scp^2 f_4(\im\remPot_\vsol,\wm)
= \mathcal O(\scp^3),\\
\label{eq:na3}
&\scp F_{1,1}(\auxobs,\im\remPot_\vsol,\wm) + \scp^2 F_{1,3}(\auxobs,\im\remPot_\vsol,\wm) = \mathcal O(\scp^3).
\end{align}

(g) We start with proving the first estimate (<ref>). For this, we need a slightly more accurate assessment of
$\im\remPot_\vsol(q)$ than developed previously. Using (<ref>), <Ref>, and
(<ref>), we have

\begin{align}
\im\remPot_\vsol(\pos)
%	&= -\langle \im\auxPot\rangle_\vsol
%	+ \tfrac\scp4\tr(\langle\nabla^2\im\auxPot\rangle_\vsol\IC^{-1})\\*[1ex]
=&-\scp f_2(\im\auxPot,\wm) - \scp^2 f_4(\im\auxPot,\wm)
% \\
%	&\quad
%	+ \scp \sum_{|\ell|=2} \frac{1}{\ell!} \langle\partial^\ell \im\auxPot\rangle_\vsol\rho_\ell(\wm)
%	+ \mathcal O(\scp^3)\\
%	&\quad
+ \scp f_2(\langle \im\auxPot\rangle_\vsol,\wm)
+ \mathcal O(\scp^3)\\
=&  - \scp^2 f_4(\im\auxPot,\wm)
%	\\
%	&\quad
+ \scp^2
%	\sum_{|\ell|=2}\sum_{|k|=2} \frac{1}{k! \ell !}\partial^{k+\ell}\im\auxPot(\pos)\rho_k(\wm)\rho_\ell(\wm)
+ \mathcal O(\scp^3).
\end{align}

Similarly, we obtain for the second term in (<ref>) that

\begin{align}
&\scp f_2(\im\remPot_\vsol,\wm)
%	= \scp\sum_{|\ell|=2}\frac{1}{\ell!}\partial^\ell \im\remPot_\vsol(\pos) \rho_\ell(\wm)\\
%	&= \scp\sum_{|\ell|=2}\frac{1}{\ell!}\left(\partial^\ell \im\auxPot(\pos)
%	- \langle\partial^\ell \im\auxPot\rangle_\vsol\right) \rho_\ell(\wm) \\
%	&
= -\scp^2
%	\sum_{|\ell|=2} \sum_{|k|=2} \frac{1}{\ell ! k!} \partial^{k+\ell}\im\auxPot(\pos)
%	\rho_\ell(\wm)\rho_k(\wm)
+ \mathcal O(\scp^3).
\end{align}

Therefore, $\im\remPot_\vsol(\pos)$ cancels both the contributions from the second and the fourth derivatives, and we have proven (<ref>).

(h) We next target the terms on the left hand side of equation (<ref>), that is,

\begin{align}
\scp F_{1,1}(\auxobs,\im\remPot_\vsol,\wm)
=-\scp^2 \sum_{|k|=2}
F_{1,1}(\auxobs,\partial^k \im\auxPot,\wm) \ \rho_k(\wm)
+ \mathcal O(\scp^3)
\end{align}

\begin{align}
F_{1,3}(\auxobs,\im\remPot_\vsol,\wm) = F_{1,3}(\auxobs,\im\auxPot,\wm).
\end{align}

In <Ref>, we provide the combinatorial argument that shows (<ref>) as a
consequence of Isserlis' theorem on the higher moments of multivariate normal distributions.
Hence, we have proven $\langle \im\remPot_\vsol\auxobs\rangle_\vsol = \mathcal O(\scp^3)$, that is,  (<ref>).

The crucial estimates of the previous proof, namely (<ref>) and (<ref>) are one order worse for the Gaussian approximation,
since it lacks the compensating averaging factors of the remainder potential. Therefore, for
the Gaussians only $\mathcal O(\scp)$ observable accuracy can be expected.

§ GAUSSIAN MOMENTS

By an application of Isserlis' theorem, the fourth order Gaussian moments can
be written as sums of products of second order moments. That is, for a
$2\dim$-dimensional Gaussian random vector

\[
(X_1,\ldots,X_{2\dim})\sim\mathcal N(0,\covm)
\]

with mean zero $0\in\R^{2\dim}$ and covariance matrix $\covm\in\R^{2\dim\times 2\dim}$, the fourth order moments satisfy

\begin{align}
&\E(X_i^4) = 3 g_{ii}^2,\\
&\E(X_i^3 X_j) = 3 g_{ii}g_{ij},\\
&\E(X_i^2X_j^2) = g_{ii}g_{jj} + 2g_{ij}^2\\
&\E(X_i^2X_jX_k) = g_{ii}g_{jk} + 2g_{ij}g_{ik}\\
&\E(X_iX_jX_kX_\ell) = g_{ij}g_{k\ell} + g_{ik}g_{j\ell} + g_{i\ell}g_{jk}
\end{align}

with $i,j,k,\ell\in\{1,\ldots,2\dim\}$. We crucially use this for proving that the fourth order summations
that appeared in the proof of <Ref> can be expressed in terms of second order

For any family $(a_{\beta,m})_{\beta,m}$ of real numbers, indexed by $m\in\N_0^{2\dim}$ and
$\beta\le m$ with $|\beta|=1$, we have

\begin{align}
&\sum_{|m|=4} \sum_{\beta\le m, |\beta|=1} \frac{1}{(m-\beta)!}\ a_{\beta,m}\,\rho_m(\wm)
= \sum_{|k|=2}\sum_{|\ell|=2}\sum_{\beta\le \ell,|\beta|=1} \frac{1}{k!} \ a_{\beta,k+\ell}\,
\rho_k(\wm)\rho_\ell(\wm).
\end{align}

We write a multi-index $m\in\N_0^{2\dim}$ of order $|m|=4$ as

\[
m = \langle j_1\rangle +  \langle j_2\rangle  +  \langle j_3\rangle  +  \langle j_4\rangle
\]

with coordinates $j_1,\ldots,j_4\in\{1,\ldots,2\dim\}$, where
the bracket $\langle j\rangle = e_j$ denotes the $j$th canonical basis vector of $\R^{2\dim}$.
We distinguish five different cases for the order four multi-index $m$.

(a) $m$ has one non-zero component, that is, $m=4\langle j\rangle$ with $j=1,\ldots,2\dim$. Then,

\begin{align}
\frac{1}{(m-\beta)!}\ a_{\beta,m}\,\rho_m(\wm)  &=
\frac{1}{3!}\, a_{\langle j\rangle,4\langle j\rangle} \, 3\rho_{2\langle j\rangle}(\wm)^2\\
&= \frac{1}{k!}\, a_{\beta,k+\ell}\,\rho_k(\wm)\rho_{\ell}(\wm)
\end{align}

with $k=2\langle j\rangle = \ell$ and $\beta=\langle j\rangle$.

(b) $m$ has two different non-zero components, that is,
$m=3\langle j_1\rangle + \langle j_2\rangle$ with $j_1\neq j_2$.
In this case, $m$ dominates two multi-indices $\beta$ of order one, and generates the terms

\begin{align}
&\left( \frac{1}{2!}\,a_{\langle j_1\rangle,m} + \frac{1}{3!}\,a_{\langle j_2\rangle,m}\right)
3\rho_{2\langle j_1\rangle}(\wm)\rho_{\langle j_1\rangle + \langle j_2\rangle}(\wm)\\
&= \frac{1}{2!} \left( a_{\langle j_1\rangle,m} + a_{\langle j_2\rangle,m}\right)
\rho_{2\langle j_1\rangle}(\wm)\rho_{\langle j_1\rangle + \langle j_2\rangle}(\wm)\\
&+\frac{1}{1!} \ a_{\langle j_1\rangle,m}\
\rho_{\langle j_1\rangle + \langle j_2\rangle}(\wm)\rho_{2\langle j_1\rangle}(\wm).
\end{align}

This amounts to the two $(k,\ell)$ pairs

\begin{align}
&k=2\langle j_1\rangle, \ \ell=\langle j_1\rangle+\langle j_2\rangle,\
\beta\in\{\langle j_1\rangle,\langle j_2\rangle\},\\
&k=\langle j_1\rangle+\langle j_2\rangle,\ \ell=2\langle j_1\rangle, \ \beta=\langle j_1\rangle.
\end{align}

In a similar manner, we show that in the cases

(c) $m$ has two identical non-zero components, that is, $m=2\langle j_1\rangle + 2\langle j_2\rangle$ with $j_1\neq j_2$,

(d) $m$ has three non-zero components, that is,
$m=2\langle j_1\rangle + \langle j_2\rangle + \langle j_3\rangle$ with pairwise distinct $j_1,j_2,j_3$,

(e) $m$ has four non-zero components, that is, $m=\langle j_1\rangle + \cdots + \langle j_4\rangle$ with  distinct $j_1,\ldots,j_4$,

we obtain the appropriate format of the resulting summands, that is,

\[
\frac{1}{k!} \sum_{\beta\le \ell,|\beta|=1} a_{\beta,k+\ell}\, \rho_k(\wm)\,\rho_\ell(\wm)
\]

with $k,\ell\in\N_0^{2d}$ such that $|k|=|\ell|=2$ and $k+\ell=m$.

For concluding the proof, we have
to verify that any possible $(k,\ell)$ pairing of order two multi-indices has appeared in one
of the five cases (a)–(e). Let $i_1,\ldots,i_4\in\{1,\ldots,2d\}$ be such that

\[
k=\langle i_1\rangle + \langle i_2\rangle,\quad
\ell = \langle i_3\rangle + \langle i_4\rangle.
\]

The combinatorics of this situation falls into the following five cases:

($\alpha$) All four coordinates agree,
that is, $i_1=\cdots=i_4=:j$.
Then, $k+\ell = 4\langle j\rangle$,
and we recognize the previous case (a).

($\beta$) Three of the four coordinates coincide with each other,

which is case (b).

($\gamma$) The four coordinates form two different pairs,

and we are in case (c).

($\delta$) Two of the four coordinates agree, while the other two are different,

which is case (d).

($\varepsilon$) All four coordinates are distinct as in case (e).

Hence, the combinatorics of the order four multi-indices and the one of pairs of order two multi-indices are the same, and we have indeed proven that the two different summation formats yield the same result as claimed.

§ ACKNOWLEDGMENTS

We thank Clotilde Fermanian-Kammerer and Didier Robert for helpful discussions on the time-dependent
Egorov theorem.