2.1 Molecular Hamiltonian

We describe a molecule in terms of nuclei and electrons. We assume that there are N nuclei, with masses $M_n$ and charge numbers $Z_n$ . The position coordinates are collected in a vector

$$ \begin{align*} x = (x_1,\ldots,x_N)\quad\text{with each } x_n\in{\mathbb R}^3. \end{align*} $$

Similarly, we collect the position coordinates of L electrons in a vector

$$ \begin{align*} y = (y_1,\ldots,y_L)\quad\text{with each } y_l\in{\mathbb R}^3. \end{align*} $$

Electronic mass and elementary charge are denoted by m and e, respectively. The kinetic energy operator

$$ \begin{align*} T = T_N + T_e \end{align*} $$

is the sum of the kinetic energy operators for the nuclei and the electrons,

$$ \begin{align*} T_N = -\sum_{n=1}^N \dfrac{\hbar^2}{2M_n} \Delta_{x_n},\quad T_e = -\sum_{l=1}^L \dfrac{\hbar^2}{2m} \Delta_{y_l}, \end{align*} $$

where $\hbar $ is the reduced Planck constant, which has the physical dimension of an action (length times momentum or energy divided by frequency). The potential energy operator

$$ \begin{align*} V = V_{NN} + V_{Ne} + V_{ee} \end{align*} $$

comprises Coulomb interactions between nuclei and nuclei, nuclei and electrons, as well as electrons and electrons. They are

$$ \begin{align*} V_{NN}(x) = \sum_{1\le m \lt n\le N} \dfrac{Z_m Z_n \,\mathrm{e}^2}{|x_m-x_n|},\quad V_{Ne}(x,y) = -\sum_{n=1}^N \sum_{l=1}^L \dfrac{Z_n \,\mathrm{e}^2}{|x_n-y_l|} \end{align*} $$

and

$$ \begin{align*} V_{ee}(y) = \sum_{1\le k\lt l\le L} \dfrac{\mathrm{e}^2}{|y_k-y_l|}. \end{align*} $$

Adding the kinetic and the potential energy operator, we obtain the molecular Schrödinger operator (or Hamiltonian)

$$ \begin{align*} H_{\mathrm{mol}} = T + V. \end{align*} $$

Before turning to the dynamics induced by the molecular Schrödinger operator, we briefly comment on its basic functional analytic properties.

In the following we distinguish between functions depending on both nuclear and electronic coordinates x and y, functions only depending on the nuclear coordinates x and functions only depending on the electronic coordinates y. We use the shorthand notation

$$ \begin{align*} L^2_{x,y} = L^2({\mathbb R}^{3N+3L}),\quad L^2_x = L^2({\mathbb R}^{3N}),\quad L^2_y = L^2({\mathbb R}^{3L}) \end{align*} $$

for the corresponding Lebesgue spaces of square-integrable complex-valued functions and analogously for Sobolev spaces. The inner products are conjugate linear in the first argument.

We view the molecular Hamiltonian $H_{\mathrm {mol}}$ as an operator on the whole molecular function space $L^2_{x,y}$ .

Kato (Reference Kato1951) proved that the molecular Hamiltonian $H_{\mathrm {mol}}$ is a self-adjoint linear operator with domain $D(H_{\mathrm {mol}}) = H^2_{x,y}$ ; see also Reed and Simon (Reference Reed and Simon1975, Theorem X.16).

2.2 Molecular Schrödinger equation

The quantum motion of a molecule is described by the time-dependent Schrödinger equation,

$$ \begin{align*} \mathrm{i}\hbar\ \dfrac{\partial\Psi}{\partial t} = H_{\mathrm{mol}} \Psi,\quad \Psi |_{t=0} = \Psi_0, \end{align*} $$

which has a unique solution

$$ \begin{align*} \Psi = \Psi(x_1,\ldots,x_N,y_1,\ldots,y_L,t) \end{align*} $$

for all times $t\in {\mathbb R}$ for all square-integrable initial data

$$ \begin{align*} \Psi_0\in L^2_{x,y}. \end{align*} $$

Using the unitary group

$$ \begin{align*} \mathrm{e}^{-\mathrm{i} H_{\mathrm{mol}}\, t/\hbar},\quad t\in{\mathbb R}, \end{align*} $$

associated with the molecular Hamiltonian, the solution of the Schrödinger equation may be written as

$$ \begin{align*} \Psi(\cdot,t) = \mathrm{e}^{-\mathrm{i} H_{\mathrm{mol}} \, t/\hbar} \, \Psi_0. \end{align*} $$

For general square-integrable initial data, it depends continuously on time. For initial data in the domain $D(H_{\mathrm {mol}}) = H^2_{x,y}$ , the solution is continuously differentiable with respect to time.

While it is reassuring that existence and uniqueness of the solution are guaranteed, the direct numerical simulation of the molecular Schrödinger equation is notoriously difficult, or rather intractable, for at least two reasons.

High-dimensionality

A molecular configuration space ${\mathbb R}^{3N+3L}$ is finite-dimensional, but typically high-dimensional. For example, a single CO $_2$ molecule has $N=3$ atoms and $L=22$ electrons. Hence its configuration space has dimension $3N+3L = 75$ . Conventional discretizations of partial differential equations by finite differences or finite element methods scale exponentially with respect to the dimension of the configuration space. Thus they are inappropriate for the molecular Schrödinger equation, and we must search for alternatives.

Multiple scales

Nuclei are much heavier than electrons, which implies that nuclei move on considerably different time and length scales than electrons. For example, electronic motion is typically associated with the order of several tens or hundreds of attoseconds, while chemical reactions take several tens or hundreds of femtoseconds. We recall that

$$ \begin{align*} 1\,\text{attosecond} = 10^{-18}\,\mathrm{s} \;\ll 1\; \text{femtosecond} = 10^{-15}\,\mathrm{s}. \end{align*} $$

Our interest here is in numerical methods that resolve the femtosecond time scale, the time scale of chemical reactions.

2.3 Electronic Schrödinger equation

The mass discrepancy of nuclei and electrons gives the motivation to clamp the nuclei, i.e. to fix a nuclear configuration $x=(x_1,\ldots ,x_N)$ , and to consider the associated electronic Hamiltonian

$$ \begin{align*} H_e(x) = T_e + V_{Ne}(x,\cdot) + V_{ee}. \end{align*} $$

The electronic Hamiltonian depends parametrically on the nuclear configuration and acts on functions of the electronic coordinates $y=(y_1,\ldots ,y_L)$ . The electronic Hamiltonian is a self-adjoint operator with domain $ H^2_y$ that is bounded from below (Kato Reference Kato1951, Theorem 1). In view of the fermionic nature of electrons, the domain $D(H_e(x))$ is usually restricted to the closed subspace of anti-symmetric functions in $H^2_y$ . The electronic Hamiltonian has a real eigenvalue $E(x)$ as the lowest point in its spectrum $\sigma (H_e(x))$ ,

$$ \begin{align*} E(x) = \min\sigma(H_e(x)), \end{align*} $$

which is called the ground state energy. Together with an eigenfunction $\Phi (x,\cdot )$ , it satisfies the stationary electronic Schrödinger equation

$$ \begin{align*} H_e(x)\Phi(x, \cdot) = E(x)\Phi(x,\cdot). \end{align*} $$

The problem of computing a solution to this high-dimensional eigenvalue problem is often referred to as the electronic structure problem.

2.4 Potential energy surfaces

The electronic ground state energy can be characterized using the Courant–Fischer min-max principle (Horn and Johnson Reference Horn and Johnson2013, Theorem 4.2.6), as

$$ \begin{align*} E(x) = \inf_{\dim(U) = 1} \ \sup_{\phi\in U,\,\|\phi\|_{L^2_y} = 1} \langle \phi,\ H_e(x)\phi \rangle_{L^2_y}, \end{align*} $$

where U ranges over the one-dimensional subspaces of $D(H_e(x))$ . For higher electronic eigenvalues a variational characterization is also convenient, and the mth eigenvalue is given by

$$ \begin{align*} E_m(x) = \inf_{\dim(U) = m} \ \sup_{\phi\in U,\,\|\phi\|_{L^2_y} = 1} \langle \phi,\ H_e(x)\phi \rangle_{L^2_y} \quad \text{{for }}m\ge 1. \end{align*} $$

Now U ranges over the m-dimensional subspaces of $D(H_e(x))$ . In this way we also obtain an enumeration of the discrete electronic spectrum that is continuous with respect to the nuclear configuration, regardless of the multiplicities of the eigenvalues that may change when varying x.

The functions that map a nuclear configuration to an electronic eigenvalue,

$$ \begin{align*} {\mathbb R}^{3N}\to{\mathbb R},\quad x\mapsto E_m(x) \end{align*} $$

or

$$ \begin{align*} {\mathbb R}^{3N}\to{\mathbb R},\quad x\mapsto E_m(x) + V_{NN}(x), \end{align*} $$

are called the potential energy surfaces of the molecule. They are important for our basic understanding of chemical reactions.

Hunziker and Günther (Reference Hunziker and Günther1980, Section 6) proved that there is a Lipschitz constant $\beta>0$ such that

$$ \begin{align*} | E_m(x) - E_m(\tilde x) | \le \beta | x-\tilde x| \end{align*} $$

for all $m\ge 1$ and all nuclear configurations $x,\tilde x\in {\mathbb R}^{3N}$ . The Lipschitz constant depends on the number of nuclei and their charge numbers as well as the number of electrons.

The general Lipschitz estimate for the functions $x\mapsto E_m(x)$ can be refined further. Hunziker (Reference Hunziker1986, Corollary of Theorem 1) proved that non-degenerate discrete eigenvalues of the electronic Hamiltonian $H_e(x)$ are analytic in $x = (x_1,\ldots , x_N) \in {\mathbb C}^{3N}$ in a neighbourhood of any nuclear configuration $x\in \mathbb {R}^{3N}$ satisfying $x_k \neq x_n$ if $k\neq n$ .

The two conditions for analyticity are not redundant. As for non-degeneracy, the standard example is provided by the $2\times 2$ matrix

$$ \begin{align*} \begin{pmatrix}\xi_1 & \xi_2\\ \xi_2 & -\xi_1\end{pmatrix}\!{,} \end{align*} $$

which is an entire function of $\xi \in {\mathbb C}^2$ , but has the eigenvalues $\pm (\xi _1^2 + \xi _2^2)^{1/2}$ that are not even differentiable at $\xi =0$ . Such degeneracies occur in many polyatomic molecules; see e.g. Domcke, Yarkony and Köppel (Reference Domcke, Yarkony and Köppel2011). They are referred to as conical intersections, since the graphs of the two eigenvalues form two intersecting, rotationally symmetric cones.

For diatomic molecules, one can adapt the coordinate system such that

$$ \begin{align*} x_1 = \biggl(-\dfrac12 r,0,0\biggr)\quad\text{and}\quad x_2 = \biggl(\dfrac12 r,0,0\biggr). \end{align*} $$

Then the electronic eigenvalues just depend on one single variable, the internuclear distance $r = |x_1-x_2|$ . The simplest molecular ion, the dihydrogen cation H $_2^+$ ( $N=2$ nuclei and $L=1$ electron), shows that analyticity may break down at $r=0$ in a rather subtle way, since its electronic ground state energy has the asymptotic expansion

$$ \begin{align*} E(r) = -2 + \dfrac23 (2r)^2 -\dfrac23 (2r)^3 + \dfrac{22}{135} (2r)^4 - \dfrac{2}{9} (2r)^5 \log(r) + O(r^5) \end{align*} $$

as $r\to 0$ ; see Klaus (Reference Klaus1983).

Remark 2.3. There are two choices for the definition of an electronic Hamiltonian. Setting

$$ \begin{align*} H_e(x) = T_e + V_{Ne}(x,\cdot) + V_{ee}, \end{align*} $$

only the nucleus–electron and electron–electron interactions are taken into account, so the electronic eigenvalues satisfy the above Lipschitz estimate. Alternatively, one may consider

$$ \begin{align*} H_e(x) + V_{NN}(x) \end{align*} $$

as an electronic Hamiltonian acting on functions of the electronic coordinates. This operator has the same spectral properties as $H_e(x)$ , and even the same eigenfunctions. However, since its eigenvalues are

$$ \begin{align*} E_m(x) +V_{NN}(x), \end{align*} $$

they are not globally continuous with respect to x.

2.5 Schrödinger equation for the nuclei via adiabatic approximation

We fix a potential energy surface $E(x)$ and let $\Phi (x,\cdot )\in L^2_y$ be a corresponding eigenfunction of unit norm that depends continuously on x. For fixed nuclear coordinates x, the solution of the time-dependent electronic Schrödinger equation

$$ \begin{align*} \mathrm{i}\hbar\, \partial_t \Psi_e = H_e(x)\Psi_e \end{align*} $$

with initial data $\psi _0(x)\Phi (x,\cdot )$ is given by

$$ \begin{align*} \Psi_e(x,y,t) = \mathrm{e}^{-\mathrm{i} E(x)t/\hbar} \psi_0(x)\cdot \Phi(x,y), \end{align*} $$

and so the solution stays for all times in the subspace of $L^2_{x,y}$ given by

(2.1) $$ \begin{align} {\mathcal V} = \{ v\in L^2_{x,y}\colon v(x,y)=\psi(x)\,\Phi(x,y), \ \psi\in L^2_x \} . \end{align} $$

This motivates the adiabatic approximation to the molecular Schrödinger equation, which is the Galerkin approximation on the subspace ${\mathcal V}$ . Given the initial value $u_0\in {\mathcal V}$ , find $u(t)=u(\cdot ,\cdot ,t)\in {\mathcal V}$ such that

(2.2) $$ \begin{align} \langle \mathrm{i} \hbar\, \partial_t u - H_{\mathrm{mol}}u \mid v \rangle = 0 \quad\text{for all } v \in {\mathcal V}, \quad u(0)=u_0. \end{align} $$

This leads to a Schrödinger equation for the nuclei on the electronic energy surface E: a calculation shows that $u(x,y,t)=\psi (x,t) \,\Phi (x,y)$ , where the nuclear wave function $\psi $ satisfies

$$ \begin{gather*} \mathrm{i}\hbar\,\partial_t \psi = H_N\psi \quad\text{with } H_N = T_N + V_{NN}+E + B_1 + B_2 , \\ B_1 = \sum_{n=1}^N {\hbar\over M_n} \,\operatorname{Im } \langle \nabla_{x_n}\Phi \mid \Phi \rangle_{L^2_y}\cdot (-\mathrm{i}\hbar\nabla_{x_n}) ,\quad B_2 = \sum_{n=1}^N {\hbar^2\over 2M_n} \, \| \nabla_{x_n} \Phi \|_{L^2_y}^2 . \end{gather*} $$

The Hamiltonian $H_N$ acts on functions of only the nuclear coordinates x, with the electronic eigenvalue E as a potential. The last two terms $B_1$ and $B_2$ contain derivatives of the electronic wave function $\Phi $ with respect to the nuclear coordinates x. They are usually neglected in computations, first on practical grounds because they are expensive to compute or simply not available, and second by the formally compelling yet factually dubious argument that they can be neglected because they have the large nuclear masses $M_n$ in the denominator and are of lower differentiation order than the kinetic energy term. The resulting simplified approximation with the Hamiltonian

(2.3) $$ \begin{align} H_{\mathrm{BO}}=T_N+V_{NN}+E \end{align} $$

is known as the time-dependent Born–Oppenheimer approximation. It describes the motion of the nuclei as driven by the potential energy surface E of the electrons. The vast majority of computations in molecular dynamics are based on this approximation.

The term $B_2$ can indeed be safely neglected: it can be shown that its omission introduces an error that is of the same magnitude as the approximation error in the adiabatic approximation.

The term $B_1$ , known as the Berry connection, vanishes for real eigenfunctions $\Phi $ and, more generally, it can be made to vanish by a transformation $\Phi (x,y) \to \mathrm {e}^{\mathrm {i} \theta (x)}\Phi (x,y)$ with $\theta $ satisfying

$$ \begin{align*} \nabla_{x_n}\theta(x) = - {\frac{\hbar}{M_n}} \operatorname{Im } \langle \nabla_{x_n}\Phi \mid \Phi \rangle_{L^2_y}. \end{align*} $$

This transformation of $\Phi $ changes $\psi (x,t)\to \mathrm {e}^{-\mathrm {i} \theta (x)}\psi (x,t)$ . The function $\theta $ is uniquely determined up to a constant if $\Phi $ is a smooth function on all of ${\mathbb R}^{3N}$ or on a simply connected domain, but else $\theta $ is only locally uniquely determined. In the latter case, $B_1$ can cause physical effects that are not captured otherwise; see Berry (Reference Berry1984) and Simon (Reference Simon1983).

2.6 Semiclassical scaling

The success of the adiabatic approximation relies on the smallness of the mass ratio of electrons and nuclei,

$$ \begin{align*} \varepsilon^2 = \dfrac m M \ll 1 \end{align*} $$

with $M=\min _n M_n$ . For ease of presentation, we assume in the following that the masses of the nuclei are all equal: $M_n=M$ for all n. In atomic units ( $\hbar =1$ , $m=1$ , $r=1$ , $e=1$ ) and with the small parameter $\varepsilon $ , the molecular Hamiltonian then takes the form

$$ \begin{align*} H_{\mathrm{mol}}^\varepsilon = -\dfrac{\varepsilon^2}2 \Delta_x + H_e(x) \quad\text{with } H_e(x)= -\dfrac12 \Delta_y - V(x,\cdot). \end{align*} $$

We are interested in solutions to the Schrödinger equation of bounded energy, and in particular of bounded kinetic energy

$$ \begin{align*} \biggl\langle \Psi \mid -\dfrac{\varepsilon^2}2 \Delta_x \Psi \biggr\rangle = \dfrac12\,\| \varepsilon\nabla_x \Psi \|^2 = O(1). \end{align*} $$

For a wave packet $\mathrm {e}^{\mathrm {i} p\cdot x} a(x)$ , with a smooth and fast decaying amplitude function a, this condition corresponds to a momentum $p \sim \varepsilon ^{-1}$ and hence to a velocity $v=p/M\sim \varepsilon $ . Motion of the nuclei over a distance $\sim 1$ can thus be expected on a time scale $\varepsilon ^{-1}$ . We therefore rescale time

$$ \begin{align*} t \to t/\varepsilon , \end{align*} $$

so that, with respect to this new time, nuclear motion over distances $\sim 1$ can be expected to occur at time $\sim 1$ . The molecular Schrödinger equation in the rescaled time then takes the form

(2.4) $$ \begin{align} \mathrm{i}\varepsilon \partial_t\Psi = H_{\mathrm{mol}}^\varepsilon\Psi . \end{align} $$

The Schrödinger equation for the nuclei becomes

(2.5) $$\begin{align} &{\mathrm{i}\varepsilon \partial_t\psi = H_N^\varepsilon\psi\quad\text{with } H_N^\varepsilon=-\dfrac{\varepsilon^2}2\Delta_x + V_{NN} +E +\varepsilon B_1+\varepsilon^2 B_2 ,} &\\ \notag& B_1= \operatorname{Im } \langle \nabla_{x}\Phi \mid \Phi \rangle_{L^2_y}\cdot \widehat p ,\quad B_2 = \dfrac{1}2\,\| \nabla_{x} \Phi \|_{L^2_y}^2 ,&\end{align} $$

with the momentum operator $\widehat p=-\mathrm {i}\varepsilon \nabla _x$ . We are interested in solutions over times $t\sim 1$ .

2.7 Error of the adiabatic approximation

We present an error bound of the adiabatic approximation that was obtained independently by Spohn and Teufel (Reference Spohn and Teufel2001) and by Martinez and Sordoni (Reference Martinez and Sordoni2002). In addition to the smallness of the scaling parameter $\varepsilon $ , a spectral gap condition is required: the eigenvalue $E(x)$ is separated from the remainder of the spectrum $\sigma (H_e(x))$ of the electronic Hamiltonian $H_e(x)$ ,

$$ \begin{align*} \mathrm{dist}( E(x),\, \sigma(H_e(x))\setminus\{ E(x)\} ) \ge \delta>0 {\quad \mbox{for }} x\in{\mathbb R}^{3N}. \end{align*} $$

The Coulomb interactions of the nuclei are mollified to smooth bounded potentials. It is assumed that the potential $V(x,y)$ has a bounded partial derivative with respect to x. Initial data are considered to lie in the approximation space ${\mathcal V}$ of (2.1),

$$ \begin{align*} \Psi_0(x,y) =\psi_0(x)\Phi(x,y) {\quad \mbox{with }} \| H_N^\varepsilon\psi_0 \| \le C_0 , \ \|\psi_0\|=1. \end{align*} $$

We consider the adiabatic approximation

$$ \begin{align*} u(x,y,t)=\psi(x,t)\Phi(x,y) , \end{align*} $$

where $\psi (x,t)$ is the solution of the nuclear Schrödinger equation (2.5) with initial data $\psi _0(x)$ . This is compared with the exact solution $\Psi (t)=\Psi (\cdot ,\cdot ,t)$ of the molecular Schrödinger equation (2.4) with the above initial data $\Psi _0$ .

Theorem 2.1 (space-adiabatic theorem (Teufel Reference Teufel2003))

Under the conditions stated above, the error of the adiabatic approximation is bounded in the $L^2$ -norm by

$$ \begin{align*} \| u(t) - \Psi(t) \| \le c\,(1+t)\,\varepsilon {\quad \mbox{for }} t\ge 0 , \end{align*} $$

where c is independent of $\varepsilon $ and t and initial data as specified above, but depends on the spectral gap $\delta $ .

Teufel (Reference Teufel2003) actually proves a more general result and provides a wealth of related theory. The twisted pseudo-differential calculus developed by Martinez and Sordoni (Reference Martinez and Sordoni2009) also includes Coulomb-type interactions. A short proof of the theorem as stated is given in Lubich (Reference Lubich2008, Section II.2.5).

In the global version stated above, Theorem 2.1 remains valid for the time-dependent Born–Oppenheimer approximation when the terms $B_1$ and $B_2$ in (2.5) are dropped. This is no longer true for local versions, where the eigenfunction $\Phi $ is defined only on a domain that is not simply connected.

The result is related in spirit (though not in proof technique) to the time-adiabatic theorem of Born and Fock (Reference Born and Fock1928) and Kato (Reference Kato1950), which states that, in a quantum system with a slowly time-varying Hamiltonian, a wave function that is initially an eigenfunction approximately remains an eigenfunction of the Hamiltonian over long times.

It is known that the adiabatic approximation generally breaks down near crossings of eigenvalues. One remedy, also covered by the space-adiabatic theory of Teufel (Reference Teufel2003), is to enlarge the approximation space by including several energy bands that are well separated from the remaining ones in the region of physical interest, using

$$ \begin{align*} {\mathcal V} = \biggl\{ u\in L^2_{x,y}\colon u(x,y)=\sum_{j=1}^m\psi_j(x)\Phi_j(x,y) , \ \psi_j\in L^2_x \biggr\}, \end{align*} $$

where $\Phi _j(x,\cdot ) (j=1,\ldots ,m)$ span an invariant subspace of the electronic Hamiltonian $H_e(x)$ . The Galerkin approximation on ${\mathcal V}$ then leads to a system of coupled Schrödinger equations with an $m\times m$ matrix-valued potential.

3.1 Variational approximation by complex Gaussians

We consider the Schrödinger equation in semiclassical scaling,

(3.1) $$ \begin{align} \mathrm{i}\varepsilon\partial_t \psi = H\psi,\quad H = -\dfrac{\varepsilon^2}{2}\Delta_x + V, \end{align} $$

as it emerges from the time-dependent Born–Oppenheimer approximation to the motion of nuclei. We let $d\ge 1$ (often $d\gg 1$ ) denote the dimension of the nuclear configuration space, and assume that the potential $V\colon {\mathbb R}^d\to {\mathbb R}$ is a smooth real-valued function. As before, the $L^2$ inner product is conjugate linear in the first argument.

We seek an approximation $u(\cdot ,t)$ (later written as $u(t)$ for short) to the solution $\psi (t)$ in the manifold of complex-valued Gaussian functions

$$ \begin{align*} & {\mathcal M} = \bigg\{\! u\in L^2({\mathbb R}^d) \mid u(x) = \exp\biggl( \dfrac{\mathrm{i}}{\varepsilon} \biggl(\dfrac{1}{2} (x-q)^T C(x-q) + p^T(x-q) + \zeta\!\biggr) \biggr),\\ & q\in{\mathbb R}^{d},\ p\in{\mathbb R}^{d},\ C=C^T\in{\mathbb C}^{d\times d}, \ \operatorname{Im } C \text{ is positive definite}, \ \zeta\in{\mathbb C}\bigg\}, \end{align*} $$

where we find by a direct computation that, for such an $u\in {\mathcal M}$ of unit $L^2$ -norm, the position and momentum averages equal the parameters q and p, respectively:

$$ \begin{align*} \langle u \mid xu \rangle =q, \quad \langle u \mid - \mathrm{i}\varepsilon \nabla_x u \rangle =p. \end{align*} $$

We invoke the Dirac–Frenkel time-dependent variational approximation principle, which determines the approximate solution

$$ \begin{align*} u(t)\in{\mathcal M} \end{align*} $$

by requiring that the residual of the Schrödinger equation be orthogonal to the tangent space ${\mathcal T}_{u(t)}{\mathcal M}$ of the manifold ${\mathcal M}$ at $u(t)$ : for all times $t\in {\mathbb R}$ ,

(3.2) $$ \begin{align}\notag &\partial_t u(t)\in{\mathcal T}_{u(t)}{\mathcal M}\ \text{ is such that}\\ &\langle v\mid -\mathrm{i}\varepsilon\partial_t u(t) + Hu(t)\rangle = 0\quad \text{for all } v\in{\mathcal T}_{u(t)}{\mathcal M}. \end{align} $$

With the orthogonal projection $P_u\colon L^2({\mathbb R}^d)\to {\mathcal T}_{u}{\mathcal M}$ onto the tangent space ${\mathcal T}_u{\mathcal M}$ , we can rewrite this equivalently as

(3.3) $$ \begin{align} \mathrm{i}\varepsilon \partial_t u = P_u H u. \end{align} $$

We start by having a closer look at the tangent space.

Lemma 3.1 (tangent space)

At every Gaussian function $u\in {\mathcal M}$ , the tangent space equals

$$ \begin{align*} {\mathcal T}_u{\mathcal M} = \{ \varphi u \mid \varphi\text{ is a complex }d\text{-variate polynomial of degree at most }2\}. \end{align*} $$

In particular, ${\mathcal T}_u{\mathcal M}$ is a complex-linear subspace of $L^2({\mathbb R}^d)$ , in the sense that $v\in {\mathcal T}_u{\mathcal M}$ implies $\mathrm {i} v\in {\mathcal T}_u{\mathcal M}$ . Moreover, for all differential operators A of order at most $2$ with constant coefficients, we have $Au\in {\mathcal T}_u{\mathcal M}$ .

Proof. As the tangent space consists of derivatives of paths on ${\mathcal M}$ passing through u, we find that every $v\in {\mathcal T}_u{\mathcal M}$ is of the form

$$ \begin{align*} v(x) = & \dfrac{\mathrm{i}}{\varepsilon} \biggl( -\dot q^T C(x-q) +\dfrac{1}{2} (x-q)^T \dot C (x-q) + \dot p^T(x-q) -p^T\dot q+ \dot\zeta \biggr)u(x) \end{align*} $$

with arbitrary $(\dot q,\dot p)\in {\mathbb R}^{2d}$ , $\dot C={\dot C}^T\in {\mathbb C}^{d\times d}$ , $\dot \zeta \in {\mathbb C}$ . Since every complex polynomial of degree $\le 2$ can be written in the form of the prefactor of the above tangent vector v, we obtain the stated characterization of ${\mathcal T}_u{\mathcal M}$ .□

As a direct corollary we obtain the following property, which is a major motivation for considering localized Gaussians to approximate the wave function $\psi $ .

Proof. If V is quadratic, then

$$ \begin{align*} Hu=-\displaystyle \frac{\varepsilon^2}{2}\Delta_x u + Vu \end{align*} $$

is a quadratic polynomial times u and hence, by Lemma 3.1, it is in the tangent space ${\mathcal T}_u{\mathcal M}$ for $u\in {\mathcal M}$ . It follows that $P_uHu=Hu$ , and hence $u(t)$ and $\psi (t)$ satisfy the same differential equation. □

We remark that this result and its proof extend directly to the case of a time-dependent quadratic potential $V(\cdot ,t)$ .

In the case of a non-quadratic potential, Lemma 3.1 simplifies (3.3) to

(3.4) $$ \begin{align} \mathrm{i}\varepsilon \partial_t u = -\dfrac{\varepsilon^2}{2}\Delta_x u + P_u V u. \end{align} $$

We will study the term $P_uVu$ in detail below, in Proposition 3.14.

3.2 Conservation properties

The variational Gaussian wave packets enjoy norm and energy conservation, as the exact Schrödinger solution does.

Proof. We use the defining relation (3.2) for $\partial _t u(t)\in {\mathcal T}_{u(t)}{\mathcal M}$ . Since the Hamiltonian H is self-adjoint, we have

$$ \begin{align*} \dfrac{\mathrm{d}}{\mathrm{d} t} \langle u(t)\mid H u(t)\rangle &= \langle \partial_t u(t)\mid Hu(t)\rangle + \langle u(t) \mid H \partial_t u(t)\rangle \\ &= 2\,\operatorname{Re }\langle \partial_t u(t) \mid Hu(t)\rangle \\ &= 2\,\operatorname{Re }\langle \partial_t u(t) \mid \mathrm{i}\varepsilon\partial_t u(t)\rangle \\ &= 0, \end{align*} $$

which proves energy conservation. Since $u(t)\in {\mathcal T}_{u(t)}{\mathcal M}$ , we also have

$$ \begin{align*} \dfrac{\mathrm{d}}{\mathrm{d} t} \|u(t)\|^2 &= 2 \,\operatorname{Re } \langle u(t)\mid \partial_t u(t)\rangle \\ &= 2\,\operatorname{Re }\langle u(t)\mid (\mathrm{i} \varepsilon)^{-1}Hu(t)\rangle \\ & = 0, \end{align*} $$

which proves norm conservation. □

The quantum-mechanical total linear momentum and total angular momentum operators, written for the position variables $x=(x_1,\ldots ,x_N)$ with $x_k\in {\mathbb R}^3$ , are

(3.5) $$ \begin{align} \begin{aligned} P & = \sum_{k=1}^N (-i\varepsilon\nabla_{x_k}) \quad\mbox{and}\quad\\ L & = \text{{Hermitian} part of }\ \sum_{k=1}^N x_k \times (-i\varepsilon\nabla_{x_k}), \end{aligned} \end{align} $$

respectively. They both commute with the kinetic energy operator $T=-\frac 12 \varepsilon ^2\Delta _x$ . For potentials that are invariant under translations, that is,

$$ \begin{align*} V(x_1,\ldots, x_N) = V(x_1+ r,\ldots, x_N+ r) \quad \mbox{for all }r\in{\mathbb R}^3, \end{align*} $$

P also commutes with multiplication with V. For potentials invariant under rotations, that is,

$$ \begin{align*} V(x_1,\ldots, x_N) = V(Rx_1,\ldots, Rx_N) \quad \mbox{for all }R\in SO(3), \end{align*} $$

L commutes with multiplication with V, and hence also with the Hamiltonian H. Their averages $\langle \psi \mid P\psi \rangle $ and $\langle \psi \mid L\psi \rangle $ are therefore conserved along solutions $\psi (t)$ to the Schrödinger equation with Hamiltonian $H=T+V$ . As we show next, this conservation property is retained under variational Gaussian wave packet dynamics. We will use the shorthand notation

$$ \begin{align*} \langle A\rangle_u = \langle u \mid A u\rangle \end{align*} $$

for self-adjoint operators A that have a well-defined average with respect to a given Gaussian wave packet $u\in {\mathcal M}$ . For a vector-valued operator (such as P and L above), $\langle A\rangle _u$ is the vector of the averages of the components of A. In the following result, $q(t)=(q_1(t),\ldots ,q_N(t))$ and $p(t)=(p_1(t),\ldots ,p_N(t))$ , where $q_k(t),p_k(t)\in {\mathbb R}^3$ denote the time-dependent position and momentum parameters, respectively, of a variational Gaussian wave packet $u(t)$ .

Proof. We formulate the proof for the angular momentum L. The self-adjoint operator L commutes with the Hamiltonian H, i.e. ${[H,L]=0}$ . Using Lemma 3.1, we find that $Lu\in {\mathcal T}_{u}{\mathcal M}$ for $u\in {\mathcal M}$ . It then follows that

$$ \begin{align*} \dfrac{\mathrm{d}}{\mathrm{d} t}\langle u \mid L u \rangle = 2\, \mathrm{Re}\, \langle Lu \mid \partial_t u \rangle = 2\, \mathrm{Re}\, \biggl\langle Lu \mid \dfrac{1}{\mathrm{i} \varepsilon}\, Hu \biggr\rangle = \biggl\langle u \mid \dfrac{1}{\mathrm{i} \varepsilon}\,[H,L] u \biggr\rangle = 0, \end{align*} $$

where we use (3.2) in the second equality. The proof for linear momentum is the same, since P has the same two properties that were stated for L at the beginning of this proof. Finally, a direct calculation shows that for a Gaussian wave packet, the quantum-mechanical linear and angular momentum equal their classical counterparts for the position parameters q and momenta p.□

3.3 Error in the $L^2$ -norm and in expectation values of observables

We now turn to error bounds which show that, for small $\varepsilon $ , the variational Gaussian approximation stays close to the solution of the Schrödinger equation over time intervals independent of $\varepsilon $ , as long as the Gaussian wave packet remains localized of width $O(\sqrt \varepsilon )$ . Over such time intervals, the error in the $L^2$ -norm is $O(t\sqrt \varepsilon )$ , and the error in observables is $O(t\varepsilon )$ .

We remark that this is an a posteriori error bound since assumption (1) about $\rho $ involves the Gaussian approximation. However, we will see in Section 3.6 that the decay of $\rho $ with time is essentially determined by the Lyapunov exponent of the classical equations of motion $\dot q=p$ , $\dot p=-\nabla V(q)$ .

For the proof of the theorem we first give some simple lemmas that will also be useful later on.

Lemma 3.6 (stability)

Suppose that $e(t)$ satisfies the Schrödinger equation up to a defect $d(t)$ ,

$$ \begin{align*} \partial_t e = \dfrac{1}{\mathrm{i} \varepsilon}He + d, \quad e(0)=e_0. \end{align*} $$

Then

$$ \begin{align*} \| e(t) \| \le \| e_0\| + \int_0^t \| d(s) \| {\,\mathrm{d}} s. \end{align*} $$

Proof. Since H is self-adjoint, $\langle e \mid He\rangle $ is real, and then the Cauchy–Schwarz inequality implies

$$ \begin{align*} \|e\| \dfrac{\mathrm{d}}{\mathrm{d} t}\|e\| &= \dfrac12 \dfrac{\mathrm{d}}{\mathrm{d} t}\|e\|^2 = \operatorname{Re }\langle e\mid \partial_te\rangle \\ &= \operatorname{Re }\biggl\langle e\mid \dfrac{1}{\mathrm{i} \varepsilon}He\biggr\rangle + \operatorname{Re }\langle e \mid d\rangle\\ & = \operatorname{Re }\langle e \mid d\rangle \le \|e\| \, \| d\|. \end{align*} $$

Integration from $0$ to t yields the result. □

It will be useful to have an explicit formula for the norm of a Gaussian wave packet.

Lemma 3.7 (norm)

For $u\in {\mathcal M}$ , we have

$$ \begin{align*} \|u\| = (\pi\varepsilon)^{d/4}\ \det(\operatorname{Im } C)^{-1/4}\ \exp\biggl(-\dfrac1\varepsilon\, \operatorname{Im }\zeta\biggr). \end{align*} $$

Therefore, if $\|u\| = 1$ , then

$$ \begin{align*} \exp\biggl(-\dfrac1\varepsilon\,\operatorname{Im }\zeta\biggr) = (\pi\varepsilon)^{-d/4} \ \det(\operatorname{Im } C)^{1/4}. \end{align*} $$

Proof. We write

$$ \begin{align*} |u(x)|^2 = \exp\biggl( -\dfrac1\varepsilon ( (x-q)^T \operatorname{Im } C(x-q) + 2\,\operatorname{Im }\zeta) \biggr) \end{align*} $$

and obtain

$$ \begin{align*} \int_{{\mathbb R}^d} |u(x)|^2 {\,\mathrm{d}} x = (\pi\varepsilon)^{d/2} \ \det(\operatorname{Im } C)^{-1/2} \ \exp\biggl(-\dfrac2\varepsilon\, \operatorname{Im }\zeta\biggr). \end{align*} $$

The following observation will also be used repeatedly.

Lemma 3.8 (moments)

For arbitrary $m,n\ge 0$ , there exists $c_{m,n}<\infty $ such that the following holds true. For any symmetric positive definite width matrix $\operatorname{Im } C$ with smallest eigenvalue bounded from below by $\rho>0$ and for all positive $\varepsilon>0$ ,

$$ \begin{align*} (\pi\varepsilon)^{-d/4}\ \det(\operatorname{Im } C)^{1/4} \bigg(\int_{{\mathbb R}^d} |x|^{2m} \bigg(1+\dfrac\rho\varepsilon |x|^2\bigg)^n \exp\bigg(-\dfrac1\varepsilon x^T \operatorname{Im } C x) {\,\mathrm{d}} x\bigg)^{1/2}\\ \le c_{m,n} \bigg( \dfrac {\varepsilon} {\rho} \bigg)^{m/2}. \end{align*} $$

Proof. We diagonalize the width matrix

$$ \begin{align*} \operatorname{Im } C = S^T \, \operatorname{diag } (\lambda_1,\ldots,\lambda_d) \ S, \end{align*} $$

with eigenvalues $\lambda _1,\ldots ,\lambda _d>\rho $ and $S\in {\mathbb R}^{d\times d}$ orthogonal. Substituting $y=Sx$ , the Gaussian integral can be rewritten as

$$ \begin{align*} \dfrac{(\lambda_1\cdots\lambda_d)^{1/4}}{(\pi\varepsilon)^{d/4}} \biggl(\int_{{\mathbb R}^d} |y|^{2m} \biggl(1+\dfrac\rho\varepsilon |y|^2\biggr)^n \exp\biggl(-\dfrac1\varepsilon (\lambda_1y_1^2 + \cdots +\lambda_d y_d^2)\biggr) {\,\mathrm{d}} y\biggr)^{1/2}. \end{align*} $$

Then we substitute $z_j=\sqrt {\lambda _j/\varepsilon }y_j$ and obtain the upper bound

$$ \begin{align*} \biggl( \dfrac {\varepsilon} {\rho} \biggr)^{m/2} \pi^{-d/4} \biggl(\int_{{\mathbb R}^d} |z|^{2m} (1+|z|^2)^n \exp(-|z|^2) {\,\mathrm{d}} z\biggr)^{1/2}, \end{align*} $$

which gives the result, since the integral depends only on m and n.□

Proof of Theorem 3.5(a): $L^2$ -error bound

Using (3.3), we have

$$ \begin{align*} \partial_t (u-\psi) &= \dfrac{1}{\mathrm{i} \varepsilon} H(u-\psi) + \dfrac{1}{\mathrm{i} \varepsilon}(P_uH - H)u \\ &= \dfrac{1}{\mathrm{i} \varepsilon}H(u-\psi) - \dfrac{1}{\mathrm{i} \varepsilon} P^\perp_uH u, \end{align*} $$

where $P^\perp _u = \mathrm {Id} - P_u$ is the projection onto the orthogonal complement, and

$$ \begin{align*} {\biggl\|\displaystyle \frac{1}{\mathrm{i} \varepsilon} P^\perp_{u(s)}Hu(s)\biggr\| = \mathrm{dist}\biggl(\displaystyle \frac{1}{\mathrm{i} \varepsilon}H u(s),{\mathcal T}_{u(s)}{\mathcal M}\biggr).} \end{align*} $$

By Lemma 3.6 we therefore have

$$ \begin{align*} \|u(t)-\psi(t)\| \le \int_0^t \mathrm{dist}\biggl(\displaystyle \frac{1}{\mathrm{i} \varepsilon}H u(s),{\mathcal T}_{u(s)}{\mathcal M}\biggr) {\,\mathrm{d}} s. \end{align*} $$

Now we evaluate this general estimate for the Gaussian manifold ${\mathcal M}$ . Let $q\in {\mathbb R}^d$ be the position of $u\in {\mathcal M}$ . We let $U_q\colon {\mathbb R}^d\to {\mathbb R}$ denote the second-order Taylor polynomial of V at q and let $W_q\colon {\mathbb R}^d\to {\mathbb R}$ be the corresponding remainder,

$$ \begin{align*} V = U_q + W_q. \end{align*} $$

By Lemma 3.1, we have $\Delta _x u, U_q u\in {\mathcal T}_u{\mathcal M}$ , and therefore

$$ \begin{align*} \mathrm{dist}\biggl(\dfrac{1}{\mathrm{i} \varepsilon} Hu,{\mathcal T}_u{\mathcal M}\biggr) &= \mathrm{dist}\biggl( \dfrac{1}{\mathrm{i} \varepsilon}W_q u,{\mathcal T}_u{\mathcal M}\biggr) \le \dfrac{1}{\varepsilon} \|W_q u\|. \end{align*} $$

Using the norm conservation of Proposition 3.3 together with Lemma 3.7, we write

$$ \begin{align*} \|W_q u\| &= (\pi\varepsilon)^{-d/4}\ \det(\operatorname{Im } C)^{1/4}\\ &\quad\, \times \biggl(\int_{{\mathbb R}^d} |W_q(x)|^2 \exp\biggl(-\dfrac1\varepsilon(x-q)^T \operatorname{Im } C (x-q)\biggr) {\,\mathrm{d}} x\biggr)^{1/2}. \end{align*} $$

Since $W_q(x)$ is the non-quadratic remainder at q, and by the assumption on the polynomial boundedness of V (of degree n, say), we have $|W_q(x)|^2 \le c_3 |x-q|^6 (1+|x-q|^{2n})$ , and hence Lemma 3.8 implies the estimate

$$ \begin{align*} \| W_q u\| \le c \ \varepsilon^{3/2}, \end{align*} $$

where the constant c depends on $\rho $ . In summary, we obtain

$$ \begin{align*} \|u(t)-\psi(t)\| \le \int_0^t \displaystyle \frac1\varepsilon \|W_{q(s)}u(s)\| {\,\mathrm{d}} s \le c\, t \sqrt\varepsilon, \end{align*} $$

which is the stated $L^2$ -error bound.□

We now turn to the error in observables. In view of Lemma 3.1, for any Gaussian $u\in {\mathcal M}$ there is a quadratic potential $U_u$ such that

$$ \begin{align*} P_u Hu = -\displaystyle \frac{\varepsilon^2}{2}\Delta_x u+ U_u u, \end{align*} $$

that is,

$$ \begin{align*} P_u Vu = U_u u \in {\mathcal T}_u{\mathcal M}. \end{align*} $$

We introduce the remainder potential $W_u=V-U_u$ , so that

(3.6) $$ \begin{align} Hu = P_uHu + W_u u = P_uHP_u u + W_u u, \end{align} $$

where we used that $u\in {\mathcal T}_u{\mathcal M}$ for all $u\in {\mathcal M}$ in the last equality. Theorem 3.5(b) (error in observables) then follows directly from the following two lemmas.

Lemma 3.9 (error in observables)

With the unitary group $U(t)= \exp (-\mathrm {i} tH/\varepsilon )$ and the notation $A(t)=U(t)^*A U(t)$ , we have

$$ \begin{align*} \langle A \rangle_{u(t)} - \langle A \rangle_{\psi(t)} = \int_0^t \biggl\langle \displaystyle \frac 1{\mathrm{i} \varepsilon} [ W_{u(t-s)}, A(s) ] \biggr\rangle_{u(t-s)}{\,\mathrm{d}} s. \end{align*} $$

Lemma 3.10 (commutator bound)

Under the conditions of Theorem 3.5, we have

$$ \begin{align*} \biggl| \biggl\langle \displaystyle \frac 1{\mathrm{i} \varepsilon} [ W_{u(t-s)}, A(s) ] \biggr\rangle_{u(t-s)} \biggr| \le c\,\varepsilon, \quad 0\le s \le t \le \bar t. \end{align*} $$

The proof of Lemma 3.10 requires techniques that are different from those of this section and is therefore deferred to Section 6.7. Lemma 3.9 is proved immediately.

Proof of Lemma 3.9

We compare the outcome of a variational expectation value with the actual one. Writing the true solution as $\psi (t) = U(t)\psi _0$ in terms of the unitary group, we can reformulate the difference of the expectation values in view of $\psi _0=u_0$ as

$$ \begin{align*} \langle \psi(t) \mid A\psi(t)\rangle - \langle u(t) \mid A u(t)\rangle = \int_0^t \displaystyle \frac{\mathrm{d}}{\mathrm{d} s} \langle u(t-s)\mid U(s)^* A U(s) u(t-s)\rangle {\,\mathrm{d}} s. \end{align*} $$

For $A(s) = U(s)^* A U(s)$ , we note that

$$ \begin{align*} \displaystyle \frac{\mathrm{d}}{\mathrm{d} s} A(s) = \displaystyle \frac{1}{\mathrm{i} \varepsilon} ( A(s)H - H A(s)). \end{align*} $$

Therefore, with $P_u^\perp =\mathrm {Id} -P_u$ ,

$$ \begin{align*} &\displaystyle \frac{\mathrm{d}}{\mathrm{d} s} \langle u(t-s)\mid A(s) u(t-s)\rangle \\ &\quad = \displaystyle \frac{1}{\mathrm{i} \varepsilon} ( \langle P_{u(t-s)}H u(t-s)\mid A(s) u(t-s)\rangle \\ &\quad\quad\, + \langle u(t-s)\mid( A(s)H - H A(s)) u(t-s)\rangle\\ &\quad\quad\, - \langle u(t-s)\mid A(s) P_{u(t-s)}H u(t-s)\rangle ) \\ &\quad = \displaystyle \frac{1}{\mathrm{i} \varepsilon} ( \langle u(t-s)\mid A(s) P^\perp_{u(t-s)}H u(t-s)\rangle\\ &\quad\quad\, -\langle P^\perp_{u(t-s)}H u(t-s)\mid A(s) u(t-s)\rangle ). \end{align*} $$

Hence

$$ \begin{align*} &\langle \psi(t) \mid A\psi(t)\rangle - \langle u(t) \mid A u(t)\rangle \\ &\quad = \displaystyle \frac{1}{\mathrm{i} \varepsilon}\int_0^t ( \langle u(t-s)\mid A(s) P^\perp_{u(t-s)}H u(t-s)\rangle \\ & \quad\quad\, -\langle P^\perp_{u(t-s)}H u(t-s)\mid A(s) u(t-s)\rangle ) {\,\mathrm{d}} s. \end{align*} $$

Since $u(t)\in {\mathcal T}_{u(t)}{\mathcal M}$ , we may write

$$ \begin{align*} P_{u(t)}H u(t) = P_{u(t)}H P_{u(t)}u(t) = H_{u(t)}u(t), \end{align*} $$

with $H_{u}=P_uHP_u$ . The above formula for the expectation values then simplifies to

$$ \begin{align*} & \langle \psi(t) \mid A\psi(t)\rangle - \langle u(t) \mid A u(t)\rangle \\[6pt] &\quad = \int_0^t \biggl\langle u(t-s)\mid\displaystyle \frac{1}{\mathrm{i} \varepsilon}[A(s),\ H-H_{u(t-s)}] u(t-s)\biggr\rangle {\,\mathrm{d}} s. \end{align*} $$

Since we have $(H-H_u)u=W_u u$ by (3.6), the result follows.□

3.4 Differential equations for the parameters

In this subsection we derive the following equations of motion for the parameters of variational Gaussian wave packet dynamics.

Theorem 3.11 (equations of motion for parameters)

For Gaussian initial data $u_0\in \mathcal M$ of unit norm, the variational Gaussian approximation (3.2) is obtained with parameters $(q(t),p(t),C(t),\zeta (t))$ that satisfy the following ordinary differential equations:

$$ \begin{align*} \dot{q} &= p,\\[3pt] \dot{p} &= -\langle \nabla_x V\rangle_u,\\[3pt] \dot{C} &= -C^2 - \langle \nabla_x^2 V\rangle_u,\\ \dot{\zeta} &= \displaystyle \frac{1}{2}|p|^2 -\langle V\rangle_u + \displaystyle \frac{\mathrm{i} \varepsilon}{2}\,\mathrm{tr}(C) + \displaystyle \frac{\varepsilon}{4} \,\mathrm{tr}( (\operatorname{Im } C)^{-1}\langle\nabla^2_x V\rangle_u). \end{align*} $$

The first observation is that the parameters change slowly, on the time scale $O(1)$ , as opposed to the highly oscillatory Gaussian wave packet defined by them, which has $O(1)$ changes on time intervals of length $O(\varepsilon )$ . We note further that the first two equations look very similar to the classical equations of motion

$$ \begin{align*} \dot q =p, \quad \dot p= - \nabla V(q), \end{align*} $$

which are obtained from $\dot q =p, \ \dot {p} = -\langle \nabla _x V\rangle _u$ in the limit $\varepsilon \to 0$ if the potential is continuous and the wave packet is localized as stated in assumption (1) of Theorem 3.5. Moreover, the last equation determines $\zeta $ , up to small perturbations, as the action integral that corresponds to the classical Lagrange function $\frac 12|p|^2 - V(q)$ . Only the differential equation for the matrix C seems not to be related to classical mechanics at first sight, but this apparent exemption will be resolved in the next subsection.

For the proof of Theorem 3.11 we first derive auxiliary results that allow us to give an explicit expression for the projected potential term $P_uVu$ in (3.4). We start by setting up an orthonormal basis of the tangent space.

Lemma 3.12 (orthonormal basis of the tangent space)

For a Gaussian $u\in {\mathcal M}$ of unit norm, let the invertible real $d\times d$ matrix Q be a Choleski factor of the inverse of the width matrix, that is,

$$ \begin{align*} \operatorname{Im } C = (QQ^T)^{-1}. \end{align*} $$

For $x\in {\mathbb R}^d$ , we denote

$$ \begin{align*} y = \displaystyle \frac{1}{\sqrt\varepsilon} \,Q^{-1}(x-q), \end{align*} $$

and define the following multivariate scaled Hermite functions $\varphi _k\colon {\mathbb R}^d\to {\mathbb C}$ for multi-indices $k\in {\mathbb N}^d$ with $|k|:=\sum _i k_i\le 2$ :

$$ \begin{align*} \varphi_k(x) = \begin{cases} u(x) & \text{if}\ |k|=0,\\ \sqrt2 y_m u(x) & \text{if}\ |k|=1\ \text{and}\ k_m=1,\\ \displaystyle \frac{1}{\sqrt2}(-1+2y_m^2)u(x) & \text{if} \ |k|=2 \ \text{and}\ k_m=2,\\ 2 y_m y_n u(x) & \text{if}\ |k|=2\ \text{and}\ k_m = k_n = 1. \end{cases} \end{align*} $$

Then $\{\varphi _k\}_{|k|\le 2}$ is an orthonormal basis of ${\mathcal T}_u{\mathcal M}$ .

Proof. By Lemma 3.1, the tangent space consists of all quadratic polynomials times the Gaussian wave packet u. To prove the orthonormality of the functions $\varphi _k$ , we write the Gaussian density as

$$ \begin{align*} |u(x)|^2 &= \exp\biggl(-\displaystyle \frac2\varepsilon\,\operatorname{Im }(\zeta)\biggr) \exp\biggl(-\displaystyle \frac1\varepsilon(x-q)^T \operatorname{Im }(C)(x-q)\biggr)\\ &=\exp\biggl(-\displaystyle \frac2\varepsilon\,\operatorname{Im }(\zeta)\biggr) \exp( -y^T y) \end{align*} $$

and calculate some Gaussian integrals. □

The orthonormal basis $\{\varphi _k\}_{|k|\le 2}$ allows us to explicitly determine projections of functions to the tangent space. We start with functions whose gradient and Hessian vanish on average.

Lemma 3.13. We consider a normalized Gaussian $u\in {\mathcal M}$ and a smooth function $W\colon {\mathbb R}^d\to {\mathbb R}$ with at most polynomial growth satisfying

$$ \begin{align*} \langle \nabla_x W\rangle_u = 0\quad\text{and}\quad \langle\nabla_x^2 W\rangle_u = 0. \end{align*} $$

Then

$$ \begin{align*} P_uWu = \langle W\rangle_u u. \end{align*} $$

Proof. We compute the projection

$$ \begin{align*} P_uWu = \sum_{|k|\le 2} \langle \varphi_k \mid Wu\rangle \varphi_k \end{align*} $$

by examining its summands order by order. The zeroth order provides

$$ \begin{align*} \langle \varphi_0 \mid Wu\rangle \varphi_0 = \langle u \mid W u\rangle u = \langle W\rangle_u u. \end{align*} $$

For the first-order terms, a partial integration shows the vector identity

$$ \begin{align*} \langle \sqrt2 yu \mid Wu\rangle & = \int_{{\mathbb R}^d} \sqrt 2 y W(x) |u(x)|^2 {\,\mathrm{d}} x\\ &=\sqrt{\displaystyle \frac{\varepsilon}{ 2}}Q^T \langle \nabla_x W\rangle_u, \end{align*} $$

since the gradient of the Gaussian density satisfies

$$ \begin{align*} \nabla_x |u(x)|^2 = -\displaystyle \frac{2}{\sqrt\varepsilon}Q^{-T} y |u(x)|^2. \end{align*} $$

Therefore

$$ \begin{align*} \sum_{|k|=1} \langle\varphi_k \mid Wu\rangle\varphi_k = 0. \end{align*} $$

For the second-order terms, we work with the complex symmetric matrix

$$ \begin{align*} \Phi_2(x) = \displaystyle \frac{1}{\sqrt2} (- I_d+2yy^T ) u(x)\in{\mathbb C}^{d\times d}. \end{align*} $$

It contains all the basis functions $\varphi _k(x)$ with $|k|=2$ up to multiplicative factors. The $d(d+1)/2$ basis functions are redundantly placed on the $d^2$ entries of the matrix $\Phi _2(x)$ as follows. The main diagonal carries the $\varphi _k(x)$ with $k_m=2$ for $m=1,\ldots ,d$ . The off-diagonal entries contain the functions $\varphi _k(x)/\sqrt 2$ for $k_m=k_n=1$ with $m\neq n$ . The off-diagonal factors $1/\sqrt 2$ compensate for listing these basis functions twice, and we have

$$ \begin{align*} \sum_{|k|=2} \langle\varphi_k \mid Wu\rangle\varphi_k = \mathrm{tr} ( \langle \Phi_2 \mid Wu\rangle \Phi_2 ). \end{align*} $$

Two partial integrations provide

$$ \begin{align*} \langle \Phi_2 \mid Wu\rangle &= \int_{{\mathbb R}^d} \displaystyle \frac{1}{\sqrt2}(- I_d+2yy^T) W(x)|u(x)|^2 {\,\mathrm{d}} x\\ &= \displaystyle \frac{\varepsilon}{2\sqrt2} Q^T \langle \nabla_x^2 W\rangle_u\,Q, \end{align*} $$

since the Hessian of the Gaussian density satisfies

$$ \begin{align*} \nabla_x^2 |u(x)|^2 = \displaystyle \frac{1}{\varepsilon} Q^{-T}(-2 I_d + 4yy^T)Q^{-1}|u(x)|^2. \end{align*} $$

Therefore

$$ \begin{align*} \sum_{|k|=2} \langle\varphi_k \mid Wu\rangle \varphi_k = 0. \end{align*} $$

For a general potential function, the tangent space projection produces the following quadratic polynomial times the Gaussian.

Proposition 3.14 (potential term projected to the tangent space)

For a normalized Gaussian $u\in {\mathcal M}$ centred at q with width matrix C and for a smooth potential $V\colon {\mathbb R}^d\to {\mathbb R}$ with at most polynomial growth, we have

$$ \begin{align*} P_uVu = \biggl( \alpha+ a^T (x-q) + \displaystyle \frac{1}{2}(x-q)^TA(x-q)\biggr) u, \end{align*} $$

where

$$ \begin{align*} \alpha &= \langle V \rangle _u-\displaystyle \frac{\varepsilon}{4}\mathrm{tr}((\operatorname{Im } C)^{-1}\langle\nabla_x^2 V\rangle_u),\\[6pt] a &= \langle \nabla_x V\rangle_u,\quad\text{and}\quad A = \langle \nabla_x^2 V\rangle_u. \end{align*} $$

Proof. We write the potential as

$$ \begin{align*} V = \langle V\rangle_u + \langle \nabla_x V\rangle_u^T (x-q) + \displaystyle \frac{1}{2}(x-q)^T\langle\nabla_x^2V\rangle_u (x-q) + W \end{align*} $$

and observe that

$$ \begin{align*} P_uVu = \biggl(\langle V\rangle_u + \langle \nabla_x V\rangle_u^T (x-q) + \displaystyle \frac{1}{2}(x-q)^T\langle\nabla_x^2V\rangle_u (x-q) \biggr)u + P_uWu. \end{align*} $$

The function W has the average

$$ \begin{align*} \langle W\rangle_u &= -\displaystyle \frac{1}{2}\langle (x-q)^T \langle\nabla_x^2 V\rangle_u (x-q)\rangle_u\\ &=-\displaystyle \frac{\varepsilon}{4} \, \mathrm{tr} ((\operatorname{Im } C)^{-1}\langle \nabla_x^2 V\rangle_u). \end{align*} $$

For the gradient and Hessian, we have

$$ \begin{align*} \nabla_x V &= \langle \nabla_x V\rangle_u + \langle\nabla_x^2 V\rangle_u (x-q) + \nabla_x W,\\ \nabla_x^2 V &= \langle \nabla_x^2 V\rangle_u + \nabla_x^2W, \end{align*} $$

and therefore

$$ \begin{align*} \langle \nabla_x W\rangle_u = 0\quad\text{and}\quad \langle\nabla_x^2 W\rangle_u = 0. \end{align*} $$

By Lemma 3.13, $P_uWu = \langle W\rangle _u u$ , so that the claimed identity for $P_uVu$ is proved. □

The projection of the potential term is the crucial ingredient for determining the equations of motion for the parameters. The time derivative and the Laplacian only require differentiation of the Gaussian wave packet.

Proof of Theorem 3.11

We start by calculating the first derivatives of the approximate solution $u=u(x,t)$ . We obtain

$$ \begin{align*} -\mathrm{i} \varepsilon\partial_t u = \biggl( -\dot{q }^T C (x-q ) + \displaystyle \frac{1}{2}(x-q )^T \dot{C } (x-q ) + \dot{p }^T(x-q ) - p ^T \dot{q } + \dot{\zeta }\biggr) u \end{align*} $$

and

$$ \begin{align*} -\mathrm{i} \varepsilon\nabla_x u = ( C (x-q ) + p ) u. \end{align*} $$

Hence

$$ \begin{align*} {-\displaystyle \frac{\varepsilon^2}{2}\Delta_x u = \biggl(\displaystyle \frac{1}{2}(x-q )^T C ^2(x-q ) + p ^TC (x-q ) + \displaystyle \frac{1}{2}|p |^2 - \displaystyle \frac{\mathrm{i}}{\varepsilon}{2}\mathrm{tr}(C )\biggr) u. } \end{align*} $$

By Proposition 3.3, the norm of the variational approximation is conserved, and Proposition 3.14 then provides

$$ \begin{align*} P_u(Vu) &= \biggl( \langle V \rangle _u-\displaystyle \frac{\varepsilon}{4}\mathrm{tr}((\operatorname{Im } C )^{-1}\langle\nabla_x^2 V\rangle_u \biggr)u\\ &\quad\, + \langle\nabla_x V\rangle_u^T (x-q )u + \displaystyle \frac{1}{2}(x-q )^T\langle\nabla_x^2V\rangle_u (x-q )u. \end{align*} $$

Hence we solve the variational equation by matching corresponding powers in the polynomials. We set $\dot {q} = p$ . Then the linear, quadratic and constant terms result in the expressions for $\dot p$ , $\dot C$ and $\dot \zeta $ as stated.□

We further note the following bound of the difference between averages and point evaluations.

Lemma 3.15 (averages)

For smooth $V\colon {\mathbb R}^d\to {\mathbb R}$ with bounded second-order derivatives and a normalized complex Gaussian $u\in {\mathcal M}$ with position centre $q\in {\mathbb R}^d$ , we have

$$ \begin{align*} | \langle V\rangle_u - V(q) | \le c\,\varepsilon, \end{align*} $$

where the constant $c>0$ depends on the second-order derivatives of V and the eigenvalues of the matrix $\operatorname{Im } C$ .

Proof. We write

$$ \begin{align*} \langle V\rangle_u &= (\pi\varepsilon)^{-d/2} (\det\operatorname{Im } C)^{1/2} \int_{{\mathbb R}^d} V(x) \exp\biggl(-\displaystyle \frac1\varepsilon (x-q)^T \operatorname{Im } C (x-q)\biggr) {\,\mathrm{d}} x\\ &= \pi^{-d/2} (\det\operatorname{Im } C)^{1/2} \int_{{\mathbb R}^d} V(q+\sqrt\varepsilon y) \exp(-y^T \operatorname{Im } C y) {\,\mathrm{d}} y. \end{align*} $$

We perform a first-order Taylor approximation of $V(q+\sqrt \varepsilon y)$ at q. We then obtain $V(q)$ due to the unit integral and no contribution from $\nabla V(q)$ due to the vanishing first-order moments. The remainder is of order $\varepsilon $ and depends on both V and $\operatorname{Im } C$ .□

3.5 Hagedorn’s parametrization of Gaussian wave packets

We now turn to a different parametrization of complex Gaussians, in which the complex symmetric width matrix C with positive definite imaginary part is factorized into two matrices with remarkable properties.

A real $2d\times 2d$ matrix Y is called symplectic if it satisfies the quadratic relation

$$ \begin{align*} {Y^TJY = J \quad\text{with } J= \begin{pmatrix} 0 & -\mathrm{Id} \\ \mathrm{Id} & 0 \end{pmatrix}\!.} \end{align*} $$

With this notion, we have the following characterization of complex symmetric matrices with positive definite imaginary part.

Lemma 3.16 (symplecticity)

Let Q and P be complex $d\times d$ matrices such that the real $2d\times 2d$ matrix

(3.7) $$ \begin{align} Y = \begin{pmatrix} \operatorname{Re } Q & \operatorname{Im } Q \\ \operatorname{Re } P & \operatorname{Im } P \end{pmatrix} \quad\text{is symplectic}. \end{align} $$

Then Q and P are invertible, and

(3.8) $$ \begin{align} C=PQ^{-1} \end{align} $$

is complex symmetric (i.e. $C=C^T$ ) with positive definite imaginary part

(3.9) $$ \begin{align} \operatorname{Im } C = (QQ^*)^{-1}. \end{align} $$

Conversely, every complex symmetric matrix C with positive definite imaginary part can be written as $C=PQ^{-1}$ with matrices Q and P satisfying (3.7) and (3.9).

We note that the symplecticity condition (3.7) is equivalent to the relations

(3.10) $$ \begin{align} {Q^T P - P^T Q = 0 , \quad Q^* P - P^* Q = 2\mathrm{i} \,\mathrm{Id} .} \end{align} $$

The factorization (3.8) is not unique, since multiplying Q and P from the right with a unitary matrix preserves the relations (3.10).

Proof. The second equation of (3.10), multiplied from the left and the right with a vector $v\in {\mathbb C}^d$ , yields

$$ \begin{align*} (Qv)^*(Pv) - (Pv)^*(Qv) = 2\mathrm{i} \, \| v \|^2 , \end{align*} $$

which shows that $v=0$ is the only vector in the null-space of Q or P. Hence these matrices are invertible. The first equation of (3.10), multiplied from the left with $(Q^{-1})^T$ and from the right with $Q^{-1}$ , gives

$$ \begin{align*} PQ^{-1} - (Q^{-1})^T P^T = 0, \end{align*} $$

which shows that $C=PQ^{-1}$ is complex symmetric. Moreover, we have

$$ \begin{align*} (\operatorname{Im } C) (QQ^*) &= \displaystyle \frac{1}{2{\mathrm{i}}}(PQ^{-1} - (Q^{-1})^*P^*)QQ^* \\ &= \displaystyle \frac{1}{2\mathrm{i}}(PQ^* - (Q^*)^{-1}(P^*Q)Q^*) \\ & ={\mathrm{Id},} \end{align*} $$

where we use the second equation of (3.10) for substituting $P^*Q=Q^*P - 2\mathrm {i} \,\mathrm {Id}$ to obtain the last equality. This gives us (3.9).

Conversely, for a complex symmetric matrix C with positive definite imaginary part, we set $Q=(\operatorname{Im } C)^{-1/2}$ and $P=CQ$ . A direct calculation shows that these matrices satisfy the relations (3.10). □

We now parametrize the variational Gaussian wave packet approximation $u(\cdot ,t)\in {\mathcal M}$ as

(3.11) $$ \begin{align} u(x,t) = \exp\biggl(\!\! \displaystyle \frac{\mathrm{i}}{ \varepsilon} \biggl(\!\!\displaystyle \frac{1}{2} (x-q(t))^T P(t)Q(t)^{-1}(x-q(t)) + p(t)^T(x-q(t)) + \zeta(t)\!\biggr) \!\!\!\biggr), \end{align} $$

where $P(t),Q(t)$ are required to satisfy the symplecticity condition (3.7). We have the following differential equations for the parameters.

Theorem 3.17 (equations of motion for Hagedorn’s parameters)

For Gaussian initial data $u_0\in \mathcal M$ of unit norm, the variational Gaussian approximation (3.2) is obtained with parameters $(q(t),p(t),Q(t),P(t),\zeta (t))$ that satisfy the following ordinary differential equations:

$$ \begin{align*} \dot{q} &= p, \ \,\, \quad\dot{p} = -\langle \nabla_x V\rangle_u,\\ \dot{Q} &= P, \quad\dot{P} = - \langle \nabla_x^2 V\rangle_u Q, \end{align*} $$

and $\zeta $ satisfies the same differential equation as in Theorem 3.11, with $C=PQ^{-1}$ and $(\operatorname{Im } C)^{-1}=QQ^*$ . Moreover, $P(t)$ and $Q(t)$ then have the symplecticity property (3.7) for all t, provided that this is satisfied for the initial values.

The equations for $q,p,\zeta $ are clearly the same as in Theorem 3.11, but the quadratic differential equation for C is replaced by equations for P and Q that are close to the linearization of the classical equations of motion $\dot q=p$ , $\dot p = -\nabla V(q)$ around the solution $(q,p)$ , that is, of the matrix differential equations $\dot Q =P$ , $\dot P = - \nabla ^2 V(q)Q$ that are linear in $(Q,P)$ .

Proof. If $\dot {Q} = P$ and $\dot {P} = - \langle \nabla _x^2 V\rangle _u Q$ , then $C=PQ^{-1}$ satisfies

$$ \begin{align*} \dot C &= - PQ^{-1} \dot Q Q^{-1} + \dot P Q^{-1} \\ &= - PQ^{-1} PQ^{-1} - \langle \nabla_x^2 V\rangle_u Q Q^{-1} \\ & = -C^2 - \langle \nabla_x^2 V\rangle_u , \end{align*} $$

which is the differential equation for C in Theorem 3.11. It remains to show that the symplecticity property of Q and P is conserved. We have

$$ \begin{align*} \dfrac{\mathrm{d}}{\mathrm{d} t}(Q^T P - P^T Q) &= \dot Q^T P + Q^T \dot P - \dot P^T Q - P^T \dot Q \\ &= P^T P -Q^T \langle \nabla_x^2 V\rangle_u Q + Q^T \langle \nabla_x^2 V\rangle_u Q -P^TP \\ & =0, \end{align*} $$

and similarly ${{({{\mathrm {d}}}/{{\mathrm {d} t}})}}(Q^* P - P^* Q) =0$ . This shows that (3.10), and hence (3.7), remains satisfied for all times. □

We refer to the matrices Q and P as position and momentum matrices, respectively, in view of their equations of motion. We will see in the following subsection and in Section 4 that their role goes well beyond merely providing more pleasing differential equations.

3.6 Ehrenfest time

The time scale on which quantum dynamics can be approximated by classical dynamics is known as the Ehrenfest time. The proof of Theorem3.5 shows that the $L^2$ -norm error of the Gaussian approximation is small on a time interval $0\le t \le \bar t$ provided that $\varepsilon ^{1/2} \ll \rho ^{3/2}$ , i.e. $\rho ^{-1} \ll \varepsilon ^{-1/3}$ , where $\rho>0$ is an $\varepsilon $ -independent lower bound of the smallest eigenvalue of the width matrices $\operatorname{Im } C(t)$ for $0\le t \le \bar t$ . By (3.9) we have $\operatorname{Im } C(t) = (Q(t)Q(t)^*)^{-1}$ , and hence the smallest eigenvalue of $\operatorname{Im } C(t)$ is equal to $1/\| Q(t) \|_2^2$ , the inverse square of the matrix $2$ -norm of $Q(t)$ . On the other hand, by Theorem 3.17, $Q(t)$ is up to $O(\varepsilon )$ terms a solution to the linearization of the classical equations of motion $\dot q =p$ , $\dot p =-\nabla V(q)$ . The growth of $\| Q(t) \|_2$ in time is therefore described by the Lyapunov exponent $\lambda _*$ of the classical dynamics: for all $\lambda = \lambda _*+\delta $ with $\delta>0$ , we have $\| Q(t) \|_2 \le c_\delta \,\mathrm {e}^{\lambda t}$ . The Lyapunov exponent thus determines the Ehrenfest time: for all $\delta>0$ ,

$$ \begin{align*} \log \rho^{-1} \le \log \max_{0\le t \le \bar t} \| Q( t) \|_2^2 \le 2(\lambda_*+\delta) \bar t + 2 \log c_\delta. \end{align*} $$

The Lyapunov exponent $\lambda _*$ is always non-negative for classical mechanics. For $\lambda _*>0$ (where we can choose $\delta =\lambda _*$ ), the condition $\rho ^{-1} \ll \varepsilon ^{-1/3}$ thus yields the condition that $\lambda _* \bar t$ be logarithmic in $\varepsilon ^{-1}$ .

The Lyapunov exponent is zero for integrable systems, for which $\| Q(t) \|_2$ grows linearly with t, i.e. $\| Q(t) \|_2\le c t$ , so that for integrable systems

$$ \begin{align*} \rho^{-1} \le \max_{0\le t \le \bar t} \| Q( t) \|_2^2 \le (c \bar t)^2. \end{align*} $$

The condition $\rho ^{-1} \ll \varepsilon ^{-1/3}$ then yields the asymptotic condition $\bar t \ll \varepsilon ^{-1/6}$ .

3.7 Gaussians evolving on a quadratic potential

We already know from Proposition 3.2 that in the case of a quadratic potential V, variationally evolving complex Gaussians are exact solutions to the Schrödinger equation (3.1). We are now in a position to give the precise form and the equations of motion of the parameters of such a Gaussian solution.

Proposition 3.18 (Gaussian solutions for a quadratic potential)

Let V be a quadratic potential, let $(q(t),p(t),Q(t),P(t))$ be a solution of the classical equations of motion and their linearization,

(3.12) $$ \begin{align} \dot q=p, \ \ \dot p=-\nabla V(q) \quad\text{and}\quad \dot Q=P,\ \ \dot P=-\nabla^2V(q)Q, \end{align} $$

and let

$$ \begin{align*} S(t)= \int_0^t \biggl(\displaystyle \frac12|p(s)|^2-V(q(s))\biggr) \,{\mathrm d}s \end{align*} $$

be the corresponding classical action. Let $Q(0)$ and $P(0)$ satisfy the symplecticity relations (3.10). Then the complex Gaussian

(3.13) $$ \begin{align} &\mathrm{e}^{\mathrm{i} S(t)/\varepsilon} \,(\pi\varepsilon)^{-d/4} (\det Q(t))^{-1/2} \\ &\quad\, \times \exp\biggl( \dfrac {\mathrm{i}}{2\varepsilon} (x-q(t))^T P(t)Q(t)^{-1} (x-q(t)) + \dfrac{\mathrm{i}}{\varepsilon} p(t)^T(x-q(t)) \biggr) \notag \end{align} $$

is a solution to the time-dependent Schrödinger equation (3.1) of unit $L^2$ -norm. Moreover, the matrices $Q(t)$ and $P(t)$ satisfy the symplecticity relations (3.10) for all t.

The branch of the square root of $\det Q(t)$ is to be chosen such that this becomes a continuous function of t.

Proof. Since a variationally evolving Gaussian is an exact solution in the case of a quadratic potential, it suffices to show that a variational Gaussian takes the stated form. We have the quadratic polynomial

$$ \begin{align*} V(x) = V(q) + \nabla V(q)^T (x-q) + \displaystyle \frac{1}{2} (x-q)^T \nabla^2V(q)(x-q). \end{align*} $$

Since for a normalized Gaussian $u\in {\mathcal M}$ we have $P_uVu=Vu$ , Proposition 3.14 shows that

$$ \begin{align*} V(q)&= \langle V \rangle _u-\displaystyle \frac\varepsilon4\mathrm{tr}((\operatorname{Im } C)^{-1}\langle\nabla_x^2 V\rangle_u), \\ \nabla V(q) &= \langle \nabla V \rangle_u, \quad \nabla^2 V(q) = \langle \nabla^2 V \rangle_u. \end{align*} $$

Theorem 3.17 therefore shows that the differential equations for the parameters become (3.12) together with

$$ \begin{align*} \dot\zeta = \displaystyle \frac{1}{2} |p|^2 - V(q) + \displaystyle \frac{1}{2} \mathrm{i} \varepsilon\, \mathrm{tr}(C )\quad\text{with } C=PQ^{-1}, \end{align*} $$

and that $Q(t),P(t)$ satisfy (3.10) for all t. The Gaussian (3.13) is of unit norm by Lemma 3.7, while the relation $\operatorname{Im } C=(QQ^*)^{-1}$ follows from Lemma 3.16. So it remains to show that

$$ \begin{align*} \mathrm{e}^{\mathrm{i} \zeta(t)/\varepsilon} = \mathrm{e}^{\mathrm{i} S(t)/\varepsilon} \, (\pi\varepsilon)^{-d/4} (\det Q(t))^{-1/2} \,\mathrm{e}^{\mathrm{i} \phi} \end{align*} $$

for some real phase $\phi $ that is independent of t. As we already know that this holds at $t=0$ , taking the logarithm and differentiating with respect to t shows that this is satisfied if, for all t,

$$ \begin{align*} \dfrac{\mathrm{i}}{\varepsilon} \dot\zeta(t)=\dfrac{\mathrm{i}}{\varepsilon} \dot S(t) - \displaystyle \frac{1}{2} \dfrac{\mathrm{d}}{\mathrm{d} t} \log \det Q(t) . \end{align*} $$

This holds true because of the differential equation for $\zeta $ and because

$$ \begin{align*} \dfrac{\mathrm{d}}{\mathrm{d} t} \log \det Q = \mathrm{tr}(\dot Q Q^{-1}) = \mathrm{tr}(PQ^{-1}) = \mathrm{tr}(C). \end{align*} $$

This yields the stated result. □

We remark that the result and its proof extend in a straightforward way to time-dependent quadratic potentials $V(\cdot ,t)$ . Replacing the averages with point evaluations in solving the equations of motion of Theorem 3.17, i.e. instead solving the differential equations (3.12) for a general smooth potential V, therefore corresponds to solving the Schrödinger equation with a locally quadratic approximation to the potential at the current position $q(t)$ .

3.8 Notes

The time-dependent variational principle was first used by Dirac (Reference Dirac1930) in the context of what is now called the time-dependent Hartree–Fock method. This was then taken up in the book by Frenkel (Reference Frenkel1934). We refer to Kramer and Saraceno (Reference Kramer and Saraceno1981) and Lubich (Reference Lubich2008) for dynamic, geometric and approximation aspects of the Dirac–Frenkel time-dependent variational principle.

Variational Gaussian wave packets became prominent through the work of Heller (Reference Heller1976b ). In the chemical literature, simplified variants of variational Gaussian wave packets are often used, in particular with a diagonal width matrix (Hartree Gaussians: see Heller Reference Heller1976b ) or with a fixed width matrix (frozen Gaussians: see Heller Reference Heller1981). Variational linear combinations of Gaussian wave packets are also used successfully; see Richings et al. (Reference Richings, Polyak, Spinlove, Worth, Burghardt and Lasorne2015) and references therein. The book chapter by Vaníček and Begušić (Reference Vaníček, Begušić, Quack and Marquardt2020) reviews the application of various non-variational thawed Gaussian approximations, in particular of the single-Hessian thawed approximation that is energy-conserving despite being non-variational; see also Begušić, Cordova and Vaníček (Reference Begušić, Cordova and Vaníček2019).

The conservation of energy of a variational Gaussian wave packet generally follows from the time-dependent variational principle, which yields a Hamiltonian system on the approximation manifold. The conservation of norm and linear and angular momentum are instances of the general result that the average of a self-adjoint operator is conserved in the variational approximation if the operator commutes with the Hamiltonian and is such that it maps every point on the approximation manifold into its tangent space. These aspects are discussed in Lubich (Reference Lubich2008, Section II.1).

Error bounds for Gaussian wave packets were first given by Hagedorn (Reference Hagedorn1980) in a non-variational setting (with point evaluations instead of averages of the potential and its gradient and Hessian). The $O(\sqrt \varepsilon )$ $L^2$ -error bound of Theorem 3.5 for variational Gaussian wave packets was given by Lubich (Reference Lubich2008, Section II.4), whereas the $O(\varepsilon )$ error bound for the averages of observables appears to be new. We note that these error bounds are not valid for Gaussians with width matrices that are restricted to be diagonal or fixed.

The equations of motion for the parameters of a multi-variate variational Gaussian wave packet, as given in Theorem 3.11, were first derived by Coalson and Karplus (Reference Coalson and Karplus1990). The differential equations for the parameters have a non-canonical Hamiltonian structure that was studied by Faou and Lubich (Reference Faou and Lubich2006) and, with different geometric interpretations, by Ohsawa and Leok (Reference Ohsawa and Leok2013); see also Ohsawa (Reference Ohsawa2015b ) for further geometrical aspects of variational Gaussian wave packets and their use in obtaining conservation laws from symmetries.

The relationship between the Ehrenfest time and the classical Lyapunov exponent was studied in a more general setting by Combescure and Robert (Reference Combescure and Robert1997), using different arguments from those given here; see e.g. Bambusi, Graffi and Paul (Reference Bambusi, Graffi and Paul1999) and Hagedorn and Joye (Reference Hagedorn and Joye2000) for further developments.

The factorization (3.8) and its use with evolving Gaussian wave packets were introduced by Hagedorn (Reference Hagedorn1980), where Proposition 3.18 was also shown. In our text we use the notation of Lubich (Reference Lubich2008, Chapter V), working with the matrices $(Q,P)$ which correspond to $(A,\mathrm {i} B)$ in Hagedorn’s work, because of their relation to symplectic matrices and because their equations of motion are close to linearized classical equations of motion. The set of complex symmetric matrices with positive definite imaginary part is known as the Siegel half-space, after Siegel (Reference Siegel1943), who studied generalized Möbius transformations of this set based on symplectic matrices. The geometry behind the factorization (3.8) and its relation to Siegel’s results and to Marsden–Weinstein reduction was explored by Ohsawa (Reference Ohsawa2015a ).

4.1 Construction of Hagedorn functions via ladder operators

We start from a normalized d-dimensional complex Gaussian in Hagedorn’s parametrization (see Sections 3.5 and 3.7):

(4.1) $$ \begin{align} &\varphi_0^\varepsilon[q,p,Q,P](x) \\[6pt] &\quad = (\pi\varepsilon)^{-d/4} (\det Q)^{-1/2} \exp\biggl( \dfrac {\mathrm{i}}{2\varepsilon} (x-q)^T PQ^{-1} (x-q) + \dfrac{\mathrm{i}}{\varepsilon} p^T(x-q) \biggr) \notag \end{align} $$

with $q,p\in {\mathbb R}^d$ and matrices $Q,P\in {\mathbb C}^{d\times d}$ that satisfy the symplecticity relations (3.10). This Gaussian is of unit $L^2$ -norm by Lemmas 3.7 and 3.16. The branch of the square root is chosen appropriately, such that $(\det Q(t))^{-1/2}$ is a continuous function of t for a continuous family of invertible matrices $Q(t)$ . We write $\varphi _0$ or $\varphi _0^\varepsilon $ instead of $\varphi _0^\varepsilon [q,p,Q,P]$ when the parameters are clear from the context.

We denote the position and momentum operators by $\widehat q=(\widehat q_j)_{j=1}^d$ and $\widehat p=(\widehat p_j)_{j=1}^d$ , respectively:

$$ \begin{align*} \widehat q \psi = x \psi,\quad \widehat p \psi = -\mathrm{i} \varepsilon \nabla_x\psi. \end{align*} $$

The hats are added to avoid confusion with the Gaussian parameters q and p. The commutator relation

(4.2) $$ \begin{align} \dfrac1{\mathrm{i} \varepsilon}\, [\widehat q_j, \widehat p_k] = \delta_{jk} \end{align} $$

is fundamental in quantum mechanics and, in particular, in the construction of the harmonic oscillator eigenfunctions (the Hermite functions) using Dirac’s ladder. Extending that construction, we will consider Hagedorn’s parameter-dependent ladder operators $A=(A_j)_{j=1}^d$ and $A^\dagger =(A_j^\dagger )_{j=1}^d$ , defined as

(4.3) $$ \begin{align} \begin{aligned} A\phantom{^\dagger} & = A[q,p,Q,P] = -\dfrac {\mathrm{i}}{\sqrt{2\varepsilon}} ( P^T(\widehat q - q) - Q^T(\widehat p - p) ){,} \\ A^\dagger & = A^\dagger[q,p,Q,P] = \dfrac {\mathrm{i}}{\sqrt{2\varepsilon}} ( P^*(\widehat q - q) - Q^*(\widehat p - p) ) . \end{aligned} \end{align} $$

We note that for $d=1$ , $\varepsilon =1$ , $q=0$ , $p=0$ , $Q=1$ , $P={\mathrm {i}}$ , these operators reduce to Dirac’s ladder operators ${{{{(\widehat q+\mathrm {i} \widehat p)}}/{{\sqrt {2}}}}}$ and ${{{{(\widehat q-\mathrm {i} \widehat p)}}/{{\sqrt {2}}}}}$ , respectively. The key properties of those operators extend as follows.

Lemma 4.1 (commutator relations)

If Q and P satisfy the symplecticity relations (3.10), then we have the commutator identities

(4.4) $$ \begin{align} [A_j, A_k^\dagger] = \delta_{jk} . \end{align} $$

Moreover, $A_j^\dagger $ is adjoint to $A_j$ on the Schwartz space ${\mathcal S}$ of functions on ${\mathbb R}^d$ that are infinitely differentiable and decay faster than any negative power:

(4.5) $$ \begin{align} \langle A_j^\dagger \varphi \mid \psi \rangle = \langle \varphi \mid A_j \psi \rangle \quad \text{for all } \varphi,\psi\in {\mathcal S} . \end{align} $$

Proof. With $Q=(Q_{jk})$ and $P=(P_{jk})$ , we have (with $q=p=0$ for simplicity)

$$ \begin{align*} [ A_j, A_k^\dagger] = \dfrac1{2\varepsilon} \Biggl[ \sum_{\ell=1}^d ( P_{\ell j} \widehat q_\ell - Q_{\ell j} \widehat p_\ell ) \, , \, \sum_{m=1}^d ( \overline P_{m k} \widehat q_m - \overline Q_{m k} \widehat p_m ) \Biggr] . \end{align*} $$

By the canonical commutator identities (4.2), this simplifies to

$$ \begin{align*} [ A_j, A_k^\dagger] = \dfrac {\mathrm{i}}2 \sum_{\ell=1}^d ( -P_{\ell j}\overline Q_{\ell k} + Q_{\ell j} \overline P_{\ell k} ) = \dfrac {\mathrm{i}}2 ( -Q^*P + P^*Q )_{k,j}, \end{align*} $$

and by the symplecticity relation (3.10), this equals $\delta _{jk}$ .

The relation (4.5) follows directly from the fact that $\widehat q_\ell $ and $\widehat p_\ell $ are symmetric operators on ${\mathcal S}$ .□

Proof. If $\varphi \in {\mathcal S}$ is in the null-space of A, then it must satisfy the linear system of partial differential equations

$$ \begin{align*} -\mathrm{i} \varepsilon\nabla\varphi(x) - p\varphi(x) = C(x-q)\varphi(x) \end{align*} $$

with the complex symmetric matrix $C=PQ^{-1}$ . Multiples of $\varphi _0$ are the only non-trivial solutions of this equation, since substituting $\varphi =\rho \varphi _0$ yields the equation $\nabla \rho =0$ , so that $\rho $ is constant. □

With the properties established in these two lemmas, we can now construct eigenfunctions of the operators $A_jA_j^\dagger $ to eigenvalues $1,2,3,\ldots $ by the same arguments that are familiar from Dirac’s construction of the eigenfunctions of the harmonic oscillator

$$ \begin{align*} \displaystyle \frac12(\widehat q^2 + \widehat p^2)= \displaystyle \frac1{\sqrt{2}}(\widehat q+\mathrm{i} \widehat p)\,\displaystyle \frac1{\sqrt{2}}(\widehat q-\mathrm{i} \widehat p)-\displaystyle \frac12. \end{align*} $$

Let $k=(k_1,\ldots ,k_d)$ be a multi-index with non-negative integers $k_j$ , and let $\langle j \rangle = (0 ,\ldots , 1 ,\ldots , 0)$ denote the jth d-dimensional unit vector.

With Lemmas 4.1 and 4.2 it follows that

$$ \begin{align*} A_jA_j^\dagger\varphi_0 = A_j^\dagger A_j\varphi_0 + \varphi_0=\varphi_0 , \end{align*} $$

and hence $\varphi _0$ is an eigenfunction of $A_jA_j^\dagger $ to the eigenvalue $1$ . Applying the operator $A_j^\dagger $ to both sides of this equation, we see that $\phi _{\langle j \rangle } :=A_j^\dagger \varphi _0$ is an eigenfunction of $A_j^\dagger A_j$ to the eigenvalue $1$ .

On the other hand, we further have by Lemmas 4.1 and 4.2 that

$$ \begin{align*} A_lA_j^\dagger\varphi_0 = A_j^\dagger A_l\varphi_0 =0 \quad \text{for }l\ne j, \end{align*} $$

and applying the operator $A_l^\dagger $ to both sides of this equation, we see that $\phi _{\langle j \rangle } =A_j^\dagger \varphi _0$ is an eigenfunction of $A_l^\dagger A_l$ to the eigenvalue $0$ .

We continue in this way to construct recursively

(4.6) $$ \begin{align} \phi_{k+\langle j \rangle}=A_j^\dagger\phi_k \end{align} $$

for every multi-index $k=(k_1,\ldots ,k_d)$ with non-negative integers $k_l$ . We thus obtain eigenfunctions $\phi _k$ to $A_l^\dagger A_l$ with eigenvalue $k_l$ and to $A_lA_l^\dagger $ with eigenvalue $k_l+1$ . These eigenfunctions are not yet normalized. To achieve this, we note that by Lemma 4.1

$$ \begin{align*} \| A_j^\dagger\phi_k \|^2 = \langle A_j^\dagger\phi_k \mid A_j^\dagger\phi_k \rangle = \langle \phi_k \mid A_jA_j^\dagger\phi_k \rangle = (k_j+1)\,\|\phi_k\|^2 . \end{align*} $$

We therefore obtain normalized eigenfunctions to $A_lA_l^\dagger $ and $A_l^\dagger A_l$ by defining the functions $\varphi _k=\varphi _k^\varepsilon [q,p,Q,P]$ recursively by

(4.7) $$ \begin{align} \varphi_{k+\langle j \rangle} = \dfrac1{\sqrt{k_j+1}}\, A_j^\dagger\varphi_k. \end{align} $$

Since

$$ \begin{align*} A_j\varphi_{k+\langle j \rangle}= \displaystyle \frac1{\sqrt{k_j+1}}\, A_jA_j^\dagger\varphi_k = \sqrt{k_j+1}\, \varphi_k, \end{align*} $$

we also have (replacing $k_j+1$ with $k_j$ )

(4.8) $$ \begin{align} \varphi_{k-\langle j \rangle} = \dfrac1{\sqrt{k_j}}\, A_j\varphi_{k}, \end{align} $$

so that $A_j^\dagger $ and $A_j$ can be viewed as raising and lowering operators, respectively, in the jth component of the multi-index. Multiplying (4.7) by $\sqrt {k_j+1}$ and (4.8) by $\sqrt {k_j}$ , summing the resulting formulas and using the definitions of $A_j$ and $A_j^\dagger $ and the fact that $QQ^*$ is a real matrix by Lemma 3.16, we obtain the three-term recurrence relation

(4.9) $$ \begin{align} Q\bigl(\sqrt{k_j+1}\,\varphi_{k+{{\langle j \rangle}}}(x)\bigr)_{j=1}^d = \sqrt{\dfrac2\varepsilon}\, (x-q)\varphi_{k}(x) - \overline Q \bigl(\sqrt{k_j}\,\varphi_{k-{{\langle j \rangle}}}(x)\bigr)_{j=1}^d . \end{align} $$

This permits us to evaluate $\varphi _k(x)$ at any given point x. It also shows that $\varphi _k$ is the product of a polynomial of degree $|k|=k_1+\cdots +k_d$ with the Gaussian $\varphi _0$ . Writing this property as

(4.10) $$ \begin{align} \varphi_k(x) = \dfrac{1}{\sqrt{2^{|k|}k!}} \ p_k\biggl(\displaystyle \frac{1}{\sqrt\varepsilon}Q^{-1}(x-q)\biggr)\ \varphi_0(x), \end{align} $$

we obtain the polynomial three-term recurrence relation

(4.11) $$ \begin{align} (p_{k+{{\langle j \rangle}}}(x))_{j=1}^d = 2x \, p_{k}(x) - 2Q^{-1}\,\overline Q (k_j\,p_{k-{{\langle j \rangle}}}(x))_{j=1}^d. \end{align} $$

This shows that if the matrix $Q^{-1}\,\overline Q$ is diagonal, then each $p_k$ is the tensor product of appropriately scaled univariate Hermite polynomials. Indeed, if $Q^{-1}\,\overline Q = \operatorname {diag}(\lambda _1,\ldots ,\lambda _d)$ , then

$$ \begin{align*} p_k(x) = \prod_{j=1}^d \lambda_j^{k_j/2} \ H_{k_j} \biggl(\dfrac{x_j}{\sqrt{\lambda_j}}\biggr) \end{align*} $$

where $(H_n)_{n\ge 0}$ are the Hermite polynomials defined by the three-term recurrence relation $H_{n+1}(y) = 2yH_n(y) - 2n H_{n-1}(y),H_0(y) = 1$ . If the matrix $Q^{-1}\,\overline Q$ is not diagonal, then the polynomials $p_k$ do not have tensor product form.

The Hagedorn functions $\varphi _k$ are eigenfunctions of the symmetric operators $A_j A_j^\dagger $ for each j. Therefore they are orthogonal,

$$ \begin{align*} (k_j+1) \langle \varphi_k,\varphi_\ell \rangle = \langle A_j A_j^\dagger \varphi_k,\varphi_\ell \rangle= \langle \varphi_k,A_j A_j^\dagger\varphi_\ell \rangle= (\ell_j+1) \langle \varphi_k,\varphi_\ell \rangle, \end{align*} $$

so that $ \langle \varphi _k,\varphi _\ell \rangle =0$ for $k\ne \ell $ . The functions $\varphi _k$ even form a complete orthonormal system, as can be shown by extending the arguments in the proof of completeness of the Hermite functions. We do not present this proof here, as our interest in the following will be in the evolution of localized wave packets with a moderate polynomial degree rather than in the expansion of general functions in $L^2({\mathbb R}^d)$ in the orthonormal basis of Hagedorn functions.

We summarize the above construction as follows.

4.2 Evolution of Hagedorn wave packets for quadratic Hamiltonians

We now consider the Hamiltonian

$$ \begin{align*} H=-\displaystyle \frac{\varepsilon^2}{2}\Delta + V= \displaystyle \frac12 \widehat p\cdot\widehat p +V \end{align*} $$

with a quadratic potential V, and take a Hagedorn function as initial data for the time-dependent Schrödinger equation (3.1). There is the following remarkable result.

Theorem 4.4 (Hagedorn functions on a quadratic potential)

Let V be a quadratic potential, let $(q(t),p(t),Q(t),P(t))$ be a solution of the classical equations of motion and their linearization,

(4.12) $$ \begin{align} \dot q=p, \ \ \dot p=-\nabla V(q) \quad\text{and}\quad \dot Q=P,\ \ \dot P=-\nabla^2V(q)Q, \end{align} $$

and let

$$ \begin{align*} S(t)= \int_0^t \biggl(\displaystyle \frac12|p(s)|^2-V(q(s))\biggr) \,{\mathrm d}s \end{align*} $$

be the corresponding classical action. Assume that $Q(0)$ and $P(0)$ satisfy the symplecticity relations (3.10). Then, for every multi-index k,

$$ \begin{align*} \mathrm{e}^{\mathrm{i} S(t)/\varepsilon} \,\varphi_k^\varepsilon[q(t),p(t),Q(t),P(t)](x) \end{align*} $$

is a solution of the time-dependent Schrödinger equation (3.1). Moreover, the symplecticity relations (3.10) are satisfied by $Q(t)$ and $P(t)$ for all t.

As a preparation for the proof, we first consider the evolution of the parameter-dependent ladder operators under the quadratic Hamiltonian. Along a solution $(q(t),p(t),Q(t),P(t))$ of the classical equations (4.12) we consider the time-dependent operators

$$ \begin{align*} A_j(t) = A_j[q(t),p(t),Q(t),P(t)],\quad A_j^\dagger(t) = A_j^\dagger[q(t),p(t),Q(t),P(t)] . \end{align*} $$

These operators evolve according to the following equations.

Lemma 4.5 (ladder evolution for quadratic potentials)

Let

$$ \begin{align*} H=-\displaystyle \frac{\varepsilon^2}{2}\Delta + V \end{align*} $$

be the Hamiltonian of (3.1) with a quadratic potential V. Then we have

$$ \begin{align*} \dot A_j = \dfrac1{\mathrm{i} \varepsilon}[A_j,H] ,\quad \dot A_j^\dagger = -\dfrac1{\mathrm{i} \varepsilon}[A_j^\dagger ,H] . \end{align*} $$

Proof. With (4.12) we obtain for $A(t)=(A_j(t))$

$$ \begin{align*} \dot A = \dfrac {\mathrm{i}}{\sqrt{2\varepsilon}} ( Q^T \nabla V(\widehat q) + P^T \widehat p ). \end{align*} $$

The same expression is obtained for ${{{{[A,H]}}/{{(\mathrm {i} \varepsilon )}}}}$ with the commutator relations (4.2), which further yield

$$ \begin{align*} \displaystyle \frac1{\mathrm{i} \varepsilon}[\widehat q_j, \widehat p_k^2] = \delta_{jk}\cdot 2 \widehat p_k \quad\text{and}\quad \displaystyle \frac1{\mathrm{i} \varepsilon}[\widehat q_j^2, \widehat p_k] = \delta_{jk}\cdot 2 \widehat q_j. \end{align*} $$

The result for $A^\dagger $ is obtained by taking complex conjugates. □

Proof of Theorem 4.4

By Proposition 3.18, the result holds true for ${k=0}$ . In view of the construction of the functions $\varphi _k$ by (4.7), the result for general k follows by induction if we can show that with a solution $\psi (\cdot ,t)$ , also $A_j^\dagger (t)\psi (\cdot ,t)$ is a solution to the Schrödinger equation (3.1). This follows from

$$ \begin{align*} \mathrm{i} \varepsilon \dfrac\partial{\partial t}(A_j^\dagger\psi)= \mathrm{i} \varepsilon {\dot A}_j^\dagger\psi + A_j^\dagger H\psi = ( \mathrm{i} \varepsilon {\dot A}_j^\dagger\psi + [A_j^\dagger,H]\psi ) + H A_j^\dagger \psi , \end{align*} $$

because $\mathrm {i} \varepsilon {\dot A}_j^\dagger \psi + [A_j^\dagger ,H]\psi =0$ by Lemma 4.5. The preservation of the symplecticity relations (3.10) is already known from Proposition 3.18.□

Like Proposition 3.18, Theorem 4.4 also remains valid, with the same proof, for a time-dependent quadratic potential $V(\cdot ,t)$ .

4.3 A Galerkin method for non-quadratic potentials

For the time-dependent Schrödinger equation (3.1) with a non-quadratic potential V, the wave function will now be approximated by a finite linear combination of Hagedorn functions whose parameters evolve according to the classical equations of motion. The coefficients of the linear combination are determined by a Galerkin condition on the time-varying approximation space that is spanned by a basis of finitely many evolving Hagedorn functions.

We search for an approximation to the wave function of the form

(4.13) $$ \begin{align} \psi_{\mathcal K}(x,t)=\mathrm{e}^{\mathrm{i} S(t)/\varepsilon} \sum_{k\in{\mathcal K}} c_k(t)\, \varphi_k(x,t) {,} \end{align} $$

where we abbreviate

$$ \begin{align*} \varphi_k(x,t)=\varphi_k^\varepsilon[q(t),p(t),Q(t),P(t)](x) \end{align*} $$

for a solution $(q(t),p(t),Q(t),P(t))$ to the classical equations (4.12), and where

$$ \begin{align*} S(t)= \int_0^t \biggl(\displaystyle \frac12|p(s)|^2-V(q(s))\biggr) \,{\mathrm d}s \end{align*} $$

is the corresponding classical action. The finite multi-index set ${\mathcal K}$ may, for example, be a cube

$$ \begin{align*} \{k\in {\mathbb N}^d\colon |k_j|\le K\} \end{align*} $$

or a simplex

$$ \begin{align*} \biggl\{ k \in {\mathbb N}^d\colon |k|= \sum_{j=1}^d k_j \le K\biggr\} \end{align*} $$

or a hyperbolically reduced set

$$ \begin{align*} \biggl\{ {k \in {\mathbb N}^d}\colon \prod_{j=1}^d (1+k_j) \le K\biggr\}. \end{align*} $$

The initial data will be assumed to be a linear combination over a smaller set ${\mathcal K}_0\subset {\mathcal K}$ .

The coefficients $c_k(t)$ are determined by the Galerkin condition

(4.14) $$ \begin{align} \biggl\langle \varphi_\ell(t) \mid \mathrm{i} \varepsilon\dfrac{\partial \psi_{\mathcal K}}{\partial t}(t) -H\psi_{\mathcal K}(t) \biggr\rangle = 0 \quad \text{for each }\ell\in{\mathcal K}\text{ and for all }t. \end{align} $$

If we write the potential as

$$ \begin{align*} V=U_{q(t)} + W_{q(t)} \end{align*} $$

with the quadratic Taylor polynomial of V at the position q,

$$ \begin{align*} U_q(x) =V(q) + \nabla V(q) ^T(x-q) +\displaystyle \frac{1}{2} (x-q)^T \nabla^2 V(q) (x-q), \end{align*} $$

and with the non-quadratic remainder $W_q$ , then we have

(4.15) $$ \begin{align} &\mathrm{i} \varepsilon\dfrac\partial{\partial t}(c_k \,\mathrm{e}^{\mathrm{i} S/\varepsilon}\varphi_k) - H(c_k \,\mathrm{e}^{\mathrm{i} S/\varepsilon}\varphi_k) = \mathrm{i} \varepsilon \dot c_k \,\mathrm{e}^{\mathrm{i} S/\varepsilon}\varphi_k \\ \notag &\quad\, + \ c_k \biggl( \mathrm{i} \varepsilon\dfrac\partial{\partial t}(\mathrm{e}^{\mathrm{i} S/\varepsilon}\varphi_k) + \dfrac{\varepsilon^2}{2}\Delta (\mathrm{e}^{\mathrm{i} S/\varepsilon}\varphi_k) - U_q \,\mathrm{e}^{\mathrm{i} S/\varepsilon}\varphi_k \biggr) - c_k \,\mathrm{e}^{\mathrm{i} S/\varepsilon} W_q \varphi_k. \end{align} $$

The term in big brackets vanishes by the version of Theorem 4.4 for time-dependent quadratic potentials, since just the quadratic part $U_q$ of the potential V determines $q(t),p(t),Q(t),P(t),S(t)$ . The Galerkin condition (4.14) then gives us the following differential equations for the coefficients $c(t)=(c_k(t))_{k\in {\mathcal K}}$ :

(4.16) $$ \begin{align} \mathrm{i} \varepsilon \dot c(t) = G(t)c(t) \quad\text{with } G(t)= ( \langle \varphi_\ell(\cdot,t) \mid W_{q(t)}\varphi_k(\cdot,t) \rangle )_{\ell,k\in{\mathcal K}} . \end{align} $$

We now easily observe norm conservation of the approximation.

With this approximation method, the differential equations (4.12) for the parameters $(q(t),p(t),Q(t),P(t))$ of the Hagedorn functions are obviously independent of the coefficients $c_k(t)$ , but the differential equations (4.16) depend on the parameters through the basis functions $\varphi _k(\cdot ,t)=\varphi _k^\varepsilon [q(t),p(t),Q(t),P(t)]$ . We remark that with a fully variational approach for parameters and coefficients such a separation would not be attained. However, we do not have exact energy conservation with the semi-variational approach considered here.

An important feature is that the coefficients are slowly varying,

$$ \begin{align*} \|\dot c(t)\| \le C\varepsilon^{1/2}, \end{align*} $$

which holds true because the norm of the Galerkin matrix $G(t)$ is bounded by $C\varepsilon ^{3/2}$ with a constant C that depends on the multi-index set ${\mathcal K}$ and on the width of the underlying Gaussian $\varphi _0(\cdot ,t)$ . This follows directly from Lemma 3.8. We even have the following stronger bounds for the entries of the matrix $G(t)$ .

Lemma 4.7 (bounds for the entries of the Galerkin matrix)

Let ${\varphi _k^\varepsilon =\varphi _k^\varepsilon [q,p,Q,P]}$ for $\varepsilon>0$ and $k\in {\mathbb N}^d$ denote the Hagedorn functions with parameters $(q,p,Q,P)$ , and let W be a real-valued function on ${\mathbb R}^d$ with arbitrarily many polynomially bounded derivatives such that W and its first two derivatives vanish at q. For $k,\ell \in {\mathbb N}^d$ we then have the bound

$$ \begin{align*} | \langle \varphi_\ell^\varepsilon \mid W \varphi_k^\varepsilon \rangle | \le c \, \varepsilon^{\mu/2} \quad\text{with } \mu= \max (|k-\ell|,3), \end{align*} $$

where c is independent of $\varepsilon $ but depends on $k,\ell $ , on the matrix 2-norm of Q and on bounds of derivatives of W.

Proof. We start from the Taylor expansion of W at q with integral remainder term, that is,

$$ \begin{align*} W(x) & =\sum_{3 \le |m| \le N} \dfrac1{|m|!} \,\partial^m W(q) (x-q)^m \\ &\quad\, + \int_0^1 \dfrac{(1-\theta)^N}{N!} \sum_{|m|=N+1} \partial^m W(q+\theta(x-q)) (x-q)^m\,{\mathrm d}\theta, \end{align*} $$

where for a multi-index $m=(m_1,\ldots ,m_d)$ , we denote the partial derivatives $\partial ^m W = \partial _1^{m_1}\cdots \partial _d^{m_d} W$ and the product $(x-q)^m=(x_1-q_1)^{m_1}\cdots (x_d-q_d)^{m_d}$ .

So we need to study the inner product $\langle \varphi _\ell ^\varepsilon \mid (x-q)^m \varphi _k^\varepsilon \rangle $ for $k,\ell ,m$ . By the normalization of the wave packets, we have

$$ \begin{align*} |\langle \varphi_\ell^\varepsilon \mid (x-q)^m \varphi_k^\varepsilon \rangle| \le \|(x-q)^m \varphi_k\|. \end{align*} $$

Next we write $\varphi _k^\varepsilon $ in the form (4.10), that is,

$$ \begin{align*} {\varphi_k^\varepsilon(x) = \dfrac{1}{\sqrt{2^{|k|}k!}} \, p_k\biggl(\displaystyle \frac{1}{\sqrt\varepsilon}Q^{-1}(x-q)\biggr) \varphi_0(x).} \end{align*} $$

Hence there exists a constant $\gamma>0$ that depends on k and $Q^{-1}\,\overline Q$ , such that

$$ \begin{align*} &\| (x-q)^m \varphi_k^\varepsilon \| \\ &\quad \le \gamma \biggl( \int_{{\mathbb R}^d} |x|^{2m} \ \biggl(1 + \displaystyle \frac{1}{\varepsilon}|Q^{-1}x|^2\biggr)^{k} \ \exp\biggl(-\displaystyle \frac1\varepsilon x^T (QQ^*)^{-1} x\biggr) {\,\mathrm{d}} x \biggr)^{1/2}. \end{align*} $$

By Lemma 3.8,

$$ \begin{align*} \langle \varphi_\ell^\varepsilon \mid (x-q)^m \varphi_k^\varepsilon\rangle = O(\varepsilon^{|m|/2}). \end{align*} $$

Here we note further that the smallest eigenvalue of the width matrix $\operatorname{Im } C=(QQ^*)^{-1}$ equals the square of the largest singular value of Q, which is the matrix 2-norm of Q.

Moreover, in view of the orthogonality of the Hagedorn functions and the three-term recurrence relation (4.9), the inner product $\langle \varphi _\ell ^\varepsilon \mid (x-q)^m \varphi _k^\varepsilon \rangle $ can be non-zero only if there exist standard unit vectors $\langle j_1 \rangle ,\ldots ,\langle j_{|m|} \rangle $ of ${\mathbb R}^d$ such that, with appropriate signs,

$$ \begin{align*} k\pm \langle j_1 \rangle \pm \cdots \pm \langle j_{|m|} \rangle = \ell, \end{align*} $$

which requires $ |k-\ell | \le |m|. $ So we obtain

(4.17) $$ \begin{align} | \langle \varphi_\ell^\varepsilon \mid (x-q)^m \varphi_k^\varepsilon \rangle | = \begin{cases} 0 &\mbox{if } |m|<|k-\ell| {,}\\ O(\varepsilon^{|m|/2}) & \mbox{if } |m| \ge |k-\ell|. \end{cases} \end{align} $$

Using these estimates after inserting the Taylor expansion of W at q in the expression $ \langle \varphi _\ell ^\varepsilon \mid W \varphi _k^\varepsilon \rangle $ gives the result. □

4.4 Approximation of higher order in $\varepsilon $ by the Hagedorn–Galerkin method

Theorem 3.5 provided an $O(\varepsilon ^{1/2})$ error bound over time $O(1)$ for variational Gaussians. For any chosen multi-index set ${\mathcal K}$ , a bound of the same approximation order $1/2$ in $\varepsilon $ can also be obtained for the Hagedorn–Galerkin method of the previous subsection. However, as we will show next, an error bound of any prescribed order in $\varepsilon $ is obtained if the initial wave packet has non-vanishing coefficients $c_k(0)$ only for k in some fixed initial multi-index set ${\mathcal K}_0$ and the multi-index set ${\mathcal K}$ of the Galerkin method is substantially larger than ${\mathcal K}_0$ . We begin by studying the particular case where ${\mathcal K}_0=\{ k_0 \}$ has a single element $k_0\in {\mathbb N}^d$ . The case of a general finite initial set ${\mathcal K}_0$ then follows by linear superposition.

Theorem 4.8 (higher-order error bound)

Let $N\ge 2$ and suppose that for some $k_0\in {\mathbb N}^d$ the multi-index set ${\mathcal K}$ of the Galerkin approximation (4.13) is such that

(4.18) $$ \begin{align} {\{ k\in {\mathbb N}^d\colon |k-k_0| < 3N \} \subset {\mathcal K}.} \end{align} $$

Let $\psi (\cdot ,t)$ denote the solution of the Schrödinger equation (3.1) with a smooth, polynomially bounded potential V, for initial data given by a single Hagedorn function with multi-index $k_0$ :

$$ \begin{align*} \psi(\cdot,0) = \varphi_{k_0}^\varepsilon[q(0),p(0),Q(0),P(0)], \end{align*} $$

where $Q(0),P(0)$ satisfy the symplecticity relations (3.10). Let $\psi _{\mathcal K}(\cdot ,t)$ be the Galerkin approximation (4.13) with parameters $(q(t),p(t),Q(t),P(t))$ satisfying the classical equations (4.12) and with coefficients $c_k(t)$ for $k\in {\mathcal K}$ determined by (4.16), with initial data $\psi _{\mathcal K}(\cdot ,0)=\psi (\cdot ,0)$ , i.e. $c_k(0)=\delta _{k,k_0}$ . Then the error is bounded by

$$ \begin{align*} \| \psi_{\mathcal K}(\cdot,t) - \psi(\cdot,t) \| \le C \varepsilon^{N/2} \quad\text{for } 0\le t\le T, \end{align*} $$

where C is independent of $\varepsilon $ , but depends on $k_0$ , ${\mathcal K}$ , on $\max _{0\le t \le T}\|Q(t)\|_2$ and T.

Proof. The proof is in three parts. In the first part we derive bounds for the coefficients $c_k(t)$ that decay geometrically with powers of $\varepsilon $ for k moving away from $k_0$ . In the second part we use this to bound the defect that results from inserting the Galerkin approximation into the Schrödinger equation. Finally, we conclude from a small defect to a small error.

(a) We begin by showing that for $0\le t \le T$

(4.19) $$ \begin{align} c_k(t) =O( \varepsilon^{|k-k_0|/6}), \end{align} $$

where the constant symbolized by the O-notation depends on $k_0$ , ${\mathcal K}$ and T but is independent of $\varepsilon $ . Let us define, for some $\alpha \ge 0$ ,

$$\begin{align*}b_k(t) = c_k(t) / \varepsilon^{\alpha |k-k_0|}. \end{align*}$$

Then $b_k(0)=\delta _{k,k_0}$ , and since $c(t)$ solves (4.16), $b(t)=(b_k(t))_{k\in {\mathcal K}}$ is a solution to the differential equation

$$ \begin{align*} \dot b(t) = F(t) b(t) \quad\text{with } F(t)=(f_{\ell k}(t)) , \quad f_{\ell k}(t) = \dfrac 1{{\mathrm{i}} \varepsilon} \varepsilon^{-\alpha |\ell-k_0|} g_{\ell k}(t) \varepsilon^{\alpha|k-k_0|}. \end{align*} $$

We choose $\alpha $ maximal such that $f_{\ell k}(t)=O(1)$ for all $\ell ,k$ . By Lemma 4.7, $g_{\ell k}(t) = O(\varepsilon ^{\max (|k-\ell |,3)/2})$ , and by the triangle inequality we have $- |\ell -k_0| + |k-k_0| \ge - |k-\ell |$ . So we require, for arbitrary $n=|k-\ell |$ ,

$$ \begin{align*} -1 - \alpha n + \displaystyle \frac{1}{2}\max(n,3) \ge 0. \end{align*} $$

For $n=3$ this yields the condition $\alpha \le 1/6$ , and it is then seen that the inequality is satisfied for all $n\in {\mathbb N}$ with $\alpha =1/6$ . For this $\alpha $ , we thus have $b(t)=O(1)$ , which proves (4.19).

(b) We show that the defect obtained on inserting the approximation $\psi _{\mathcal K}$ into the Schrödinger equation,

$$ \begin{align*} d(t) := \dfrac{\partial \psi_{\mathcal K}}{\partial t}(t) - \dfrac 1 {\mathrm{i} \varepsilon} H \psi_{\mathcal K}(t) , \end{align*} $$

is bounded, uniformly for $0\le t \le T$ , by

(4.20) $$ \begin{align} \| d(t) \| = O(\varepsilon^{N/2}). \end{align} $$

By (4.15) and (4.16), the defect equals

$$ \begin{align*} d(t) = \dfrac 1 {\mathrm{i} \varepsilon} P_{\mathcal K}(t)^\perp W_{q(t)} \psi_{\mathcal K}(t), \end{align*} $$

where $P_{\mathcal K}(t)^\perp ={\mathrm {Id}} - P_{\mathcal K}(t)$ and $P_{\mathcal K}(t)$ is the orthogonal projection from $L^2({\mathbb R}^d)$ onto the space spanned by

$$ \begin{align*} \varphi_k(t)=\varphi_k^\varepsilon[q(t),p(t),Q(t),P(t)]\quad\text{for } k\in{\mathcal K}. \end{align*} $$

We split $W_{q(t)}$ into its $(N+1)$ st degree Taylor polynomial $W_{N+1}(\cdot ,t)$ at $q(t)$ and the remainder $R_{N+1}(\cdot ,t)$ . With the integral formula for the remainder term, Lemma 3.8 shows that

(4.21) $$ \begin{align} \biggl\| \displaystyle \frac1\varepsilon R_{N+1}(t) \psi_{\mathcal K}(t)\biggr\| = O(\varepsilon^{(N+2)/2-1}) = O(\varepsilon^{N/2}). \end{align} $$

On the other hand, $W_{N+1}(t)\psi _{\mathcal K}(t)$ is a linear combination of Hagedorn functions $\varphi _\ell (t)$ with $|\ell | \le K +N+1$ , where $K=\max _{k\in {\mathcal K}}|k|$ . We thus have the finite sum

$$ \begin{align*} \dfrac 1 {\mathrm{i} \varepsilon} P_{\mathcal K}(t)^\perp W_{N+1}(t) \psi_{\mathcal K}(t) = \dfrac 1 {\mathrm{i} \varepsilon}\, \mathrm{e}^{\mathrm{i} S(t)/\varepsilon} \sum_{\ell \notin {\mathcal K}} \sum_{k\in{\mathcal K}} c_k(t) \langle \varphi_\ell(t), W_{N+1}(t) \varphi_k(t) \rangle, \end{align*} $$

where the first sum is over all $\ell \in {\mathbb N}^d$ with $|\ell | \le K +N+1$ and $\ell \notin {\mathcal K}$ . By Lemma 4.7,

$$ \begin{align*} \langle \varphi_\ell(t), W_{N+1}(t) \varphi_k (t)\rangle = O(\varepsilon^{\max(|k-\ell|,3)/2}), \end{align*} $$

and together with (4.19) this gives us

$$ \begin{align*} \dfrac 1 {\mathrm{i} \varepsilon} \,c_k(t) \langle \varphi_\ell(t), W_{N+1}(t) \varphi_k(t) \rangle = O(\varepsilon^{-1+|k-k_0|/6 + \max(|k-\ell|,3)/2}). \end{align*} $$

By the triangle inequality, we have $|k-k_0|\ge |\ell -k_0|-|k-\ell |$ , and by condition (4.18) we have $|\ell -k_0|\ge 3N$ for $\ell \notin {\mathcal K}$ . Therefore

$$ \begin{align*} -1+\dfrac{|k-k_0|}6 + \dfrac{\max(|k-\ell|,3)}2 \ge -1 + \dfrac{N}2 -\dfrac{|k-\ell|}6+ \dfrac{\max(|k-\ell|,3)}2 \ge \dfrac N2, \end{align*} $$

so that

$$ \begin{align*} \dfrac 1 {\mathrm{i} \varepsilon} \,c_k(t) \langle \varphi_\ell(t), W_N(t) \varphi_k(t) \rangle = O(\varepsilon^{N/2}), \end{align*} $$

which together with (4.21) yields (4.20).

(c) Finally, with Lemma 3.6 we conclude from a small defect to a small error: for $0\le t \le T$ ,

$$ \begin{align*} \| \psi_{\mathcal K}(t) - \psi(t) \| \le \int_0^t \| d(s)\| {\,\mathrm{d}} s = O(\varepsilon^{N/2}). \end{align*} $$

This gives the stated result. □

Remark 4.1. For $N=2$ , condition (4.18) can be replaced by the weaker condition

$$ \begin{align*} {\{ k\in {\mathbb N}^d\colon |k-k_0| \le 3 \} \subset {\mathcal K}.} \end{align*} $$

We first obtain $c_k(t)=O(\varepsilon ^{1/2})$ for $k\ne k_0$ . As in (b) above, using this bound, the above condition on ${\mathcal K}$ and Lemma 4.7 for $k=k_0$ and $|\ell -k_0|\ge 4$ , we obtain the defect bound $\|d(t)\|=O(\varepsilon )$ , which yields the $O(\varepsilon )$ error bound.

By the above proof and by linearity, we immediately get the following more general result where the coefficients of the initial data decay sufficiently fast away from some finite multi-index set ${\mathcal K}_0$ .

Theorem 4.9 (higher-order error bound)

Let $N\ge 2$ and suppose that for some finite multi-index set ${\mathcal K}_0\subset {\mathbb N}^d$ , the multi-index set ${\mathcal K}$ of the Galerkin approximation (4.13) is such that, with the distance $\mathrm {dist}(k,{\mathcal K}_0) := \min _{k_0\in {\mathcal K}_0} {|k-k_0|}$ ,

$$ \begin{align*} \{ k\in {\mathbb N}^d\colon \mathrm{dist}(k,{\mathcal K}_0) < 3N \} \subset {\mathcal K}. \end{align*} $$

Let $\psi (\cdot ,t)$ denote the solution of the Schrödinger equation (3.1) with a smooth, polynomially bounded potential V, with initial data given by a linear combination of Hagedorn function with multi-indices $k\in {\mathcal K}$ :

$$ \begin{align*} \psi(\cdot,0) = \sum_{k\in{\mathcal K}} c_k(0) \varphi_{k_0}^\varepsilon[q(0),p(0),Q(0),P(0)], \end{align*} $$

where $Q(0),P(0)$ satisfy the symplecticity relations (3.10) and the coefficients are bounded by

$$ \begin{align*} |c_k(0)| \le C_0 \,\varepsilon^{\mathrm{dist}(k,{\mathcal K}_0)/6}, \quad k\in{\mathcal K}. \end{align*} $$

Let $\psi _{\mathcal K}(\cdot ,t)$ denote the Galerkin approximation (4.13) with the parameters $(q(t),p(t),Q(t),P(t))$ satisfying the classical equations (4.12) and with coefficients $c_k(t)$ for $k\in {\mathcal K}$ determined by (4.16), with initial data $\psi _{\mathcal K}(\cdot ,0)= \psi (\cdot ,0)$ . Then the error is bounded by

$$ \begin{align*} \| \psi_{\mathcal K}(\cdot,t) - \psi(\cdot,t) \| \le C \varepsilon^{N/2} \quad\text{for } 0\le t\le T, \end{align*} $$

where C is independent of $\varepsilon $ but depends on ${\mathcal K}_0$ , ${\mathcal K}$ , $C_0$ , on $\max _{0\le t \le T}\|Q(t)\|_2$ and T.

The above proof shows that as long as the position matrix satisfies a bound $\| Q(t) \|_2 \le M$ and the positions remain in a fixed compact set, the error is bounded by

$$ \begin{align*} \| \psi_{\mathcal K}(\cdot,t) - \psi(\cdot,t) \| \le C_M \,t\,\varepsilon^{N/2}, \end{align*} $$

with $C_M$ independent of t, so that potentially the approximation remains accurate over very long time scales. The growth of $\| Q(t) \|_2$ is determined by the linearized classical equations (4.12). In the case of a positive Lyapunov exponent, the time scale cannot exceed a constant times $\log \varepsilon ^{-1}$ , but in more favourable cases, such as integrable and near-integrable systems, where $\| Q(t) \|_2$ has only linear growth for all (or extremely long) times, the time scale of approximation can be longer; see the discussion of the Ehrenfest time in Section 3.6.

4.5 Notes

The construction and results in this section are due to Hagedorn (Reference Hagedorn1998b , Reference Hagedorn1981). Beyond the estimates of finite order in $\varepsilon $ shown here, Hagedorn and Joye (Reference Hagedorn and Joye1999b , Reference Hagedorn and Joye2000) derived exponentially small estimates in the case of an analytic potential.

Hagedorn (Reference Hagedorn1998b ) also showed that his functions behave well under the scaled Fourier transform

$$ \begin{align*} {{\mathcal F}_\varepsilon \varphi(\xi) = (2\pi\varepsilon)^{-d/2} \int_{{\mathbb R}^d} \varphi(x)\, \mathrm{e}^{-\mathrm{i} \xi\cdot x/\varepsilon}\, \mathrm{d} x ,\quad \xi \in \mathbb{R}^d.} \end{align*} $$

That is, for every k,

$$ \begin{align*} {\mathcal F}_\varepsilon \varphi_k^\varepsilon [q,p,Q,P](\xi) = \mathrm{e}^{-\mathrm{i} p\cdot q/\varepsilon} \varphi_k^\varepsilon [p,-q,P,-Q](\xi). \end{align*} $$

Lasser and Troppmann (Reference Lasser and Troppmann2014) gave transformation properties of the Hagedorn functions under Wigner and FBI (Fourier–Bros–Iagolnitzer) transforms and derived remarkable properties of the polynomial factor of the Hagedorn functions. The generating function for these polynomials was given by Hagedorn (Reference Hagedorn2015) and Dietert, Keller and Troppmann (Reference Dietert, Keller and Troppmann2017).

Hagedorn’s semiclassical wave packets are closely related to, or coincide with, generalized coherent states or generalized squeezed states, as studied by Combescure (Reference Combescure1992), Robert (Reference Robert2007) and Combescure and Robert (Reference Combescure and Robert2012). The precise relationship was expounded by Lasser and Troppmann (Reference Lasser and Troppmann2014) and more recently by Ohsawa (Reference Ohsawa2019).

Ohsawa (Reference Ohsawa2018) derived differential equations for the variational approximation by a single Hagedorn function $\varphi _k^\varepsilon $ of arbitrary index k, for an approximate Hamiltonian. A fully variational approximation by semiclassical wave packets (4.13) appears not to have been considered in the literature, not least because a separation of the motion of the position and momentum parameters $[q,p,Q,P]$ from that of the coefficients $c_k$ is not feasible for a fully variational approximation, as opposed to the semi-variational approximation considered here.

Hagedorn wave packets were proposed as a computational tool for semiclassical quantum dynamics by Faou, Gradinaru and Lubich (Reference Faou, Gradinaru and Lubich2009); see also Gradinaru and Hagedorn (Reference Gradinaru and Hagedorn2014).

5.1 Continuous superpositions of thawed and frozen Gaussians

We consider semiclassical approximation for the Schrödinger equation (1.1) with a smooth, subquadratic potential function $V\colon {\mathbb R}^d\to {\mathbb R}$ and general initial data $\psi _0$ . For convenience we take the initial wave function in the class of complex-valued Schwartz functions ${\mathcal S}({\mathbb R}^d)$ , but we note that our error bounds extend directly to general $L^2$ initial data by density. We recall that a Schwartz function on ${\mathbb R}^d$ is an infinitely differentiable function that together with its derivatives decays faster than the inverse of any polynomial.

Our main tool for the representation of wave functions will be the wave packet transform. We again let $\langle \cdot \mid \cdot \rangle $ denote the $L^2$ inner product on ${\mathbb R}^d$ and let $\|\cdot \|$ denote the $L^2$ -norm.

Proposition 5.1 (wave packet transform)

For any Schwartz function $g\colon {\mathbb R}^d\to {\mathbb C}$ of unit norm, we define the corresponding wave packetFootnote ¹

$$ \begin{align*} g_z(x) = \varepsilon^{-d/4}\, g \biggl(\dfrac{x-q}{\sqrt\varepsilon}\biggr)\, \mathrm{e}^{\mathrm{i} p\cdot (x-q)/\varepsilon},\quad x\in{\mathbb R}^d \end{align*} $$

for $z=(q,p)\in {\mathbb R}^{2d}$ . Then, for every Schwartz function $\psi \in {\mathcal S}({\mathbb R}^d)$ ,

$$ \begin{align*} \psi (x) &= (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi\rangle\, g_z (x) {\,\mathrm{d}} z,\quad x\in{\mathbb R}^d,\\ \|\psi\|^2 &= (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} |\langle g_z \mid \psi\rangle|^2 {\,\mathrm{d}} z. \end{align*} $$

Proof. We use the inversion formula for the Fourier transform: a Schwartz function f is reconstructed from its scaled Fourier transform

$$ \begin{align*} {\mathcal F}_\varepsilon f(\xi) = (2\pi\varepsilon)^{-d/2} \int_{{\mathbb R}^d}f(x) \,\mathrm{e}^{-\mathrm{i} \xi\cdot x/\varepsilon}{\,\mathrm{d}} x \end{align*} $$

by the inversion formula

$$ \begin{align*} f(x) = (2\pi\varepsilon)^{-d/2} \int_{{\mathbb R}^d} {\mathcal F}_\varepsilon f(\xi) \,\mathrm{e}^{\mathrm{i} \xi\cdot x/\varepsilon}{\,\mathrm{d}} \xi . \end{align*} $$

For a Schwartz function $\psi $ and for $x\in {\mathbb R}^d$ , this yields

$$ \begin{align*} &(2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi\rangle\, g_z(x) {\,\mathrm{d}} z \\ &\quad = (2\pi\varepsilon)^{-d} \varepsilon^{-d/2} \int_{{\mathbb R}^{3d}} \overline g\biggl(\displaystyle \frac{y-q}{\sqrt\varepsilon}\biggr) \, g\biggl(\displaystyle \frac{x-q}{\sqrt\varepsilon}\biggr)\, \mathrm{e}^{\mathrm{i} p\cdot(x-y)/\varepsilon}\,\psi(y){\,\mathrm{d}} (y,q,p)\\ &\quad = \varepsilon^{-d/2} \int_{{\mathbb R}^d} (2\pi\varepsilon)^{-d/2} \int_{{\mathbb R}^d} {\mathcal F}_\varepsilon\biggl( \overline g\biggl(\displaystyle \frac{\cdot-q}{\sqrt\varepsilon}\biggr) \, \psi \biggr)(p)\, \mathrm{e}^{\mathrm{i} p\cdot x/\varepsilon} {\,\mathrm{d}} p \ g\biggl(\displaystyle \frac{x-q}{\sqrt\varepsilon}\biggr) {\,\mathrm{d}} q \\ &\quad = \varepsilon^{-d/2} \int_{{\mathbb R}^{d}} \biggl|g\biggl(\displaystyle \frac{x-q}{\sqrt\varepsilon}\biggr)\biggr|^2\,\psi(x){\,\mathrm{d}} q \\ & \quad = \psi(x), \end{align*} $$

since the normalization of the function g implies

$$ \begin{align*} \varepsilon^{-d/2} \int_{{\mathbb R}^{d}} \biggl|g\biggl(\displaystyle \frac{x-q}{\sqrt\varepsilon}\biggr)\biggr|^2\,{\,\mathrm{d}} q = \int_{{\mathbb R}^d} |g(y)|^2 {\,\mathrm{d}} y = 1. \end{align*} $$

Moreover, by the inversion formula proved above,

$$ \begin{align*} \|\psi\|^2 &= \biggl\langle \psi\mid (2\pi\varepsilon)^{-d}\int_{{\mathbb R}^{2d}}\langle g_z \mid \psi\rangle g_z {\,\mathrm{d}} z\biggr\rangle \\ &= (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi\rangle \langle \psi \mid g_z\rangle {\,\mathrm{d}} z \\ &= (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} |\langle g_z \mid \psi\rangle|^2 {\,\mathrm{d}} z.\end{align*} $$

We note that if $\psi $ is a Schwartz function, then $z\mapsto \langle g_z\mid \psi \rangle $ is also a Schwartz function. This property ensures the existence of many integrals appearing below.

We use the inversion formula of Proposition 5.1 for the standard normalized Gaussian profile

$$ \begin{align*} g(x) = \pi^{-d/4} \exp\biggl(-\displaystyle \frac{1}{2} |x|^2\biggr),\quad x\in{\mathbb R}^d \end{align*} $$

and represent the initial data for the Schrödinger evolution as

$$ \begin{align*} \psi_0 = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\, g_z {\,\mathrm{d}} z. \end{align*} $$

We will now explore two related possibilities of transforming the identity

$$ \begin{align*} \mathrm{e}^{-\mathrm{i} t H/\varepsilon}\psi_0 = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\, \mathrm{e}^{-\mathrm{i} t H/\varepsilon} g_z {\,\mathrm{d}} z \end{align*} $$

into semiclassical approximations. Both approaches use the same basic ingredients from classical mechanics.

• • The flow $\Phi ^t\colon {\mathbb R}^{2d}\to {\mathbb R}^{2d}$ of the classical equations of motion

  (5.1) $$ \begin{align} \dot q = p,\quad \dot p = -\nabla V(q), \end{align} $$
  which associates the solution $\Phi ^t(z)=(q(t,\,z),p(t,\,z))$ with the initial datum $z\in {\mathbb R}^{2d}$ .

• • The linearized equations of motion, linearized at $q(t)$ ,

  (5.2) $$ \begin{align} \dot Q = P,\quad \dot P = -\nabla^2 V(q) Q{.} \end{align} $$

• • The action integral

  (5.3) $$ \begin{align} S(t,\,z) = \int_0^t \biggl(\displaystyle \frac{1}{2}|p(s,z)|^2 - V(q(s,z))\biggr) {\,\mathrm{d}} s \end{align} $$
  for the trajectory $\Phi ^t(z) = (q(t,\,z),p(t,\,z))$ with the initial datum $z\in {\mathbb R}^{2d}$ .

5.1.1 Thawed Gaussian approximation

We approximate an individual time-evolved Gaussian wave packet by the solution of the time-dependent Schrödinger equation with locally quadratic approximation to the potential:

$$ \begin{align*} \mathrm{i} \varepsilon\partial_t\varphi_z^{\mathrm{th}} = -\displaystyle \frac{\varepsilon^2}{2}\Delta\varphi_z^{\mathrm{th}}+ U_{q(t,\,z)}\varphi_z^{\mathrm{th}}, \quad \varphi_z^{\mathrm{th}}(0) = g_z, \end{align*} $$

where $U_q$ denotes the second-order Taylor polynomial of V expanded around the point q. This motivates us to consider the thawed Gaussian superposition

$$ \begin{align*} {{\mathcal I}}_{\mathrm{th}}(t)\psi_0 = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\, \varphi_z^{\mathrm{th}}(t) {\,\mathrm{d}} z \end{align*} $$

as an approximation to the Schrödinger solution $\psi (t)$ . From Proposition3.18 (in the version for a time-dependent quadratic potential) we know that $\varphi _z^{\mathrm {th}}(t)$ is a Gaussian wave packet of the form

$$ \begin{align*} \varphi_z^{\mathrm{th}}(t) = \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} \,(g[C(t,\,z)])_{\Phi^t(z)} \end{align*} $$

with a thawed Gaussian profile

$$ \begin{align*} g[C](x) = \pi^{-d/4} \det(\operatorname{Im } C)^{1/4} \exp\biggl(\displaystyle \frac{{\mathrm{i}}}{2} x^T C x\biggr), \quad x\in{\mathbb R}^d. \end{align*} $$

The width matrix $C(t,\,z)$ is determined by the matrix Riccati equation

(5.4) $$ \begin{align} \begin{aligned} \dot C(t,\,z) &= - C(t,\,z)^2 - \nabla^2 V(q(t,\,z)), \\ C(0,z) &= \mathrm{i} \,\mathrm{Id}_d, \end{aligned} \end{align} $$

or equivalently in Hagedorn’s parametrization (see Section 3.5) as

$$ \begin{align*} C(t,\,z)=P(t,\,z)Q(t,\,z)^{-1}, \end{align*} $$

where $Q(t,\,z)$ and $P(t,\,z)$ are the solution to the linearized classical equations of motion (5.2) corresponding to the linearization at $q(t,\,z)$ and with initial data

$$ \begin{align*} Q(0,z)=\mathrm{Id}, \quad P(0,z)=\mathrm{i} \, \mathrm{Id}. \end{align*} $$

We may thus rewrite the thawed approximation as

$$ \begin{align*} {{\mathcal I}}_{\mathrm{th}}(t)\psi_0 = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\, \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} \,(g[C(t,\,z)])_{\Phi^t(z)} {\,\mathrm{d}} z. \end{align*} $$

Let us compare the classically thawed wave packet $\varphi _z^{\mathrm {th}}(t)$ with its variational cousin $u(t)$ .

• • Both functions provide approximations to the Schrödinger solution $\psi (t)$ of order $\sqrt \varepsilon $ in the $L^2$ -norm and of order $\varepsilon $ for expectation values. (Analogous arguments to those for the proof of Theorem 3.5 apply to $\varphi _z^{\mathrm {th}}(t)$ as well.)

• • Both functions conserve norm, while only the variational Gaussian $u(t)$ is energy-conserving.

• • The classical equations of motion require point evaluations of the potential and its derivatives and are thus computationally less demanding than the variational equations of motion, which are built on averages.

We will find that the oscillations of the different Gaussians in the continuous superposition ${{\mathcal I}}_{\mathrm {th}}(t)\psi _0$ improve accuracy, and we will prove an error estimate of order $\varepsilon $ in the $L^2$ -norm (Theorem 5.2).

5.1.2 Frozen Gaussian approximation

The Herman–Kluk propagator is based on a different approach for the dynamics of the individual building blocks. They are still Gaussian wave packets with classically moving phase space centres, but their width matrices are kept frozen. The low accuracy of such a frozen Gaussian

$$ \begin{align*} \varphi_z^{\mathrm{fr}}(t) = \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} g_{\Phi^t(z)} \end{align*} $$

is compensated by an amazing reweighting in phase space: using the linearized classical motion we construct a smooth function $a_\natural \colon {\mathbb R}\times {\mathbb R}^{2d}\to {\mathbb C}$ , known as the Herman–Kluk prefactor, such that

$$ \begin{align*} {{\mathcal I}}_\natural(t)\psi_0 = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\, a_\natural(t,\,z)\,\mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon}\, g_{\Phi^t(z)} {\,\mathrm{d}} z \end{align*} $$

approximates the Schrödinger solution $\psi (t)$ with the same asymptotic accuracy as the superposition of thawed Gaussians, that is, with an error of order $\varepsilon $ with respect to the $L^2$ -norm (Theorem 5.3).

5.2 Accuracy of the Gaussian superposition and numerical algorithm

5.2.1 Thawed Gaussians

We consider an origin-centred Gaussian $g[C]_0$ with complex symmetric width matrix $C\in {\mathbb C}^{d\times d}$ with positive definite imaginary part. The semiclassical Fourier transform of such a function is a Gaussian with width matrix $-C^{-1}$ , that is,

$$ \begin{align*} {\mathcal F}_\varepsilon g[C]_0 = g[-C^{-1}]_0. \end{align*} $$

A short calculation analogous to the one in Lemma 3.16 shows that Hagedorn’s parametrization of $C=PQ^{-1}$ allows us to write the imaginary part of the matrix $-C^{-1}$ as

$$ \begin{align*} \operatorname{Im } (-C^{-1}) = (PP^*)^{-1}. \end{align*} $$

We thus have for the minimal eigenvalues that

$$ \begin{align*} \lambda_{\mathrm{min}}(\operatorname{Im } C) &= \|Q\|^{-2},\\ \lambda_{\mathrm{min}}(\operatorname{Im }(-C^{-1})) &= \|P\|^{-2}. \end{align*} $$

Both the width of a thawed Gaussian in position and in momentum space play a crucial role when analysing the accuracy of the continuous thawed superposition. The following spectral parameter will be important for estimating the approximation error.

Definition 5.1 (spectral parameter)

Let $C\in {\mathbb C}^{d\times d}$ be a complex symmetric matrix with positive definite imaginary part. We set $\rho _Q = \lambda _{\mathrm {min}}(\operatorname{Im } C)$ and $\rho _P = \lambda _{\mathrm {min}}(\operatorname{Im } (-C^{-1}))$ and call the positive number

$$ \begin{align*} \rho_* = \dfrac{\rho_Q\rho_P}{2\rho_Q+2\rho_P} = \displaystyle \frac{1}{2} (\|Q\|^2+\|P\|^2)^{-1} \end{align*} $$

the spectral parameter of the matrix C.

A uniform lower bound on the spectral parameter is enough to prove that the continuous thawed superposition is accurate of order $\varepsilon $ for arbitrary initial data that are Schwartz functions.

Theorem 5.2 (thawed Gaussian approximation)

We now consider the thawed Gaussian approximation

$$ \begin{align*} {{\mathcal I}}_{\mathrm{th}}(t)\psi_0 = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\, \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} \,(g[C(t,\,z)])_{\Phi^t(z)} {\,\mathrm{d}} z \end{align*} $$

for an arbitrary Schwartz function $\psi _0\colon {\mathbb R}^d\to {\mathbb C}$ . We assume that the potential function $V\colon {\mathbb R}^d\to {\mathbb R}$ is smooth and its derivatives of order $\ge 2$ are all bounded. We also assume that for the solutions $C(t,\,z)$ of the Riccati equation (5.4), there exists a positive lower bound $\rho>0$ such that the spectral parameters satisfy

$$ \begin{align*} \rho_*(t,\,z) \ge \rho \quad \text{for all }\ (t,\,z)\in[\,0,\bar t ]\times{\mathbb R}^{2d}. \end{align*} $$

Then the solution $\psi (t)$ of the time-dependent Schrödinger equation with potential V and initial data $\psi _0$ satisfies

$$ \begin{align*} \|\psi(t) - {{\mathcal I}}_{\mathrm{th}}(t)\psi_0\| \le c \,t\,\varepsilon\,\|\psi_0\|,\quad 0\le t\le \bar t, \end{align*} $$

where the constant $c<\infty $ depends on the spectral bound $\rho $ , on derivative bounds of V and on $\bar t$ but is independent of $\varepsilon $ , $\psi _0$ and $t\le \bar t$ .

The proof proceeds in three steps and is given later on. We will establish the following.

First step: basic norm bounds

We consider oscillatory integral operators of the form

$$ \begin{align*} ({{\mathcal I}}\psi)(x) = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi\rangle\, a(z)\, (x-\Phi_q(z))^m \, (g[C(z)])_{\Phi(z)} {\,\mathrm{d}} z, \end{align*} $$

that have the following building blocks:

• • $a\colon {\mathbb R}^{2d}\to {\mathbb C}$ is a bounded function,

• • $\Phi =(\Phi _q,\Phi _p)\colon {\mathbb R}^{2d}\to {\mathbb R}^{2d}$ is a volume-preserving map,

• • $\{C(z) \mid z\in {\mathbb R}^{2d}\}$ is a family of complex symmetric matrices with positive definite imaginary part having a positive lower bound for its spectral parameters.

For a multi-index $m=(m_1,\ldots ,m_d)\in {\mathbb N}_0^d$ and for $x=(x_1,\ldots ,x_d)\in {\mathbb R}^d$ we write $x^m=x_1^{m_1}\cdots x_d^{m_d}$ and we denote $|m|=m_1+\cdots +m_d$ .

We derive norm bounds that are of order $\varepsilon ^{|m|/2}$ and depend on the norm of $\psi $ , the supremum of $|a(z)|$ , and spectral bounds for the matrices $C(z)$ .

Second step: norm bounds accounting for the collective oscillation

We analyse the phase of the oscillatory integral operator ${{\mathcal I}}_{\mathrm {th}}(t)$ taking into account the interplay between the action integral $S(t,\,z)$ and the Hamiltonian flow $\Phi ^t(z) =(q(t,\,z),\,p(t,\,z))$ . By an elegant integration by parts, we obtain norm bounds for integral operators of the form

$$ \begin{align*} ({{\mathcal I}}(t)\psi)(x) = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi\rangle\, (x-q(t,\,z) )^m \, \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} (g[C(t,\,z)])_{\Phi^t(z)} {\,\mathrm{d}} z, \end{align*} $$

that depend on the parity of $|m|$ and are of order $\varepsilon ^{\lceil |m|/2\rceil }$ , where $\lceil |m|/2\rceil $ denotes the smallest integer $\ge |m|/2$ .

Third step: defect calculation

We calculate the defect and apply the stability Lemma 3.6. Together with the norm bounds we will have proved Theorem 5.2.

5.2.2 Frozen Gaussians

The frozen Gaussian approximation relies on the same information from the classical dynamics as the thawed one. It uses the flow, the linearized flow, and the action integral. However, the linearization now defines a reweighting factor that allows us to keep the individual Gaussians of frozen unit width. Let us define this reweighting factor.

Definition 5.2 (Herman–Kluk prefactor)

Let $Q_\natural (t,\,z)$ and $P_\natural (t,\,z)$ denote the solution to the linearized equations of motion (5.2) subject to the initial conditions

$$ \begin{align*} Q_\natural(0,z) = \mathrm{Id}\quad\text{and}\quad P_\natural(0,z) = -\mathrm{i} \,\mathrm{Id}. \end{align*} $$

Define the Herman–Kluk matrix as the complex $d\times d$ matrix

$$ \begin{align*} M_\natural(t,\,z) = Q_\natural(t,\,z) + \mathrm{i} P_\natural(t,\,z), \end{align*} $$

whose invertibility will be proved in Lemma 5.6. Then the smooth complex-valued function $a_\natural \colon {\mathbb R}\times {\mathbb R}^{2d}\to {\mathbb C}$ ,

$$ \begin{align*} a_\natural(t,\,z) = \sqrt{2^{-d}\det(M_\natural(t,\,z))}, \end{align*} $$

is called the Herman–Kluk prefactor. The branch of the square root is determined by continuity with respect to time.

We note that, with $D\Phi ^t(z)$ denoting the Jacobian matrix of the flow map,

$$ \begin{align*} M_\natural(t,\,z) &= (\mathrm{Id}, \mathrm{i} \, \mathrm{Id}) \, D\Phi^t(z) \begin{pmatrix} \mathrm{Id} \\ - \mathrm{i} \, \mathrm{Id} \end{pmatrix} \\ &=\partial_q q(t,\,z)-\mathrm{i} \partial_p q(t,\,z) +\partial_p p(t,\,z) + \mathrm{i} \partial_q p(t,\,z) . \end{align*} $$

In contrast to the thawed evolution, which uses the linearized flow for multiplying its components to build the thawed width matrices, the frozen dynamics sums the corresponding matrices and takes a determinant. Either way, we obtain continuous superpositions that are first-order accurate with respect to the semiclassical parameter $\varepsilon $ .

Theorem 5.3 (Herman–Kluk propagator)

We consider the Herman–Kluk propagator

$$ \begin{align*} ({{\mathcal I}}_\natural(t)\psi_0)(x) = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\,a_\natural(t,\,z)\, \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} \,g_{\Phi^t(z)}(x) {\,\mathrm{d}} z, \quad x\in {\mathbb R}^d, \end{align*} $$

for arbitrary initial data $\psi _0\colon {\mathbb R}^d\to {\mathbb C}$ that are Schwartz functions. We assume that the potential function $V\colon {\mathbb R}^d\to {\mathbb R}$ is smooth and its derivatives of order $\ge 2$ are all bounded. Then the solution $\psi (t)$ of the time-dependent Schrödinger equation with initial data $\psi _0$ satisfies

$$ \begin{align*} \|\psi(t) - {{\mathcal I}}_\natural(t)\psi_0\| \le c \,t\,\varepsilon\,\|\psi_0\|,\quad 0\le t\le \bar t, \end{align*} $$

where the constant $c <\infty $ depends on derivative bounds of V and on $\bar t$ but is independent of $\varepsilon $ , $\psi _0$ and $t\le \bar t$ .

The proof strategy for the frozen approximation is analogous to the thawed one. Let us briefly comment on its three main steps.

First step: basic norm bounds

We establish basic norm bounds for frozen integral operators that contain polynomial terms of the form $(x-\Phi _q(z))^m$ . Again these bounds are of order $\varepsilon ^{|m|/2}$ . However, their proof is considerably simpler than in the thawed case, since the frozen Gaussian profiles allow for a direct Fourier inversion argument that is unfortunately not applicable for the thawed operators.

Second step: norm bounds accounting for the collective oscillations

The second step of the proof consists again in analysing the phase of the oscillatory integral operator and provides the improved norm bound that is of order $\lceil |m|/2\rceil $ . The calculations are slimmer than in the thawed case, since there are no quadratic polynomials generated by derivatives of the width matrices.

Third step: defect calculation

The third step calculates the defect of the frozen approximation and performs an integration by parts for determining the Herman–Kluk prefactor. In contrast to the first two steps of the proof, here the analysis of the frozen approximation is a bit more demanding than that of the thawed one.

Remark 5.1 (mass and energy conservation)

Both the thawed and the frozen Gaussian approximation fail to be mass- or energy-conserving unless the potential is a polynomial of degree $\le 2$ .

5.3 Norm bounds for continuous superpositions of frozen Gaussians

Since the first step for proving the error estimate for the frozen Gaussian approximation is considerably less demanding than for the thawed one, we start by deriving basic norm bounds for oscillatory integral operators that generalize the Herman–Kluk propagator. This first basic estimate requires neither that the profile function is Gaussian nor that the time evolution stems from a classical Hamiltonian flow.

Proposition 5.4 (norm bound, general profile)

We assume the following.

1. (1) Let $g\in {\mathcal S}({\mathbb R}^d)$ be of unit norm.

2. (2) Let $a\colon {\mathbb R}^{2d}\to {\mathbb C}$ be a measurable and bounded function.

3. (3) Let $\Phi =(\Phi _q,\Phi _p)\colon {\mathbb R}^{2d}\to {\mathbb R}^{2d}$ be a volume-preserving diffeomorphism.

For $m\in {\mathbb N}_0^d$ we define, for $\psi \in {\mathcal S}({\mathbb R}^d)$ ,

$$ \begin{align*} ({{\mathcal I}} \psi)(x) = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi\rangle\, a(z)\,(x-\Phi_q(z))^m\, g_{\Phi(z)}(x) {\,\mathrm{d}} z. \end{align*} $$

Then ${{\mathcal I}}\psi $ is square-integrable and satisfies

$$ \begin{align*} \|{{\mathcal I}}\psi\| \le c_m \,\varepsilon^{|m|/2} \sup_{z\in{\mathbb R}^{2d}}|a(z)| \, \|\psi\|, \end{align*} $$

where $c_m>0$ depends on the mth moment of g. In particular, $c_0 = 1$ .

Proof. We use the inner products and the associated norms in the Hilbert spaces $L^2({\mathbb R}^d)$ and $L^2({\mathbb R}^{2d})$ . We distinguish them by subscripts x and z, respectively. Note that both inner products are conjugate linear in the first argument. Let $\varphi \in L^2({\mathbb R}^d)$ . Then

$$ \begin{align*} \langle \varphi \mid {{\mathcal I}}\psi\rangle_x = \langle (2\pi\varepsilon)^{-d/2} \langle (x-q)^m g_{z} \mid \varphi \rangle_x\circ\Phi \mid a(z)\, (2\pi\varepsilon)^{-d/2}\langle g_z \mid \psi\rangle_x \rangle_z. \end{align*} $$

Hence, by the Cauchy–Schwarz inequality,

$$ \begin{align*} |\langle \varphi \mid {{\mathcal I}}\psi\rangle_x| \le \|a\|_\infty\ (2\pi\varepsilon)^{-d/2} \|\langle (x-q)^m\,g_{z}\mid \varphi\rangle_x\circ\Phi \|_z (2\pi\varepsilon)^{-d/2}\| \langle g_z \mid \psi\rangle_x\|_z. \end{align*} $$

Since the profile function g has unit norm, the norm formula of Proposition 5.1 implies

$$ \begin{align*} (2\pi\varepsilon)^{-d}\|\langle g_z \mid \psi\rangle_x\|_z^2 = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} |\langle g_z \mid \psi\rangle_x|^2 {\,\mathrm{d}} z = \|\psi\|_x^2. \end{align*} $$

Since $\Phi $ is volume-preserving, we also have

$$ \begin{align*} \|\langle (x-q)^m \, g_{z} \mid \varphi\rangle_x\circ\Phi\|_z^2 = \|\langle (x-q)^m \, g_z \mid \varphi\rangle_x\|_z^2. \end{align*} $$

We observe that

$$ \begin{align*} {(x-q)^m g_z (x) = \varepsilon^{|m|/2} (x^m g)_z (x)} \end{align*} $$

and apply Proposition 5.1 to the normalized profile function $x^m g/\|x^m g\|$ . We obtain

$$ \begin{align*} (2\pi\varepsilon)^{-d} \|\langle (x-q)^m \, g_{z} \mid \varphi\rangle_x\circ\Phi\|_z^2 = \varepsilon^{|m|}\, \|x^m g\|_x^2\, \|\varphi\|_x^2. \end{align*} $$

In summary, setting

$$ \begin{align*} c_m = \|x^m g\|_x, \end{align*} $$

we have proved

$$ \begin{align*} |\langle \varphi \mid {{\mathcal I}}\psi\rangle_x| \le c_m\, \varepsilon^{|m|/2} \sup_{z\in{\mathbb R}^{2d}}|a(z)| \, \|\varphi\|_x\ \|\psi\|_x. \end{align*} $$

Proposition 5.4 ensures the well-definedness of the oscillatory Herman–Kluk integral. With

$$ \begin{align*} g(x) &= \pi^{-d/4} \exp\biggl(-\displaystyle \frac{1}{2}|x|^2\biggr),\quad m=0,\\ a(z) &= \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} a_\natural(t,\,z) \quad\text{and}\quad \Phi(z) = \Phi^t(z), \end{align*} $$

it implies that ${{\mathcal I}}_\natural (t)\psi _0$ is a square-integrable function with

$$ \begin{align*} \|{{\mathcal I}}_\natural(t)\psi_0\| \le \sup_{z\in{\mathbb R}^{2d}} |a_\natural(t,\,z)| \, \|\psi_0\|. \end{align*} $$

With a slight adjustment, the proof of Proposition 5.4 extends to a larger class of integral operators, that incorporates not just powers $(x-\Phi _q(z))^m$ but functions of the more general form

$$ \begin{align*} (x-\Phi_q(z))^m \, b(\Phi_q(z),x), \end{align*} $$

where $b\colon {\mathbb R}^{2d}\to {\mathbb R}$ is a measurable and bounded function. We will need this extension later on, when estimating the non-quadratic contributions of the potential to the approximation error.

Corollary 5.5 (norm bound, general profile)

Let $b\colon {\mathbb R}^{2d}\to {\mathbb C}$ be a measurable and bounded function. Under the assumptions of Proposition 5.4 define for any square-integrable $\psi \colon {\mathbb R}^d\to {\mathbb C}$

$$ \begin{align*} ({{\mathcal I}} \psi)(x) = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi\rangle\, a(z)\,(x-\Phi_q(z))^m\,b(\Phi_q(z),x)\, g_{\Phi(z)}(x) {\,\mathrm{d}} z. \end{align*} $$

Then ${{\mathcal I}}\psi $ is square-integrable and satisfies

$$ \begin{align*} \|{{\mathcal I}}\psi\| \le c_m \ \varepsilon^{|m|/2} \sup_{z\in{\mathbb R}^{2d}}|a(z)|\ \sup_{(x,z)\in{\mathbb R}^{3d}}|b(q,x)| \, \|\psi\|, \end{align*} $$

where $c_m>0$ depends on the moments of the profile function g. In particular, $c_0 = 1$ .

Proof. We literally repeat the proof of Proposition 5.4, adding the following observation. Denote $f(x) = x^m g(x)$ . By Fourier inversion we have

$$ \begin{align*} &\|\langle (x-q)^m \,b(q,x)\, g_z \mid \varphi\rangle_x\|_z^2\\ &\quad = \varepsilon^{|m|-d/2} \int_{{\mathbb R}^{4d}} \overline b(q,y) b(q,x) \overline f\biggl(\displaystyle \frac{y-q}{\sqrt\varepsilon}\biggr) f\biggl(\displaystyle \frac{x-q}{\sqrt\varepsilon}\biggr)\\ &\quad\quad\, \times \overline\varphi(y) \varphi(x) \,\mathrm{e}^{\mathrm{i} p\cdot(x-y)/\varepsilon}{\,\mathrm{d}} (x,y,z)\\ &\quad = (2\pi\varepsilon)^{d}\, \varepsilon^{|m|-d/2} \int_{{\mathbb R}^{2d}} |b(q,x)|^2 \, \biggl|f\biggl(\displaystyle \frac{x-q}{\sqrt\varepsilon}\biggr)\biggr|^2\, |\varphi(x)|^2 {\,\mathrm{d}} (x,q)\\ &\quad \le (2\pi\varepsilon)^{d}\, \varepsilon^{|m|}\, \sup_{(q,x)\in{\mathbb R}^{2d}}|b(q,x)|^2 \ \|f\|^2 \ \|\varphi\|^2. \end{align*} $$

So far we have not used the fact that the profile function of the Herman–Kluk propagator is Gaussian and we have not used the properties of the Hamiltonian flow, nor have we incorporated the phase contributions from the action integral. However, a more detailed analysis will allow us to recognize that the contributions from monomials $(x-q(t,\,z))^m$ can be smaller than expected, depending on the parity of $|m|$ . For this analysis, the second step for the proof of Theorem 5.3, we open the inner product integral in ${{\mathcal I}}_\natural (t)\psi _0$ and write

$$ \begin{align*} \langle g_z \mid \psi_0\rangle\,\mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} \,g_{\Phi^t(z)}(x) = (\pi\varepsilon)^{-d/2} \int_{{\mathbb R}^{d}} \,\mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon}\, \psi_0(y) {\,\mathrm{d}} y, \end{align*} $$

using a complex-valued phase function $\Psi (t,x,y,z)$ , that is quadratic with respect to $y-q$ and $x-q(t,\,z)$ . This function has remarkable properties.

Lemma 5.6 (phase function, frozen Gaussians)

The phase function

$$ \begin{align*} &\Psi\colon {\mathbb R}\times{\mathbb R}^d\times{\mathbb R}^d\times{\mathbb R}^{2d}\to{\mathbb C},\\ &\Psi(t,x,y,z) = \displaystyle \frac{{\mathrm{i}}}{2}(|y-q|^2+|x-q(t,\,z)|^2) \\ &\qquad\qquad\qquad- p\cdot(y-q) + p(t,\,z)\cdot(x-q(t,\,z)) + S(t,\,z) \end{align*} $$

satisfies for all $(t,x,y,z)$

$$ \begin{align*} (\mathrm{i} \partial_q + \partial_p)\Psi(t,x,y,z) = M(t,\,z)^T (x-q(t,\,z)), \end{align*} $$

where $M(t,\,z)$ denotes the complex $d\times d$ matrix

$$ \begin{align*} M(t,\,z) = \partial_q q(t,\,z) -\mathrm{i} \partial_p q(t,\,z) + \partial_p p(t,\,z) + \mathrm{i} \partial_q p(t,\,z). \end{align*} $$

The matrix $M(t,\,z)$ has the following properties.

1. (1) $M(t,\,z)$ is invertible with

   $$ \begin{align*} |\textrm{det}\, M(t,\,z)| = \sqrt{\det (\mathrm{Id}_{2d} + D\Phi^t(z)^T D\Phi^t(z))}. \end{align*} $$

2. (2) Its time derivative satisfies

   $$ \begin{align*} \partial_t M(t,\,z) = \partial_q p(t,\,z) -\mathrm{i} \partial_p p(t,\,z) - \nabla^2 V(q(t,\,z))(\partial_p q(t,\,z)+\mathrm{i} \partial_q q(t,\,z)) . \end{align*} $$

3. (3) If V is a polynomial of degree $\le 2$ , then $M(t)$ is independent of z.

Proof. We start by differentiating the action integral. By the classical equations of motion, we have

$$ \begin{align*} \partial_q S(t,\,z) &= \int_0^t (\partial_q p(s,z)^T p(s,z) - \partial_q q(s,z)^T\nabla V(q(s,z))) {\,\mathrm{d}} s\\ &= \int_0^t (\partial_q \dot q(s,z)^T p(s,z) + \partial_q q(s,z)^T\dot p(s,z)) {\,\mathrm{d}} s\\ &= \int_0^t \displaystyle \frac{\mathrm{d}}{\mathrm{d} s}(\partial_q q(s,z)^T p(s,z)) {\,\mathrm{d}} s \\ & = \partial_q q(t,\,z)^T p(t,\,z)-p \end{align*} $$

since $\partial _q q(0,z) = \mathrm {Id}_d$ . Analogously,

$$ \begin{align*} \partial_p S(t,\,z) = \int_0^t \displaystyle \frac{\mathrm{d}}{\mathrm{d} s}(\partial_p q(s,z)^T p(s,z)) {\,\mathrm{d}} s = \partial_p q(t,\,z)^T p(t,\,z). \end{align*} $$

Now we compute the derivatives of the phase function,

$$ \begin{align*} \partial_q\Psi(t,x,y,z) & = -{\mathrm{i}}(y-q) -\mathrm{i} \,\partial_q q(t,\,z)^T(x-q(t,\,z)) + p \\ &\quad\, + \partial_q p(t,\,z)^T(x-q(t,\,z))\\ &\quad\, - \partial_q q(t,\,z)^Tp(t,\,z) + (\partial_q q(t,\,z)^T p(t,\,z)-p)\\ & = -{\mathrm{i}}(y-q) +(\partial_q p(t,\,z)-\mathrm{i} \,\partial_q q(t,\,z))^T(x-q(t,\,z)) \end{align*} $$

and

$$ \begin{align*} \partial_p\Psi(t,x,y,z) &= -\mathrm{i} \,\partial_p q(t,\,z)^T(x-q(t,\,z)) -(y-q) \\ &\quad\, + \partial_p p(t,\,z)^T(x-q(t,\,z))\\ &\quad\, - \partial_p q(t,\,z)^Tp(t,\,z) + \partial_p q(t,\,z)^Tp(t,\,z)\\ &= (\partial_p p(t,\,z)-\mathrm{i} \,\partial_p q(t,\,z))^T(x-q(t,\,z)) -(y-q). \end{align*} $$

This implies

$$ \begin{align*} (\mathrm{i} \partial_q + \partial_p)\Psi(t,x,y,z) = M(t,\,z)^T(x-q(t,\,z)). \end{align*} $$

To prove invertibility, we work with the Jacobian matrix of the flow,

$$ \begin{align*} D\Phi^t(z) = \begin{pmatrix} \partial_q q(t,\,z) & \partial_p q(t,\,z)\\ \partial_q p(t,\,z) & \partial_p p(t,\,z)\end{pmatrix}\!{,} \end{align*} $$

and the matrix $F = F(t,\,z)$ defined by

$$ \begin{align*} F :&= \mathrm{Id}_{2d} + D\Phi + \mathrm{i} J(\mathrm{Id}_{2d}-D\Phi) \\ &= \begin{pmatrix}\mathrm{Id}_d + \partial_q(q+\mathrm{i} p) & \partial_p(q+\mathrm{i} p) -\mathrm{i} \,\mathrm{Id}_d\\ \partial_q(p-\mathrm{i} q) + \mathrm{i} \,\mathrm{Id}_d & \mathrm{Id}_d + \partial_p(p-\mathrm{i} q)\end{pmatrix}\!. \end{align*} $$

Since the two left blocks of F commute, its determinant satisfies

$$ \begin{align*} \det(F) &= \det(F_{11}F_{22}-F_{21}F_{12})\\ &=\det(\partial_q(q+\mathrm{i} p) + \partial_p(p-\mathrm{i} q) + \mathrm{i} \partial_q(p-\mathrm{i} q) -\mathrm{i} \partial_p(q+\mathrm{i} p))\\ & = \det(2M) \\ & = 2^d \det(M). \end{align*} $$

We calculate

$$ \begin{align*} F^*F &= (\mathrm{Id}_{2d} + D\Phi)^T(\mathrm{Id}_{2d} + D\Phi) + (\mathrm{Id}_{2d}-D\Phi)^T(\mathrm{Id}_{2d}-D\Phi)\\ &= 2\mathrm{Id}_{2d} + 2D\Phi^T D\Phi, \end{align*} $$

where we have used that the cross-terms involving J vanish,

$$ \begin{align*} &{\mathrm{i}}(\mathrm{Id}_{2d} + D\Phi)^T J (\mathrm{Id}_{2d}-D\Phi) - \mathrm{i} (\mathrm{Id}_{2d}-D\Phi)^TJ^T(\mathrm{Id}_{2d} + D\Phi)\\ &\quad = \mathrm{i} (-J D\Phi+ D\Phi^TJ) + {\mathrm{i}}(JD\Phi-D\Phi^T J) \\ &\quad = 0, \end{align*} $$

due to the symplecticity of $D\Phi ^t(z)$ . Hence

$$ \begin{align*} |\textrm{det}\, M|^2 = 4^{-d} \det(F^* F) = \det(\mathrm{Id}_{2d} + D\Phi^T D\Phi). \end{align*} $$

For the time derivative we simply calculate

$$ \begin{align*} \partial_t M(t,\,z) = (\partial_q -\mathrm{i} \partial_p) p(t,\,z) - \nabla^2 V(q(t,\,z)) (\partial_p +\mathrm{i} \partial_q) q(t,\,z). \end{align*} $$

If the potential V is a polynomial of degree $\le 2$ , then the Jacobian $D\Phi ^t$ and thus $M(t)$ are independent of z.□

The key observation of Lemma 5.6 is that the Wirtinger derivative

$$ \begin{align*} (\mathrm{i} \partial_q+\partial_p)\Psi(t,x,y,z) \end{align*} $$

is independent of the variable y and linear in $x-q(t,\,z)$ . This allows for an integration by parts that turns monomial terms $(x-q(t,\,z))^m$ into powers of the form $\varepsilon ^{\lceil |m|/2\rceil }$ , where $\lceil |m|/2\rceil $ denotes the smallest integer $\ge |m|/2$ . This explains the high accuracy of appropriate continuous superpositions of frozen Gaussians, which is of order $\varepsilon $ instead of $\sqrt \varepsilon $ .

Proposition 5.7 (norm bound, frozen Gaussians)

Let the functions $a,b \colon {\mathbb R}^{2d}\to {\mathbb C}$ be smooth and bounded. Let

$$ \begin{align*} g(x) = \pi^{-d/4}\exp\biggl(-\displaystyle \frac{1}{2}|x|^2\biggr), \quad x\in{\mathbb R}^d, \end{align*} $$

be the standard normalized Gaussian and let $m\in {\mathbb N}_0^d$ . For any square-integrable function $\psi \colon {\mathbb R}^d\to {\mathbb C}$ define

$$ \begin{align*} & {{\mathcal I}}(t)\psi \\ & \quad = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi\rangle \,a(t,\,z)\, (x-q(t,\,z))^m \,b(q(t,\,z),x) \, \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} g_{\Phi^t(z)} {\,\mathrm{d}} z {.} \end{align*} $$

Then ${{\mathcal I}}(t)\psi $ is square-integrable and satisfies

$$ \begin{align*} \|{{\mathcal I}}(t)\psi\| \le \gamma_{m}\, c_{|m|}(a,b,\Phi^t)\ \varepsilon^{\lceil |m|/2\rceil} \|\psi\| {,} \end{align*} $$

where the constant $\gamma _{m}>0$ depends on the polynomial degree m and

$$ \begin{align*} & c_{|m|}(a,b,\Phi^t) \\ & \quad = \sup_{|\alpha|\le |m|,(x,z)\in{\mathbb R}^{3d}} \|\partial^\alpha_z (a(z) b(q,x) M(t,\,z)^{-T})\| \, \sup_{|\alpha|\le |m|, z\in{\mathbb R}^{2d}} |\partial^\alpha_z\Phi^t(z)|. \end{align*} $$

Proof. We perform an inductive proof over $|m|$ . The case $|m| = 0$ is already covered by the norm bound of Corollary 5.5. We therefore start with $|m|=1$ , where $m=e_j$ for some $j=1,\ldots ,d$ . We use the derivative formula of Lemma5.6,

$$ \begin{align*} (x-q(t,\,z))\, \mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon} = \displaystyle \frac{\varepsilon}{{\mathrm{i}}} \,M(t,\,z)^{-T}(\mathrm{i} \partial_q + \partial_p) \,\mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon}. \end{align*} $$

Integration by parts yields

$$ \begin{align*} &{{\mathcal I}}(t)\psi(x) \\ &\quad = \dfrac{1}{2^{d} (\pi\varepsilon)^{3d/2}} \int_{{\mathbb R}^{3d}} a(t,\,z) (x-q(t,\,z))_j\, b(q(t,\,z),x)\, \mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon} \psi(y) {\,\mathrm{d}}(y,z)\\ &\quad = \dfrac{\varepsilon{\mathrm{i}}}{2^{d} (\pi\varepsilon)^{3d/2}} \sum_{k=1}^d \int_{{\mathbb R}^{3d}} (\mathrm{i} \partial_q+\partial_p)_k (a(t,\,z)\,M(t,\,z)^{-T}_{jk}\, b(q(t,\,z),x)) \\ &\quad\quad\, \times\mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon} \psi(y) {\,\mathrm{d}}(y,z). \end{align*} $$

Using Corollary 5.5, we obtain some constant $\gamma _{e_j}>0$ such that

$$ \begin{align*} \|{{\mathcal I}}(t)\psi\| \le \gamma_{e_j}\, c_1(a,b,\Phi^t)\,\varepsilon\, \|\psi\|. \end{align*} $$

For the inductive step, we consider

$$ \begin{align*} &{{\mathcal I}}(t)\psi(x) \\ &\quad =\dfrac{1}{2^{d} (\pi\varepsilon)^{3d/2}} \int_{{\mathbb R}^{3d}} \! a(t,\,z) (x-q(t,\,z))^{m+e_j}\, b(q(t,\,z),x)\,\\ &\quad\quad\, \times \mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon} \psi(y) {\,\mathrm{d}}(y,z)\\ &\quad =\dfrac{\varepsilon{\mathrm{i}}}{2^{d} (\pi\varepsilon)^{3d/2}} \sum_{k=1}^d \int_{{\mathbb R}^{3d}} a_{jk}(t,x,z) \,\mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon} \psi(y) {\,\mathrm{d}}(y,z) \end{align*} $$

with

$$ \begin{align*} a_{jk}(t,x,z) &= (\mathrm{i} \partial_q+\partial_p)_k\, ( a(t,\,z)\,M(t,\,z)^{-T}_{jk}\, b(q(t,\,z),x))(x-q(t,\,z))^m\\ &\quad\, + a(t,\,z)\,M(t,\,z)^{-T}_{jk}\, b(q(t,\,z),x)\ (\mathrm{i} \partial_q+\partial_p)_k (x-q(t,\,z))^m. \end{align*} $$

We therefore write

$$ \begin{align*} {{\mathcal I}}(t) = \varepsilon({{\mathcal I}}_0(t) + {{\mathcal I}}_1(t)), \end{align*} $$

where the integral operator ${{\mathcal I}}_0(t)$ contains monomials in $x-q(t,\,z)$ of order $|m|$ , while those in ${{\mathcal I}}_1(t)$ are of order $|m|-1$ . By the inductive hypothesis, there exist some constants $\widetilde \gamma _0,\widetilde \gamma _1>0$ such that

$$ \begin{align*} \|{{\mathcal I}}_0(t)\psi\| \le \widetilde \gamma_0\, c_{|m|+1}(a,b,\Phi^t) \,\varepsilon^{\lceil |m|/2\rceil} \|\psi\| \end{align*} $$

and

$$ \begin{align*} \|{{\mathcal I}}_1(t)\psi\| \le \widetilde \gamma_1\, c_{|m|}(a,b,\Phi^t) \,\varepsilon^{\lceil (|m|-1)/2\rceil} \|\psi\|. \end{align*} $$

We have

$$ \begin{align*} \lceil(|m|-1)/2\rceil+1 = \begin{cases} |m|/2+1 = \lceil (|m|+1)/2\rceil & \text{for }m\text{ even,}\\ (|m|+1)/2 = \lceil (|m|+1)/2\rceil & \text{for }m\text{ odd,} \end{cases} \end{align*} $$

and therefore

$$ \begin{align*} \|{{\mathcal I}}(t)\psi\| \le \gamma_{m+e_j}\,c_{|m|}(a,b,\Phi^t) \,\varepsilon^{\lceil (|m|+1)/2\rceil} \, \|\psi\| \end{align*} $$

for some constant $\gamma _{m+e_j}>0$ .□

5.4 The Herman–Kluk propagator

We next calculate the defect of a frozen Gaussian superposition that is built on an arbitrary, not yet determined weight factor $a(t,\,z)$ .

Lemma 5.8 (defect calculation)

For a smooth and bounded function $a \colon {\mathbb R}\times {\mathbb R}^{2d}\to {\mathbb C}$ and square-integrable $\psi _0 \colon {\mathbb R}^d\to {\mathbb C}$ , we consider

$$ \begin{align*} {{\mathcal I}}(t)\psi_0 = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\,a(t,\,z)\, \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} \,g_{\Phi^t(z)}(x) {\,\mathrm{d}} z. \end{align*} $$

Then the defect

$$ \begin{align*} d(t) = \biggl(\displaystyle \frac{1}{\mathrm{i} \varepsilon} H-\partial_t\biggr){{\mathcal I}}(t)\psi_0 \end{align*} $$

satisfies

$$ \begin{align*} d(t) = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\,\delta(t)\, \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} \,g_{\Phi^t(z)} {\,\mathrm{d}} z, \end{align*} $$

where

$$ \begin{align*} &\delta(t,x,z) = -\partial_t a(t,\,z) + \displaystyle \frac{1}{\mathrm{i} \varepsilon} a(t,\,z)\biggl(-\displaystyle \frac{1}{2}|x-q(t,\,z)|^2 +\displaystyle \frac{\varepsilon d}{2}\biggr) \\ &\quad\, + \displaystyle \frac{1}{\mathrm{i} \varepsilon} a(t,\,z)\biggl( \displaystyle \frac{1}{2} (x-q(t,\,z))^T\, \nabla^2 V(q(t,\,z)) (x-q(t,\,z)) + W_{q(t,\,z)}(x)\biggr), \end{align*} $$

and

$$ \begin{align*} W_q = V - U_q \end{align*} $$

denotes the non-quadratic remainder of the potential V expanded around a point $q\in {\mathbb R}^d$ .

Proof. We start by calculating the two time derivatives, namely

$$ \begin{align*} \partial_t \,\mathrm{e}^{\mathrm{i} S(t)/\varepsilon} = -\displaystyle \frac{1}{\mathrm{i} \varepsilon}\biggl( \displaystyle \frac{1}{2}|p(t)|^2-V(q(t))\biggr) \,\mathrm{e}^{\mathrm{i} S(t)/\varepsilon} \end{align*} $$

and

$$ \begin{align*} \partial_t g_{\Phi^t} &= \partial_t\biggl(-\displaystyle \frac{1}{2\varepsilon}|x-q(t)|^2 + \displaystyle \frac{\mathrm{i}}{ \varepsilon} p(t)^T (x-q(t)) \biggr) g_{\Phi^t}\\ &= \displaystyle \frac{1}{\mathrm{i} \varepsilon}( (\mathrm{i} p(t)-\dot p(t))^T(x-q(t)) +|p(t)|^2) g_{\Phi^t}. \end{align*} $$

This implies that

$$ \begin{align*} &\partial_t( \mathrm{e}^{\mathrm{i} S(t)/\varepsilon}g_{\Phi^t}) \\ &\quad = \displaystyle \frac{1}{\mathrm{i} \varepsilon}\biggl( \displaystyle \frac{1}{2}|p(t)|^2 + V(q(t)) + (\mathrm{i} p(t)-\dot p(t))^T(x-q(t)) \biggr) \,\mathrm{e}^{\mathrm{i} S(t)/\varepsilon}g_{\Phi^t} \end{align*} $$

and consequently

$$ \begin{align*} &\partial_t \,{{\mathcal I}}(t)\psi_0 \\ &\quad = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\, \displaystyle \frac{1}{\mathrm{i} \varepsilon}\biggl( \mathrm{i} \varepsilon\partial_t a(t) + \biggl( \displaystyle \frac{1}{2}|p(t)|^2 + V(q(t)) \biggr) a(t) \biggr)\\ &\quad\quad\, \times \,\mathrm{e}^{{\mathrm{i} S(t)/\varepsilon}} g_{\Phi^t} {\,\mathrm{d}} z\\ &\quad\quad\, + (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\, \displaystyle \frac{1}{\mathrm{i} \varepsilon}\, a(t) \, (\mathrm{i} p(t)-\dot p(t))^T(x-q(t)) \,\mathrm{e}^{\mathrm{i} S(t)/\varepsilon}g_{\Phi^t} {\,\mathrm{d}} z. \end{align*} $$

We also calculate the space derivatives,

$$ \begin{align*} \partial_j g_z(x) = \biggl( -\displaystyle \frac1\varepsilon(x-q)_j + \displaystyle \frac{\mathrm{i}}{\varepsilon} p_j\biggr) g_z(x),\quad j=1,\ldots,d, \end{align*} $$

and

$$ \begin{align*} -\displaystyle \frac{\varepsilon^2}{2} \Delta_x g_z(x) = \biggl( -\displaystyle \frac{1}{2}|x-q|^2 + \mathrm{i} p^T(x-q) + \displaystyle \frac{1}{2}|p|^2 +\displaystyle \frac{\varepsilon d}{2}\biggr)g_z(x). \end{align*} $$

Adding the time and space derivatives together with a Taylor expansion of the potential, we obtain the claimed form of the defect. □

In view of the frozen norm bounds determined in Proposition 5.7, the calculated form of the defect is encouraging. The leading-order terms are

• • the time derivative of the reweighting function $a(t,\,z)$ ,

• • the contributions that are quadratic in $x-q(t,\,z)$ , one stemming from the Laplacian, the other from the quadratic expansion of the potential.

The non-quadratic remainder of the potential will contribute to the overall error, and we have a first indication of an integral representation of the dynamics that is exact for quadratic potentials.

Next, we explicitly calculate the terms generated by the integration by parts used for the improved norm bound in Proposition 5.7. We will obtain an explicit formula for the reweighting function, which solely depends on the classical flow map and provides an error of order $\varepsilon $ .

Proof of Theorem 5.3

We let $A(t,\,z)$ denote the real symmetric $d\times d$ matrix

$$ \begin{align*} A(t,\,z) = -\mathrm{Id}+\nabla^2 V(q(t,\,z)) \end{align*} $$

and write the defect determined in Lemma 5.8 as

$$ \begin{align*} \delta(t,x,z) = -\partial_t a(t,\,z) + \delta_2(t,x,z) + \displaystyle \frac{1}{\mathrm{i} \varepsilon}a(t,\,z) \biggl( \displaystyle \frac{\varepsilon d}{2} + W_{q(t,\,z)}(x)\biggr), \end{align*} $$

where the quadratic part is given by

$$ \begin{align*} \delta_2(t,x,z) = \displaystyle \frac{1}{2\mathrm{i} \varepsilon}\, a(t,\,z)\,(x-q(t,\,z))^T\,A(t,\,z) (x-q(t,\,z)) \end{align*} $$

and the non-quadratic remainder by

$$ \begin{align*} W_q(x) = \displaystyle \frac{1}{2} \sum_{|m|=3} (x-q)^m \int_0^1 (1-\theta)^2\ \partial^m V((1-\theta)q+\theta x) {\,\mathrm{d}}\theta. \end{align*} $$

We perform an integration by parts for the quadratic contributions. By the derivative formula of Lemma 5.6,

$$ \begin{align*} (x-q(t,\,z))\, \mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon} = \displaystyle \frac{\varepsilon}{{\mathrm{i}}} \,M(t,\,z)^{-T}(\mathrm{i} \partial_q + \partial_p) \,\mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon}, \end{align*} $$

we have

$$ \begin{align*} &\int_{{\mathbb R}^{2d}} \delta_2(t,x,z) \,\mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon} {\,\mathrm{d}} z \\[6pt] &\quad = -\int_{{\mathbb R}^{2d}} a(t,\,z)\,\displaystyle \frac{1}{2}(x-q(t,\,z))^T\,A(t,\,z) M(t,\,z)^{-T} (\mathrm{i} \partial_q + \partial_p)\,\mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon}{\,\mathrm{d}} z\\[6pt] &\quad = \int_{{\mathbb R}^{2d}} \widetilde \delta_2(t,x,z)\,\mathrm{e}^{\mathrm{i} \Psi(t,x,y,z)/\varepsilon} {\,\mathrm{d}} z \end{align*} $$

with

$$ \begin{align*} \widetilde \delta_2(t,x,z) = \displaystyle \frac{1}{2} \sum_{k,\ell=1}^d (\mathrm{i} \partial_q + \partial_p)_\ell ( a(t,\,z)\,(x-q(t,\,z))_k \,(A(t,\,z) M(t,\,z)^{-T})_{k,\ell}). \end{align*} $$

We calculate the derivative

$$ \begin{align*} (\mathrm{i} \partial_q + \partial_p)_\ell (x-q(t,\,z))_k = -(\mathrm{i} \partial_q q(t,\,z)^T + \partial_p q(t,\,z)^T)_{k,\ell}, \end{align*} $$

and obtain

$$ \begin{align*} &\displaystyle \frac{1}{2}\sum_{k,\ell=1}^d (A(t,\,z) M(t,\,z)^{-T})_{k,\ell} (\mathrm{i} \partial_q + \partial_p)_\ell (x-q(t,\,z))_k \\ &\quad = -\displaystyle \frac{1}{2}\ \mathrm{tr}((\mathrm{i} \partial_q q(t,\,z)^T + \partial_p q(t,\,z)^T) A(t,\,z)\, M(t,\,z)^{-T}). \end{align*} $$

By the phase function Lemma 5.6,

$$ \begin{align*} &- (\mathrm{i} \partial_q q^T + \partial_p q^T ) A - \mathrm{i} M^T \\ &\quad = (\mathrm{i} \partial_q q^T + \partial_p q^T ) (\mathrm{Id}-\nabla^2 V ) - ( \mathrm{i} \partial_qq^T + \mathrm{i} \partial_p p^T - \partial_q p^T + \partial_p q^T )\\ &\quad = \partial_t M^T. \end{align*} $$

We thus obtain

$$ \begin{align*} &{\widetilde\delta_2(t,x,z) + \displaystyle \frac{{d}}{2{\mathrm{i}}} = \displaystyle \frac{1}{2}\,\mathrm{tr} (\partial_t M(t,\,z)^T\ M(t,\,z)^{-T} ) \,a(t,\,z)} \\ &\quad\, + \displaystyle \frac{1}{2} \sum_{k,\ell=1}^d (x-q(t,\,z))_k (\mathrm{i} \partial_q + \partial_p )_\ell ( a(t,\,z)\,(A(t,\,z) M(t,\,z)^{-T})_{k,\ell} ). \end{align*} $$

This suggests choosing $a(t,\,z)$ as the solution of the ordinary differential equation

$$ \begin{align*} \partial_t a(t,\,z) = \displaystyle \frac{1}{2}\,\mathrm{tr} (\partial_t M(t,\,z)^T\ M(t,\,z)^{-T} ) a(t,\,z),\quad a(0,z) = 1. \end{align*} $$

By Liouville’s formula,

$$ \begin{align*} \partial_t \det(M(t,\,z)) = \det(M(t,\,z))\ \mathrm{tr}(\partial_t M(t,\,z) M(t,\,z)^{-1}). \end{align*} $$

Therefore

$$ \begin{align*} a(t,\,z) = \sqrt{2^{-d}\det(M(t,\,z))} = a_\natural(t,\,z). \end{align*} $$

The corresponding defect

$$ \begin{align*} d_\natural(t) = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\,\delta_\natural(t)\, \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} \,g_{\Phi^t(z)}{\,\mathrm{d}} z \end{align*} $$

is determined by the function

$$ \begin{align*} \delta_\natural(t,x,z) &= \displaystyle \frac{1}{2} \sum_{k,\ell=1}^d (x-q(t,\,z))_k (\mathrm{i} \partial_q + \partial_p )_\ell ( a_\natural(t,\,z)\,(A(t,\,z) M(t,\,z)^{-T})_{k,\ell} ) \\ &\quad\, + \displaystyle \frac{1}{\mathrm{i} \varepsilon}a_\natural(t,\,z) W_{q(t,\,z)}(x). \end{align*} $$

The function $\delta _\natural $ consists of two parts. The first one is linear in $x-q(t,\,z)$ . The frozen norm bound of Proposition 5.7 results in an upper bound for the corresponding oscillatory integral of order $\varepsilon ^{\lceil 1/2\rceil } = \varepsilon $ . The second part involving the non-quadratic remainder of the potential is cubic in $x-q(t,\,z)$ but divided by $\varepsilon $ , resulting in an upper bound that is of order

$$ \begin{align*} \varepsilon^{\lceil 3/2\rceil-1} = \varepsilon^{2-1} = \varepsilon \end{align*} $$

as well. Altogether we obtain

$$ \begin{align*} \|d_\natural(t)\| \le C\, \varepsilon\, \|\psi_0\| {,} \quad 0\le t\le \bar t, \end{align*} $$

where the constant $C>0$ is independent of $\varepsilon $ and $\psi _0$ . By the stability Lemma3.6, the error then satisfies

$$ \begin{align*} \|\psi(t) - {{\mathcal I}}_\natural(t)\psi_0\| \le \int_0^t \|d_\natural(s)\| {\,\mathrm{d}} s \le C\, t\, \varepsilon\,\|\psi_0\|. \end{align*} $$

As for the variational Gaussian approximation and the Hagedorn wave packets, we observe exactness for quadratic potentials.

Corollary 5.9 (exactness for quadratic potentials)

If the potential function V is quadratic, then the Herman–Kluk propagator is exact, that is, $\psi (t) = {{\mathcal I}}_\natural (t)\psi _0$ for all square-integrable initial data $\psi _0$ and all times t.

Proof. Quadratic potentials have a constant Hessian $\nabla ^2 V$ . Therefore the Jacobian of the flow, the matrices $A(t)$ and $M(t)$ , and the Herman–Kluk prefactor $a_\natural (t)$ do not depend on z. Consequently, $d_\natural (t) = 0$ . □

5.5 First norm bounds for continuous thawed Gaussian superpositions

The basic norm bounds we have used to analyse the frozen Gaussian superposition, that is, Proposition 5.4 and Corollary 5.5, do not directly apply to the thawed Gaussian superposition

$$ \begin{align*} {{\mathcal I}}_{\mathrm{th}}(t)\psi_0 = (2\pi\varepsilon)^{-d} \int_{{\mathbb R}^{2d}} \langle g_z \mid \psi_0\rangle\, \mathrm{e}^{\mathrm{i} S(t,\,z)/\varepsilon} g(t)_{\Phi^t(z)} {\,\mathrm{d}} z, \end{align*} $$

due to the z-dependence of the Gaussian width matrix $C(t,\,z)$ . It is natural to assume that the imaginary parts of the time-evolved family $C(t,\,z)$ are bounded from below in the sense that

$$ \begin{align*} \operatorname{Im } C(t,\,z) \ge \operatorname{Im } C_0 \quad {\text{for all}}\ (t,\,z)\in{\mathbb R}\times {\mathbb R}^{2d}, \end{align*} $$

where $C_0$ is a complex symmetric $d\times d$ matrix with $\operatorname{Im } C_0>0$ . Then it is tempting to write