How do wave packets spread? Time evolution on Ehrenfest time scales

The semiclassical propagation of wave packets is of fundamental importance in many applications of quantum mechanics and has therefore been studied extensively in the literature [1–9]. Time-dependent WKB approximations apply to extended initial states of the form

$$\begin{equation} \psi (x)=A(x)\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S(x)}, \end{equation} \tag{ 1 }$$

where the phase function S(x) is real valued and smooth, and the amplitude A(x) is smooth and non-oscillating. The time evolution of such states can be expressed in terms of transport along a family of classical trajectories determined by initial positions x and initial momenta p = ∇S(x); see, e.g., [2, 5].

A different class of initial states is given by coherent states, see, e.g., [3, 7], which are strongly localized, and therefore, their propagation can be described up to certain times using only one trajectory and the linearized motion around it. A Gaussian coherent state is a function of the form

$$\begin{equation} \psi _Z^B(x)=\frac{(\det\, {\rm Im}\,B)^{1/4}}{(\pi \hbar )^{n/4}}\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }[p\cdot (x-q)+(x-q)\cdot B (x-q)/2]}, \end{equation} \tag{ 2 }$$

where $Z=(p,q)\in \mathbb{R}^{2n}$ is a set of parameters, which determine where the state is localized in phase space and B is a complex symmetric n × n matrix with Im B > 0. The latter condition ensures that the state is localized around x = q. The roles of B and Z become clearer if we look at the Wigner function of such a state, which is a Gaussian of the form

$$\begin{equation} W(z)=\frac{1}{(\pi \hbar )^n}\,\mathrm{e}^{-(z-Z)\cdot G (z-Z)/\hbar }, \end{equation} \tag{ 3 }$$

with phase space coordinates $z\in \mathbb{R}^n\times \mathbb{R}^n$ and where G is a positive-definite symmetric and symplectic matrix determined by B. As was described by Hepp and Heller [10–13], the time evolution of a coherent state can in the semiclassical limit be described by using just the classical trajectory Z(t) = (p(t), q(t)) through Z and the linearized flow around it which is a time-dependent 2n × 2n symplectic matrix ${\mathcal S}(t)$. This simple description of the semiclassical propagation of coherent states made them a very useful tool in many applications, in particular in chemistry due to the work of Heller and coworkers and their appearance in initial value representations like the Herman–Kluk propagator [3, 4, 13–18]

We will be using time-dependent WKB methods developed for states of the form (1) to understand the propagation of the states (2) for long times. In order to identify the relevant time scales, it is useful to consider the Wigner function of the time-evolved coherent state, which in leading order is given by

$$\begin{equation} W(t,z)\approx \frac{1}{(\pi \hbar )^n}\,\mathrm{e}^{-(z-Z(t))\cdot G(t)(z-Z(t))/\hbar }, \end{equation} \tag{ 4 }$$

where $G(t)= {\mathcal S}^{\dagger }(t)G{\mathcal S}(t)$ and ${\mathcal S}(t)$ is the linearized flow around the central trajectory Z(t). The basic idea leading to this approximation is that since a coherent state is localized around a point in phase space, one can approximate the Hamiltonian by its quadratic Taylor expansion around the classical trajectory Z(t). The resulting Schrödinger equation can then be solved explicitly. This approximation can only be expected to be accurate as long as the propagated state stays localized, and from (4), we see that this is the case as long as λ_min[G(t)]/ℏ ≫ 1, where we denote by λ_min/max[G(t)] the smallest or the largest eigenvalue of the matrix $G(t)={\mathcal S}^{t}(t)G{\mathcal S}(t)$, respectively. Since G(t) is symmetric and positive, we have λ_min[G(t)] > 0, and that G(t) is symplectic implies that 1/λ_min[G(t)] = λ_max[G(t)], so the localization condition becomes

$$\begin{equation} \lambda _{\rm max}[G(t)]\hbar \ll 1 . \end{equation} \tag{ 5 }$$

In order to relate this more explicitly to the properties of the classical dynamics, we use the estimate $\lambda _{\rm max}[G(t)]=\Vert G(t)\Vert \le C \Vert {\mathcal S}(t)\Vert ^2$ and so the condition becomes $\Vert {\mathcal S}(t)\Vert \sqrt{\hbar }\ll 1$. Here, we use for a matrix A the norm induced by the Euclidean norm $\Vert x\Vert =\sqrt{ x\cdot x}$, i.e. $\Vert A\Vert :=\sup _{\Vert x\Vert =1} \frac{\Vert Ax\Vert }{\Vert x\Vert }$. The time at which the width of a propagated coherent state reaches order 1 is called the Ehrenfest time T_E(ℏ), and by the previous discussion, we see that it is determined by the condition

$$\begin{equation} \Vert {\mathcal S}(T_E(\hbar ))\Vert \sqrt{\hbar }= 1 . \end{equation} \tag{ 6 }$$

Hence, the Ehrenfest time depends on the dynamical properties of the classical system near the trajectory Z(t). If the trajectory Z(t) is hyperbolic, with the largest Liapunov exponent λ > 0, then ||S(t)|| ∼ e^λt and so

$$\begin{equation} T_E(\hbar )=\frac{1}{2\lambda }\ln {\frac{1}{\hbar }} . \end{equation} \tag{ 7 }$$

On the other hand, if the dynamics is integrable, we typically have $\Vert {\mathcal S}(t)\Vert \sim t$ for large t and then

$$\begin{equation} T_E(\hbar )=\frac{1}{\sqrt{\hbar }} . \end{equation} \tag{ 8 }$$

Ehrenfest showed in his classical paper [19] how one can recover classical mechanics from quantum dynamics by considering expectation values in propagated localized states. His construction breaks down if the state becomes delocalized, and is hence named Ehrenfest time for the time scale at which this occurs. In more recent years, the behaviour of time-evolved coherent states for long times has attracted the attention of mathematicians and rigorous estimates for long times have been derived; in particular, in [20, 21], it was shown that the remainder terms remain small for times satisfying $\vert t\vert < \frac{1}{3}T_E$ in the hyperbolic case. One should emphasize however that the breakdown of Ehrenfest's argument at the Ehrenfest time does not mean that the quantum to classical correspondence breaks down, it only means that the approximate propagation of a wave packet based on one central trajectory breaks down. If one includes more trajectories, one can use semiclassical propagation well beyond the Ehrenfest time as has been demonstrated in [22, 23], although from a mathematical point of view this is still a challenging problem.

The principal aim of this work is to present a propagation scheme for coherent states, which works at times of the order of the Ehrenfest time t ∼ T_E and is able to describe the transition from a localized state (2) to an extended state (1) at the Ehrenfest time in a uniform way. In a previous paper [25], it was demonstrated that at the Ehrenfest time a coherent state that is transported along a hyperbolic trajectory becomes effectively a WKB state of the type (1) associated with the unstable manifold of that trajectory; see also [24]. Here, we will give the theoretical foundation for that claim by presenting a new propagation scheme, which is a combination of a time-dependent WKB approximation and a metaplectic correction.

The qualitative change in the nature of the state at the Ehrenfest time was described in [26] for the cat map and in [27] for the propagation of coherent states centred on unstable fixed points in one-dimensional multiple well potentials. In contrast to the present approach, that work was not based on direct propagation of the state in the position representation but on the Wigner function and the Egorov theorem from [28]. That means in particular that we can give a much more detailed description of the state at and beyond the Ehrenfest time.

The plan of the paper is as follows. In section 2, we recall time-dependent WKB theory for the propagation of states of the form (1) and rewrite it in an exact form where we explicitly separate the classical transport and the quantum dispersion. The main idea of this paper is then developed in section 3 where we show how to approximate the action of the dispersive part on a coherent state using a simple metaplectic operator and combine this with the time-dependent WKB approximation to obtain a method that allows us to describe the propagation of coherent states at and beyond the Ehrenfest time. In section 4, we develop the phase space geometry underlying the metaplectically extended WKB method. This allows us to determine its range of validity, in particular its dependence on dynamical properties of the classical system. We then develop the geometric picture even further in section 5 to allow more general Hamiltonians and the inclusion of caustics after the Ehrenfest time; the price we have to pay is that remainder estimates become less explicit now. The last two sections are devoted to examples. In section 6, we consider a couple of simple explicitly or almost explicitly solvable examples, which demonstrate in some detail how our method works in practice. We show in particular that for a potential barrier we can describe the transition from transmission to reflection near the critical energy in a uniform way. In section 7, we reconsider the situation from [25] and demonstrate how our method reproduce the results. As a particular application, this shows using [29] that in a strongly chaotic system coherent states become semiclassically equidistributed beyond the Ehrenfest time. Finally, we summarize our results in the conclusions. In the appendix, we collect some more technical results from the semiclassical analysis that we use in the main body of the work.

Our method is an extension of the well-known time-dependent WKB method for real-valued phase functions. So we start by recalling this method, and to keep the discussion simple, we restrict ourselves first to the case of the standard Schrödinger equation with the potential V in $\mathbb{R}^n$, which in appropriate units reads

$$\begin{equation} \mathrm{i}\hbar \partial _t \psi =-\frac{\hbar ^2}{2}\Delta \psi +V\psi . \end{equation} \tag{ 9 }$$

In section 5, we will discuss a more general setting. We want to solve this equation for an initial state of the form $\psi =A_0\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S_0}$, where S₀ is real valued. Inserting the usual WKB ansatz

$$\begin{equation} \psi (t,x)=A(t,x)\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S(t,x)} , \end{equation} \tag{ 10 }$$

with S(t, x) being real valued, into the Schrödinger equation gives an expression that we separate, following the standard treatment [5], into the following two equations:

$$\begin{eqnarray} \partial _tS+\frac{1}{2}\vert \nabla S\vert ^2+V&=0 \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \mathrm{i}\partial _tA+\mathrm{i}\nabla S \cdot \nabla A +\frac{1}{2}\mathrm{i}(\Delta S) A&= \frac{\hbar }{2}\Delta A . \end{eqnarray} \tag{ 12 }$$

If A does not depend on ℏ, then this is a splitting into different powers of ℏ. Equation (11) is the Hamilton–Jacobi equation and its solutions are described using the propagation of the Lagrangian manifolds defined by p = ∇S(t, q):

$$\begin{equation} \Lambda _t:=\lbrace (\nabla S(t,x),x)\rbrace \subset \mathbb{R}^n\times \mathbb{R}^n , \end{equation} \tag{ 13 }$$

which are transported by the Hamiltonian flow associated with the classical Hamiltonian $H(p,q)=\frac{1}{2}p^2+V(q)$. We will analyse this further in section 4. The phase function S(t, x) often exists only on a finite time interval and develops singularities if caustics appear, then one has to use a refined ansatz. Because of the association with Lagrangian submanifolds WKB states of the form (10) are often called Lagrangian states in the mathematical literature [30, 31].

If we ignore in (12) the term of order ℏ on the right-hand side, then this is a transport equation that describes how the amplitude A is transported along the Lagrangian manifold Λ_t. Let us define the transport operator T(t) as the solution to the transport equation

$$\begin{equation} \mathrm{i}\partial _t T(t)=-\big[\mathrm{i}\nabla S\cdot \nabla +\mathrm{i}\case{1}{2} \Delta S\big]T(t), \qquad T(0)=I, \end{equation} \tag{ 14 }$$

where the operator in brackets on the right-hand side is self-adjoint, and hence T(t) is unitary. Hence, if we make for A(t) an ansatz

$$\begin{equation} A(t)=T(t)D(t)A_0, \end{equation} \tag{ 15 }$$

with an operator D(t) that is to be determined, and insert this into (12), then we obtain for D(t) the equation

$$\begin{equation} \mathrm{i}\partial _t D(t)=-\frac{\hbar }{2}\Delta (t)D(t) , \end{equation} \tag{ 16 }$$

where

$$\begin{equation} \Delta (t):=T^*(t)\Delta T(t), \end{equation} \tag{ 17 }$$

and with the initial condition D(0) = I. Since Δ(t) is self-adjoint, D(t) is unitary as well. From a more detailed analysis of the operator T(t) in section 4.2, we will learn that Δ(t) is a generalized Laplacian of the form

$$\begin{equation} \Delta (t)=\sum \alpha _{ij}(t,x)\partial _i\partial _j+\sum \beta _i(t,x)\partial _i, \end{equation} \tag{ 18 }$$

see (49). Collecting all terms we have a solution to the original Schrödinger equation of the form

$$\begin{equation} \psi (t)=[T(t)D(t)A_0]\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S(t)}. \end{equation} \tag{ 19 }$$

So the time evolution is described by three parts: (i) propagation of the Lagrangian manifold Λ_t, (ii) classical transport of the amplitude by the unitary operator T(t) and (iii) quantum dispersion described by D(t). Note that both T(t) and D(t) depend on the initial state via the initial phase function S₀.

The main condition we need in order that (19) holds on a time interval [0, T] is that the projection of Λ_t to position space is smooth for all t ∈ [0, T], i.e. Λ_t does not develop caustics. Under this condition, the representation (19) is exact, and S(t) and T(t) are both determined by transport along classical trajectories in phase space, the only contribution not yet linked to classical dynamics is the dispersive part D(t). We obtain the usual WKB result by neglecting the dispersion, i.e. by replacing D(t) by the identity I. To see what kind of error we make by doing this, we can use Duhamel's principle that we will also state in the general form for later use [20]. Let $\hat{H}(t)$ and $\hat{H}_1(t)$ be two time-dependent self-adjoint operators, and U(t) and U₁(t) be the time evolution operators generated by them with the initial conditions U(0) = U₁(0) = I; then,

$$\begin{equation} U(t)-U_1(t)=-\frac{\mathrm{i}}{\hbar }\int _0^t U(t)U^*(s)[\hat{H}(s)-\hat{H}_1(s)]U_1(s)\, \mathrm{d}s. \end{equation} \tag{ 20 }$$

This formula allows us to estimate how close the time evolutions generated by two different Hamiltonians are to each other. If we choose $\hat{H}(t)=-\frac{\hbar ^2}{2}\Delta (t)$ and $\hat{H}_1=0,$ then U(t) = D(t) and U₁(t) = I, and (20) gives

$$\begin{equation} D(t)A_0=A_0+\frac{\mathrm{i}\hbar }{2}\int _0^t D(t)D^*(t^{\prime })\Delta (t^{\prime })A_0 \, \mathrm{d}t^{\prime }. \end{equation} \tag{ 21 }$$

Since D is unitary, we then directly obtain

$$\begin{equation} \Vert D(t)A_0-A_0\Vert \le \frac{\hbar }{2}\int _0^t\Vert \Delta (t^{\prime })A_0\Vert \, \mathrm{d}t^{\prime }. \end{equation} \tag{ 22 }$$

So if the second-order derivatives of A₀ are bounded and the coefficients of Δ(t) are not growing with t, then the solution to Schrödinger's equation satisfies:

$$\begin{equation} \psi (t)=(T(t)A_0)\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S(t)}+O(\vert t\vert \hbar ), \end{equation} \tag{ 23 }$$

which is the standard time-dependent WKB result. Furthermore, Duhamel's formula can be iterated and gives an expansion with an error term of order O((ℏ|t|)^N) if the first 2N derivatives of A₀ are bounded and the coefficients of Δ(t) do not grow with t [20].

We would like to apply (19) to a coherent state, i.e. a state of the form (2). With S₀(x) = p · (x − q) and $A_0(x)=(\pi \hbar )^{-n/4}\,\mathrm{e}^{-\frac{1}{2\hbar }(x-q)\cdot B(x-q)}$, this state is of the form that we can use in (19), but now

$$\begin{equation} \Vert \Delta (t)A_0\Vert \sim 1/\hbar , \end{equation} \tag{ 24 }$$

so (22) does not allow us to use the standard WKB approximation D(t) ≈ I.

The way to solve this problem is to approximate the action of D(t) on A₀ not by the identity, as in the WKB method, but to borrow from the standard propagation of coherent states and approximate the generator of D(t) by its Taylor expansion around the centre of A₀ up to second order. Since D(t) acts on a state that is concentrated in position space at x = q and in momentum space at p = 0, this means that we freeze the coefficients of Δ(t) at x = q and approximate D(t) by the operator generated by it.

To formalize this idea, we introduce for $q\in \mathbb{R}^n$ the unitary operator L_q as

$$\begin{equation} (L_qa)(x):=\hbar ^{-n/4}a((x-q)/\sqrt{\hbar }), \end{equation} \tag{ 25 }$$

which shifts by q and then rescales in ℏ, so that (L_qa)(x) is concentrated around x = q, e.g., if a(x) = π^−n/4 e^−xBx/2, then A₀ = L_qa is the amplitude of the coherent state (2) considered above. Now inserting A = L_qa into (16) gives for a the equation

$$\begin{equation} \mathrm{i}\partial _t a=\hbar L_q^*\Delta (t)L_q a, \end{equation} \tag{ 26 }$$

and we will approximate the operator on the right-hand side by

$$\begin{equation} \Delta _q(t):=\lim _{\hbar \rightarrow 0}\hbar L_q^*\Delta (t)L_q . \end{equation} \tag{ 27 }$$

From (18), one finds that Δ_q(t) = ∑α_ij(t, q)∂_i∂_j and

$$\begin{equation} \hat{Q}_q(t):=\hbar L_q^*\Delta (t)L_q-\Delta _q(t)=O_t(\sqrt{\hbar }) \end{equation} \tag{ 28 }$$

in the sense that $\hat{Q}_q(t)a(x)=O_t(\sqrt{\hbar })$ if a is smooth and has bounded derivatives. So the operator Δ_q(t) is indeed obtained from Δ(t) by freezing the coefficients at x = q and furthermore discarding the first-order terms. In (52), we will obtain more explicit bounds on (28) in terms of the classical flow.

Therefore, if we define the operator D_q(t) by

$$\begin{equation} \mathrm{i}\partial _t D_q(t)=-\case{1}{2}\Delta _q(t)D_q(t),\qquad D_q(0)=I, \end{equation} \tag{ 29 }$$

then we expect that D(t)L_qa ≈ L_qD_q(t)a holds. To check this, we use again Duhamel's principle (20)

$$\begin{equation} D(t)L_qa= L_qD_q(t)a+\frac{\mathrm{i}}{2} \int _0^t D(t)D^*(s)L_q\hat{Q}_q(s) D_q(s) a\, \mathrm{d}s \end{equation} \tag{ 30 }$$

and, since D(t) and L_q are unitary, we obtain from (28)

$$\begin{equation} \Vert D(t)L_qa-L_qD_q(t)a\Vert \le \int _0^t \Vert \hat{Q}_q(s)D_q(s)a\Vert \, \mathrm{d}s =O_t(\sqrt{\hbar }). \end{equation} \tag{ 31 }$$

The t-dependence in the remainder term will be governed by the behaviour of $\hat{Q}_q(t)$. In the next section, we give conditions on the initial state and on the classical dynamics under which the coefficients of $\hat{Q}_q(t)$ stay bounded, see (51) and (52).

To summarize our results so far, if $\psi _0=A_0\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S_0}$, where A₀ = L_qa is concentrated around q, then the time-evolved state is given by

$$\begin{equation} \psi (t)=[T(t)L_qD_q(t)a]\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S(t)}+O_t(\sqrt{\hbar }), \end{equation} \tag{ 32 }$$

so we can use the standard time-dependent WKB to propagate coherent states if we include the additional operator D_q(t). This operator has a quadratic generator and is therefore a metaplectic operator, see [3, 32, 33], and so we call it the metaplectic correction to the standard time-dependent WKB method. In order to understand the behaviour of D_q(t) in some more detail, we have to look at the propagation of the Lagrangian manifold Λ by the classical dynamics, and the associated geometric interpretation of our extended time-dependent WKB scheme. This will be the subject of the next section.

We finally note that it is sometimes useful to work with M_q(t) := L_qD_q(t)L*_q, because then L_qD_q(t)a = M_q(t)L_qa = M_q(t)A₀ and we can allow for A₀ to be of a more general form than the one induced by the specific scaling with L_q, as long as it is concentrated around x = q. Note that M_q satisfies the equation

$$\begin{equation} \mathrm{i}\hbar \partial _t M_q(t)=-\frac{\hbar ^2}{2}\Delta _q(t) M_q(t),\quad M_q(0)=I, \end{equation} \tag{ 33 }$$

and the time evolution of an initial state $\psi _0=A_0\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S_0}$, where A₀(x) is concentrated around x = q, can be approximated by

$$\begin{equation} \psi (t)= [T(t) M_q(t)A_0]\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S(t)}+O_t(\sqrt{\hbar }). \end{equation} \tag{ 34 }$$

In order to understand the scope and the accuracy of the metaplectic extension of time-dependent WKB theory, which we outlined in the previous section, we have to understand some of the geometry underlying it. We will see that in particular the time dependence of the remainder term is governed by the way in which the classical flow transports the initial Lagrangian manifold Λ. Furthermore, a proper understanding of the geometry will allow us to extend the scheme beyond caustics.

4.1. Transport of the Lagrangian manifold

Let us first return to the solutions of the Hamilton–Jacobi equation (11), which we will discuss for a general Hamiltonian H(ξ, x):

$$\begin{equation} \partial _t S(t,x)+H(\nabla S(t,x),x)=0. \end{equation} \tag{ 35 }$$

As is well known, there are two main ingredients involved in the solution; see, e.g., [5, 30, 31]. The first is the transport of the initial Lagrangian manifold

$$\begin{equation} \Lambda :=\lbrace (\nabla S(x),x)\, :\, x\in U\rbrace , \end{equation} \tag{ 36 }$$

where $U\subset \mathbb{R}^n$ is an open set that contains the support of the amplitude A(x). Let us denote by Φ^t(ξ, x), where $(\xi ,x)\in \mathbb{R}^n\times \mathbb{R}^n$, the Hamiltonian flow generated by H, i.e. Φ^t(ξ, x) = (p(t), q(t)), where (p(t), q(t)) are the solutions to Hamilton's equation $\dot{p}=-\nabla _qH(p,q)$ and $\dot{q}=\nabla _p H(p,q)$ with initial conditions p(t = 0) = ξ and q(t = 0) = x, respectively. Using this flow, we can transport the Lagrangian manifold Λ and obtain a family of Lagrangian manifolds

$$\begin{equation} \Lambda _t=\Phi ^t(\Lambda )=\lbrace \Phi ^t(\nabla S(x),x)\, :\, x\in U\rbrace . \end{equation} \tag{ 37 }$$

We will call the initial Lagrangian manifold Λ non-contracting (with respect to the flow Φ^t) if there exists a C > 0 such that for all t > 0

$$\begin{equation} \Vert \mathrm{d}\Phi ^t|_{\Lambda }(z)\Vert \ge C ,\quad {\rm for\ all}\quad z\in \Lambda . \end{equation} \tag{ 38 }$$

Here, $\mathrm{d}\Phi ^t|_{\Lambda }(z):T_z\Lambda \rightarrow T_{\Phi ^t(z)}\Lambda _t$ is the restriction of the linearized flow at z = (p, q) to Λ, which is given in local coordinates by the matrix of the derivatives of the components of Φ^t with respect to the coordinates on Λ. This condition means that, roughly speaking, trajectories starting nearby on Λ do not coalesce into each other. An example for a non-contracting submanifold is the unstable manifold of a hyperbolic trajectory, and more generally, if a system is hyperbolic, then any submanifold that is transversal to the stable foliation is non-contracting. Stable manifolds are then of course examples of manifolds which are not non-contracting. In an integrable system, any manifold that lies in the regular part of the foliation of phase space into invariant tori is non-contracting.

The second ingredient needed to determine S(t, x) is the projection of Λ_t to position space along the momentum directions, i.e. we take the projection $\pi :\mathbb{R}_p^n\times \mathbb{R}_q^n\rightarrow \mathbb{R}_q^n$ defined by π(p, q) = q and restrict it to Λ_t. If this map,

$$\begin{equation} \pi _{\Lambda _t}:\Lambda _t\rightarrow \mathbb{R}^n, \end{equation} \tag{ 39 }$$

has no singularities, i.e. π(Λ_t) = U_t does not contain any caustics of Λ_t, then there exists, at least locally, a phase function S(t, x) such that

$$\begin{equation} \Lambda _t=\lbrace (\nabla S(t,q), q)\, :\, q\in U_t\rbrace , \end{equation} \tag{ 40 }$$

for some open set $U_t\in \mathbb{R}^n$. By assumption, $\pi _{\Lambda _0}:\Lambda _0\rightarrow U_0\subset \mathbb{R}^n$ is smooth so we can invert it, with inverse given by $\pi _{\Lambda _0}^{-1}(x)=(\nabla S(0,x),x)$, and hence, we can form

$$\begin{equation} \phi _{\Lambda }(t,x):=\pi _{\Lambda _t}\Phi ^t\big (\pi _{\Lambda _0}^{-1}(x)\big ), \end{equation} \tag{ 41 }$$

which is a map from U → U_t, and if $\pi _{\Lambda _t}$ is smooth, it is a smooth and invertible map. But note that ϕ_Λ(t) is not a flow. Analogously to the above notion of non-contracting, (38), we will say that ϕ_Λ(t) is non-contracting on $U\subset \mathbb{R}^n$ if there exist C > 0 such that

$$\begin{equation} \Vert \phi _{\Lambda }^{\prime }(t,x)\Vert \ge C,\quad {\rm for\ all}\quad x\in U, \end{equation} \tag{ 42 }$$

independent of t ⩾ 0, where ϕ'_Λ(t, x) denotes the matrix of first-order derivatives of the components of ϕ_Λ(t, x). A necessary condition for ϕ_Λ to be non-contracting is that Λ₀ is non-contracting, but in addition we need that the projection $\pi :\Lambda _t\rightarrow \mathbb{R}^n$ is not singular. This implies that we have to stay away from caustics.

4.2. The transport operator

Another description of the map ϕ_Λ(t, x) is as follows: take the trajectory (p(t), q(t)) = Φ^t(ξ, x), which starts at t = 0 at q(t = 0) = x with the initial momentum ξ = ∇S(x); then, the map ϕ_Λ(t, x) is the q component of this trajectory, ϕ_Λ(t, x) = q(t). Since by Hamilton's equations $\dot{q}=\nabla _p H(p,q)$ and on Λ_t we have p = ∇S(t, q), we find that ϕ_Λ(t, x) is the unique solution to

$$\begin{equation} \frac{\mathrm{d}\phi _{\Lambda }(t,x)}{\mathrm{d}t}=\nabla _pH(\nabla S(t,\phi _{\Lambda }(t,x)),\phi _{\Lambda }(t, x)),\quad {\rm with}\quad \phi _{\Lambda }(0,x)=x. \end{equation} \tag{ 43 }$$

Note that we allowed here for a general Hamiltonian, and not only one of the simple form kinetic plus potential energy.

We will now show that the transport operator T(t) defined in (14) can be written using this map as

$$\begin{equation} (T(t)A)(x)=\big[\det \phi _{\Lambda }^{\prime }\big(t,\phi _{\Lambda }^{-1}(t,x)\big)\big]^{-1/2} A\big(\phi _{\Lambda }^{-1}(t,x)\big). \end{equation} \tag{ 44 }$$

By comparing this formula with appendix A.2 and using (43), we see that the operator in (44) satisfies

$$\begin{equation} \mathrm{i}\hbar \partial _t T(t)=\hat{K}(t) T(t), \end{equation} \tag{ 45 }$$

where $\hat{K}(t)$ is the Weyl quantization of the function

$$\begin{equation} K(t;p,q)=\nabla _pH(\nabla S(t,q),q)\cdot p. \end{equation} \tag{ 46 }$$

In particular, we see that for the special case $H=\frac{1}{2}p^2+V(q),$ we have from appendix A.2

$$\begin{equation} \hat{K}(t)=\frac{\hbar }{\mathrm{i}}\bigg [\nabla S(t,x)\nabla +\frac{1}{2}\Delta S(t,x) \bigg ] \end{equation} \tag{ 47 }$$

and comparing with the transport equation (12) proves our claim.

4.3. The metaplectic correction

We can use this observation to give a more detailed description of the operator

$$\begin{equation} {-}\frac{\hbar ^2}{2}\Delta (t)=-\frac{\hbar ^2}{2} T^*(t)\Delta T(t). \end{equation} \tag{ 48 }$$

This is an application of Egorov's theorem and the technical details are given in appendix A.3.2, where we show that the Weyl symbol of this operator is given by

$$\begin{equation} \delta H(t,\xi ,x)=\frac{1}{2}\xi \cdot {\mathcal A}(t,x) \xi +\frac{\hbar ^2}{16}\sum _i \mathop {\mathrm{Tr}}\nolimits \big[\big(\phi _{\Lambda }^{-1}\big)_i^{\prime \prime }(t,x)\big]^2, \end{equation} \tag{ 49 }$$

where ${\mathcal A}(t,x)$ is a symmetric matrix given by

$$\begin{equation} {\mathcal A}(t,x)=(\phi _{\Lambda }^{\prime }(t,x) \phi _{\Lambda }^{\prime }(t,x)^{\dagger })^{-1}, \end{equation} \tag{ 50 }$$

and (ϕ⁻¹_Λ)_i'' is the matrix of second derivatives of the ith component of ϕ⁻¹_Λ. Bounds on the coefficients of the operator Δ(t) play an important role in estimating the accuracy of the time-dependent WKB propagation and the metaplectic extension, because Δ(t) appears in the error terms (22) and (31). This is where the non-contraction condition that we introduced above becomes important; from (42), we have that if ϕ_Λ(t) is non-contracting, then

$$\begin{equation} \Vert {\mathcal A}(t,x)\Vert \le C, \end{equation} \tag{ 51 }$$

and hence, the operator Δ(t) has bounded coefficients. We can now also give a more explicit description of the operator Δ_q(t). A short calculation shows that for a phase space function H(ξ, x), which is quadratic in the momentum ξ, we have $\hbar L_q^*\hat{H}L_q =\hat{H}_q$, with $H_q(\xi , x)=H(\xi , q+\sqrt{\hbar }\, x)$ and so $H_q(\xi , x)=H_q(\xi , q)+O(\sqrt{\hbar })$. Applying this to $\delta H(t,\xi ,x)=\frac{1}{2}\xi {\mathcal A}(t,x) \xi +O(\hbar ^2)$ gives that

$$\begin{equation} \Delta _q(t)=\nabla _x\cdot {\mathcal A}(t,q)\nabla _x\qquad {\rm and}\qquad \hbar L^*_q\Delta (t)L_q=\Delta _q(t)+O_t(\sqrt{\hbar }). \end{equation} \tag{ 52 }$$

If ϕ_Λ(t) is non-contracting, then the remainder is bounded uniformly in time.

Since we have now an explicit expression for Δ_q(t), we can compute the action of the metaplectic operator D_q(t), defined in (29), on a function a. Set

$$\begin{equation} C_t:=\int _0^t {\mathcal A}(s,q)\, \mathrm{d}s \end{equation} \tag{ 53 }$$

and let $\hat{a}(\xi )=\int \mathrm{e}^{-\mathrm{i}x\xi } a(x)\, \mathrm{d}x$ be the Fourier transform of a, and then,

$$\begin{equation} (D_q(t)a)(x)=\frac{1}{(2\pi )^n}\int \mathrm{e}^{-\frac{\mathrm{i}}{2} \xi \cdot C_t\xi }\hat{a}(\xi )\, \mathrm{e}^{\mathrm{i}x\cdot \xi }\, \mathrm{d}\xi . \end{equation} \tag{ 54 }$$

Similarly, we find for the rescaled operator M_q(t) that integrating (33) gives

$$\begin{equation} (M_q(t)A)(x)=\frac{1}{(2\pi \hbar )^n}\int \mathrm{e}^{-\frac{\mathrm{i}}{\hbar }\frac{1}{2} \xi \cdot C_t\xi }\hat{A}_{\hbar }(\xi )\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar } x\cdot \xi }\, \mathrm{d}\xi , \end{equation} \tag{ 55 }$$

where $\hat{A}_{\hbar }(\xi )=\int \mathrm{e}^{-\frac{\mathrm{i}}{\hbar } x\xi } A(x)\, \mathrm{d}x$. The operators D_q(t) and M_q(t) act as the Fourier multipliers with Gaussian functions. For an initial amplitude of the form $A(x)=(\pi \hbar )^{-n/4}\,\mathrm{e}^{-\vert x\vert ^2/(2\hbar )}$, one finds in particular

$$\begin{equation} M_q(t)A(x)=\frac{1}{(\pi \hbar )^{n/4}}\frac{1}{\sqrt{\det (I+\mathrm{i}C_t)}}\, \mathrm{e}^{-\frac{1}{2\hbar }x\cdot (I+\mathrm{i}C_t)^{-1}x}. \end{equation} \tag{ 56 }$$

4.4. The role of the initial manifold

The localized initial states that we consider are of the form $\psi =(L_q a)\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar } S}$, and the phase function determines the initial Lagrangian manifold, which is crucial for the semiclassical propagation scheme that we presented. The question we want to consider now is if the state ψ determines the Lagrangian manifold uniquely, or in other words, if there might be other functions $\tilde{a}$ and $\tilde{S}$, such that $\psi =(L_q\tilde{a})\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }\tilde{S}}$. In that case, we would have some freedom in the choice of the initial Lagrangian manifold we use for propagation.

The condition $(L_q a)\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar } S}=(L_q\tilde{a})\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }\tilde{S}}$ gives after multiplication by $\mathrm{e}^{-\frac{\mathrm{i}}{\hbar }\tilde{S}}$ and application of L*_q that

$$\begin{equation} \tilde{a}(x)=a(x)\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }[S(q+\sqrt{\hbar }\, x)-\tilde{S}(q+\sqrt{\hbar }\, x) ]}, \end{equation} \tag{ 57 }$$

and by the Taylor expansion, we see that

$$\begin{eqnarray} &&\fl \frac{\mathrm{i}}{\hbar }[S(q+\sqrt{\hbar }\, x)-\tilde{S}(q+\sqrt{\hbar }\, x)]=\frac{\mathrm{i}}{\hbar }[S(q)-\tilde{S}(q)]+\frac{\mathrm{i}}{\sqrt{\hbar }}[\nabla S(q)-\nabla \tilde{S}(q)]x +\mathrm{i}R(\hbar , q,x),\nonumber\\ \end{eqnarray} \tag{ 58 }$$

where R(ℏ, q, x) is a smooth real-valued function. The first term on the right-hand side gives just a constant phase factor, so if

$$\begin{equation} \nabla S(q)=\nabla \tilde{S}(q), \end{equation} \tag{ 59 }$$

then $\tilde{a}$ defines a nice smooth function. But this condition means that the two Lagrangian submanifolds Λ and $\tilde{\Lambda }$, generated by S and $\tilde{S}$, respectively, intersect at (p, q) = (∇S(q), q). We therefore conclude that in order to propagate a state localized around a phase space point (p, q), we can use any Lagrangian manifold through that point that is non-contracting. This freedom of choice can be used to select, for instance, initial manifolds for which the propagation becomes particularly easy, e.g., some for which no caustics develop.

In the case that the trajectory through (p, q) is hyperbolic, any Lagrangian submanifold that is transversal to the stable directions is non-contracting, and hence can be used for propagation. Although it seems most natural to take the unstable manifold, any other manifold that is sufficiently transversal to stable direction will converge exponentially fast to the unstable manifold and hence should work with comparable efficiency. We test this with the example in section 7.

4.5. Relation to the standard coherent state propagation and time scales

For times that are short compared to the Ehrenfest time, our approach should reproduce the standard results on coherent state propagation. For the general class of initial amplitudes of the form (25), the propagation is reviewed in [34]. A propagated coherent state is of the form

$$\begin{equation} \psi (t,x)=\mathrm{e}^{\frac{\mathrm{i}}{\hbar }l(t)} \frac{1}{\hbar ^{n/4}} a_{\hbar }\bigg (t,\frac{x-q(t)}{\sqrt{\hbar }}\bigg )\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar } p(t)\cdot (x-q(t)}, \end{equation} \tag{ 60 }$$

where (p(t), q(t)) is a classical trajectory, l(t) is an action-type phase that also includes the Maslov terms, and $a_{\hbar }(t,x)=a_0(t,x)+\sqrt{\hbar }\, a_{1}(t,x)+\cdots$ with leading term given as the solution to

$$\begin{equation} \mathrm{i}\partial _ta_0(t,x)=\big[{-}\case{1}{2}\Delta + \case{1}{2} x\cdot V^{\prime \prime }(q(t)) x\big]a_0(t,x),\qquad a_0(0,x)=a_0. \end{equation} \tag{ 61 }$$

This means that the centre of the state is propagated along the trajectory (p(t), q(t)) and its shape changes according to the linearized dynamics around the centre.

If we expand S(t, x) around x = q(t) up to second order, then we find that we can write our approximation (32) in the form (60) with

$$\begin{equation} \fl a_{\hbar }(t,x)= b(t,x)\,\mathrm{e}^{\frac{\mathrm{i}}{2}x\cdot S^{\prime \prime }(t,q(t))x+O(\sqrt{\hbar })},\qquad {\rm where}\quad b(t,x)=[L^*_{q(t)} T(t)L_qD_q(t)a](x). \end{equation} \tag{ 62 }$$

Now, using the operator identities

$$\begin{eqnarray} &&\mathrm{i}\partial _t L^*_{q(t)}= L^*_{q(t)}\mathrm{i}\partial _t+L^*_{q(t)} \mathrm{i}\dot{q}\cdot \nabla , \end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray} &&\sqrt{\hbar }L^*_{q(t)}\nabla L_{q(t)}=\nabla , \end{eqnarray} \tag{ 64 }$$

$$\begin{eqnarray} &&L^*_{q(t)} \nabla S(t,x)L_{q(t)}=\nabla S(t,q(t)) +\sqrt{\hbar } S^{\prime \prime }(t, q(t)) x+O(\hbar ) , \end{eqnarray} \tag{ 65 }$$

$$\begin{eqnarray} &&L_{q(t)}^*T(t)L_q \Delta _qL^*_qT^*(t)L_{q(t)}=\Delta +O(\sqrt{\hbar }) , \end{eqnarray} \tag{ 66 }$$

which all follow by direct computation, we obtain with $\nabla S(t,q(t))=p(t)=\dot{q}(t)$ that

$$\begin{equation} \mathrm{i}\partial _t b(t,x)=-\bigg (\mathrm{i}x\cdot S^{\prime \prime }\nabla +\frac{\mathrm{i}}{2}\Delta S +\frac{1}{2}\Delta \bigg ) b(t,x)+O(\sqrt{\hbar }), \end{equation} \tag{ 67 }$$

where S and its derivatives are all evaluated at x = q(t). If we furthermore use $\Delta ( b(t,x)\mathrm{e}^{\frac{\mathrm{i}}{2}x\cdot S^{\prime \prime }x})=[\Delta b(t,x)+2\mathrm{i}\nabla b(t,x)\cdot S^{\prime \prime } x + (-\vert S^{\prime \prime }x\vert ^2+\mathrm{i}\Delta S)b(t,x)]\,\mathrm{e}^{\frac{\mathrm{i}}{2}x\cdot S^{\prime \prime }x}$ and combine this with (67), we find

$$\begin{equation} \mathrm{i}\partial _t a_h(t,x)=\big[{-}\case{1}{2}\Delta -\case{1}{2}x\cdot \dot{S}^{\prime \prime } x-\case{1}{2}\vert S^{\prime \prime } x\vert ^2\big]a_h(t,x)+O(\sqrt{\hbar }). \end{equation} \tag{ 68 }$$

Finally, from the Hamilton–Jacobi equation (11), we obtain $-x\cdot \dot{S}^{\prime \prime } x-\vert S^{\prime \prime } x\vert ^2/2=x\cdot V^{\prime \prime } x$, and hence, the leading term a₀(t, x) satisfies (61).

This was a formal computation that showed that our result reproduces the previously existing results, but did not gave much insight as to where the standard approximation might break down. To see this more clearly let us just look at the amplitude T(t)L_qb with b = D_qa. By the Taylor expansion around x = q(t) = ϕ_Λ(t, q), we have ϕ⁻¹_Λ(x) − q ≈ [ϕ_Λ'(t, q(t))]⁻¹(x − q(t)), and hence,

$$\begin{equation} T(t)L_q b(x)\approx \frac{1}{\hbar ^{\frac{n}{4}}\sqrt{\det \phi _{\Lambda }^{\prime }(t,q(t))}} \, b( [\sqrt{\hbar }\, \phi _{\Lambda }^{\prime }(t,q(t))]^{-1}(x-q(t))). \end{equation} \tag{ 69 }$$

From this expression, we see that the state will be localized around x = q(t) if

$$\begin{equation} \sqrt{\hbar }\Vert \phi _{\Lambda }^{\prime }(t,q(t))\Vert \ll 1, \end{equation} \tag{ 70 }$$

and since ϕ_Λ(t, x) is a projection of the flow from phase space, this is satisfied if t ≪ T_E(ℏ). On the other hand, if $\sqrt{\hbar }\Vert \phi _{\Lambda }^{\prime }(t,q(t))\Vert \approx 1$, then the amplitude is no longer localized around a point, and the state becomes a Lagrangian or a WKB state.

So far we have avoided the discussion of caustics, but although there are situations where no caustics occur, in many interesting situations we expect caustics to develop. A caustic is the image in position space $\mathbb{R}^n$ of the singularities of the projection $\pi _{\Lambda _t}:\Lambda _t\rightarrow \mathbb{R}^n$, which means that at a caustic the tangent plane T_zΛ_t to Λ_t at a point z = (p, q) ∈ Λ_t becomes vertical in the sense that T_zΛ_t∩V_z ≠ {0}, where $V_z=\lbrace (\xi ,q)\,\, ;\, \xi \in \mathbb{R}^n\rbrace$ denotes the vertical subspace. In a neighbourhood of a caustic, the manifold Λ_t can no longer be represented by a generating function S(t, x) as in (40), and so the simple time-dependent WKB propagation we presented above can no longer work. This problem is usually solved by switching to a different representation of the quantum state, e.g., taking the Fourier transform switches from the position representation to the momentum representation, and in this new representation, we have to consider instead the projection of Λ_t to momentum space.

We want to use a method that allows for a bit of flexibility so that we can accommodate all possible different representations. To this end, we look at how we can in general approximate the full Hamiltonian H by a simpler Hamiltonian H₁ in a way that the propagation of Lagrangian states initially associated with Λ can be approximated using H₁. The crucial observation that will guide our conditions on H₁ is that if B(p, q) and A(x) are sufficiently smooth functions and $\hat{B}$ is the Weyl quantization of B, then

$$\begin{equation} \hat{B}(A\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S})=(B_{\Lambda }A+O(\hbar ))\, \mathrm{e}^{\frac{\mathrm{i}}{\hbar }S}, \end{equation} \tag{ 71 }$$

where B_Λ(x) = B(∇S(x), x); this is a classical result, see, e.g, [31]. This means that the Lagrangian state is concentrated in phase space on Λ; in particular, if B_Λ = 0, then $\hat{B} \psi =O(\hbar )$. More generally, if B vanishes of order N on Λ, i.e. (∂^αB)|_Λ = 0 for |α| ⩽ N, then

$$\begin{equation} \hat{B}\psi =O(\hbar ^N). \end{equation} \tag{ 72 }$$

Since a propagated Lagrangian state is concentrated near the propagated Lagrangian manifold Λ_t = Φ^t(Λ), where Φ^t is the Hamiltonian flow generated by H, we expect that if H₁ is close to H near Λ_t, then it will provide an accurate approximation for the purpose of propagating states concentrated near Λ_t. The conditions for a function H₁ to be a general first-order approximation of H near Λ_t are

$$\begin{equation} (H-H_1(t))|_{\Lambda _t}=0\qquad {\rm and}\qquad (\nabla H-\nabla H_1(t))|_{\Lambda _t}=0 . \end{equation} \tag{ 73 }$$

These conditions ensure that the Hamiltonian vector fields generated by H and H₁ agree on Λ_t; hence, if we denote by Φ^t₁ the time t map generated by H₁, then

$$\begin{equation} \Phi ^t|_{\Lambda }=\Phi _1^t|_{\Lambda }, \end{equation} \tag{ 74 }$$

i.e. the classical dynamics agree when restricted to Λ.

Let us now turn to the quantizations of H and H₁. Let U(t) and U₁(t) be the time evolution operators generated by $\hat{H}$ and $\hat{H}_1$, respectively; then, we expect from (74) that U(t)ψ ≈ U₁(t)ψ for a Lagrangian state associated with Λ. To verify this, we make an ansatz

$$\begin{equation} U(t)=U_1(t)V(t), \end{equation} \tag{ 75 }$$

and inserting this into the Schrödinger equation for U(t) gives for V(t) the equation

$$\begin{equation} \mathrm{i}\hbar \partial _t V(t)=\widehat{\delta H}(t) V(t),\qquad {\rm with}\quad V(0)=I, \end{equation} \tag{ 76 }$$

where

$$\begin{equation} \widehat{\delta H}(t)=U_1^*(t)[\hat{H}-\hat{H}_1(t)]U_1(t). \end{equation} \tag{ 77 }$$

By Egorov's theorem, see appendix A.3, the operator $\widehat{\delta H}(t)$ is, in leading order in ℏ, the quantization of H − H₁ transported along the map Φ^t₁:

$$\begin{equation} \delta H(t)=[H-H_1(t)]\circ \Phi _1^t+O_t(\hbar ^2). \end{equation} \tag{ 78 }$$

This implies from (73) that

$$\begin{equation} \delta H(t)|_{\Lambda _0}=O_t(\hbar ^2) \qquad {\rm and}\qquad \nabla \delta H(t)|_{\Lambda _0}=O_t(\hbar ^2), \end{equation} \tag{ 79 }$$

and therefore, from (72), we have for an initial state of the form $\psi =A_0\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S_0}$, where A₀ has bounded derivatives, that

$$\begin{equation} \Vert \widehat{\delta H}(t)\psi \Vert =O_t(\hbar ^2). \end{equation} \tag{ 80 }$$

Thus, the generator of V(t) is small when applied to ψ, and therefore, Duhamel's principle (20) gives

$$\begin{equation} V(t)\psi =\psi +O_t(\hbar ), \end{equation} \tag{ 81 }$$

and hence, U₁(t)ψ is a good approximation for U(t)ψ if ψ is a Lagrangian state.

This is a generalization of the standard time-dependent WKB method presented in section 2. The choice of H₁ corresponding to the results in section 2 is

$$\begin{equation} H_1(t;\xi ,x):=H(\nabla S(t,x),x)+\nabla _{\xi }H(\nabla S(t,x),x)\cdot \xi , \end{equation} \tag{ 82 }$$

which we obtain from H(p, q) = H(∇S(t, x) + ξ, x) as a first-order approximation in the shifted momentum ξ = p − ∇S(t, x). For a Hamiltonian of the form $H=\frac{1}{2}p^2+V(q),$ one finds using appendix A.2 that $\hat{H}_1=H(\nabla S(t,x),x)-\mathrm{i}\hbar \nabla S(t,x)\cdot \nabla -\frac{\mathrm{i}\hbar }{2}\Delta S(t,x)$, and so by comparing this with (14) and (11), we see that in this case we have

$$\begin{equation} U_1=\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S(t)}T(t)\mathrm{e}^{-\frac{\mathrm{i}}{\hbar }S_0}\,\, \quad {\rm and}\quad V(t)=\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S_0}D(t)\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar }S_0}. \end{equation} \tag{ 83 }$$

So how does a caustic affect this relation? In defining H₁ we used coordinates (ξ, x) near Λ_t defined by p = ∇S(t, x) + ξ and q = x, so that Λ_t = {ξ = 0}, but at a caustic this coordinate system becomes singular. In order to repair this we can just use a different canonical coordinate system near Λ_t, one which does not become singular, and choose for H₁ the first-order Taylor approximation in the transversal coordinates.

The only problem we might encounter with this more general version of time-dependent WKB theory is that the control of the time dependence of the remainder term (80) becomes more involved. In section 2, we were able to take advantage of the very explicit version of Egorov's theorem from appendix A.3.2 to reduce this problem to the non-contractiveness of the transport map ϕ_Λ(t). For different choices of H₁, this problem can become much harder and we reserve a more thorough investigation for the future. But this question concerns our ability to obtain explicit error estimates; it does not prevent us from using, after suitable testing, the method for explicit numerical propagation in concrete problems.

Now, if the initial state ψ₀ is localized in phase space around a point (p, q) ∈ Λ₀ so that the estimate (80) no longer holds, then it is natural to approximate V(t) by taking only the behaviour of its generator δH(t, ξ, x) near (p, q) into account. From (73), we have δH(t, p, q) = 0 and ∇δH(t, p, q) = 0, and hence, the simplest nontrivial approximation is quadratic, with z = (ξ, x) and z₀ = (p, q), and we have

$$\begin{equation} \delta H(t,x,\xi )\approx \delta H_2(t,z)=\case{1}{2}(z-z_0)\delta H^{\prime \prime }(t,z_0)(z-z_0). \end{equation} \tag{ 84 }$$

The time evolution $\tilde{M}_{p,q}(t)$ generated by $\widehat{\delta H}_2(t)$ via $\mathrm{i}\hbar \partial _t M_{p,q}(t)=\widehat{\delta H}_2(t) M_{p,q}(t)$ and M_{p, q}(0) = I is a metaplectic operator since the generator $\widehat{\delta H}(t)$ is the quantization of a quadratic function on phase space. Therefore, for ψ strongly localized around (p, q), we expect V(t)ψ ≈ M_{p, q}ψ, and hence, we have arrived now at the more general version of (34), which reads

$$\begin{equation} U(t)\psi = U_1(t) M_{p,q}(t)\psi +O_t(\sqrt{\hbar }) \end{equation} \tag{ 85 }$$

if ψ is concentrated in phase space at (p, q) as in (25).

The relation of M_{p, q}(t) to M_q(t) follows from (83); if H₁ is of the form (82), we find that

$$\begin{equation} M_{p,q}(t)=\mathrm{e}^{\frac{\mathrm{i}}{\hbar } S_0^{(2)}} M_q(t)\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar } S_0^{(2)}}, \end{equation} \tag{ 86 }$$

where p = ∇S₀(q) and $S_0^{(2)}(x)=p(x-q)+\frac{1}{2}(x-q)S_0^{\prime \prime }(q)(x-q)$ is the quadratic part of the Taylor expansion of S₀ around x = q.

Since M_{p, q}(t) has a quadratic generator, it is the quantization of a linear symplectic map P_{p, q}(t) on $T_{p,q}(\mathbb{R}^n\times \mathbb{R}^n),$ and by construction, this map can be expressed in terms of the linearizations of the maps Φ^t and Φ^t₁. Since V(t) = U⁻¹₁(t)U(t), we can view V(t) as a quantization of the map (Φ^t₁)⁻¹○Φ^t, and M_{p, q}(t) is then the quantization of the linearization of that map around p, q:

$$\begin{equation} P_{p,q}(t)=\big(\mathrm{d}\Phi _1^{t}\big)^{-1}\,\mathrm{d}\Phi ^t. \end{equation} \tag{ 87 }$$

Since Φ^t and Φ^t₁ are identical on Λ, we have that $P_{p,q}(t)|_{T_{p,q}\Lambda _0}=I$, and hence, P_{p, q} is a shear relative to the tangent space T_{p, q}Λ₀ of the initial Lagrangian manifold at (p, q).

In case that H₁ is given in (82), the map P_{p, q}(t) can be described in some more detail. Note that since H₁ in (82) contains only a linear term in the momentum, we find that Hamilton's equations generated by H₁ give for x(t) the equation

$$\begin{equation} \dot{x}=\nabla _pH(\nabla S(t,x),x), \end{equation} \tag{ 88 }$$

and since the right-hand side does not depend on the momentum ξ, the trajectory x(t) does not depend on the initial momentum ξ. Hence, Φ^t₁ maps vertical subspaces into vertical subspaces, i.e. $\mathrm{d}\Phi ^t_1 (V_{z_0})\subset V_{z(t)}$, where for z = (p, q), we set $V_z=\lbrace (\xi ,q)\, ,\xi \in \mathbb{R}^n\rbrace$. Therefore, the map P_{p, q}(t) satisfies

$$\begin{equation} P_{p,q}(t)^{-1}|_{T_{p,q}\Lambda _0}=I\quad {\rm and}\quad P_{p,q}(t)^{-1}(V_{(p,q)})=(\mathrm{d}\Phi ^t)^{-1}(V_{\Phi ^t(p,q)}), \end{equation} \tag{ 89 }$$

and as we show in appendix B, this implies that the knowledge of the Lagrangian subspace $(\mathrm{d}\Phi ^t)^{-1}(V_{\Phi ^t(p,q)})$ determines the map P_{p, q}(t)⁻¹ and hence P_{p, q}(t) uniquely.

Let us look at two examples where we can determine the long time limit of $(\mathrm{d}\Phi ^t)^{-1}(V_{\Phi ^t(p,q)})$ from dynamical conditions.

• (a)  
  Assume we have a one-dimensional system with an x-independent Hamiltonian H(p), which satisfies ∂²H/∂p² > 0. Then, the velocity increases with p and hence $(\mathrm{d}\Phi ^t)^{-1}(V_{\Phi ^t(p,q)})\rightarrow \lbrace (0,q), q\in \mathbb{R}\rbrace$ for t → ∞. Therefore, if $T_{(p,q)}\Lambda _0\ne \lbrace (0,q), q\in \mathbb{R}\rbrace ,$ we can use the result from appendix B and see that for large t the map P_{p, q}(t) tends to a limit P^∞_{p, q}, and hence, the corresponding metaplectic correction M_{p, q}(t) tends to a limit, too.
• (b)  
  If the trajectory through (p, q) is hyperbolic and the vertical subspaces V_z are transversal to the unstable subspaces V^u_z along the trajectory, then $(\mathrm{d}\Phi ^t)^{-1}(V_{\Phi ^t(p,q)})$ tends to the stable subspace V^s_{p, q} for large t, and at an exponential rate if the hyperbolicity is uniform. So in this situation, we find that for large t the map P_{p, q}(t) tends to a limit P^∞_{p, q}, which is uniquely defined by the conditions that P^∞_{p, q}|_{Tp, qΛ0} = I and P^∞_{p, q}(V^s_{p, q}) = V_{p, q}. Therefore, in this situation, the metaplectic correction tends for large t to a limit, too, and exponentially fast if the hyperbolicity is uniform.

In order to illustrate how the theory that we developed in the previous sections works, it is very instructive to look at a couple of simple examples. They will allow us to understand in more detail the interplay among the classical dynamics, the choice of an initial Lagrangian manifold and the analytical constructions we developed.

6.1. The free particle

The first example, we look at is the free particle in one dimension, with the Hamilton operator $-\frac{\hbar ^2}{2}\Delta$. We want to propagate a state that is initially localized at $(p,q)\in \mathbb{R}^2$, where p ⩾ 0. As an initial phase function, we choose

$$\begin{equation} S_0(x)=p(x-q)+\frac{\alpha }{2}(x-q)^2, \end{equation} \tag{ 90 }$$

where $\alpha \in \mathbb{R}$ is a real parameter. The corresponding Lagrangian manifold is

$$\begin{equation} \Lambda _0=\lbrace (p+\alpha x, q+x), x\in \mathbb{R}\rbrace , \end{equation} \tag{ 91 }$$

which is a line through (p, q) with slope α. The corresponding classical Hamilton function is $H(\xi )=\frac{1}{2}\xi ^2$ and the Hamiltonian flow it generates is given by Φ^t(ξ, x) = (ξ, x + tξ). Applying the flow to Λ₀ gives $\Lambda _t=\Phi ^t(\Lambda _0)=\lbrace (p+\alpha x, q+x +t(p+\alpha x)), x\in \mathbb{R}\rbrace ,$ and by replacing x with x/(1 + αt), this can be rewritten as

$$\begin{equation} \fl \Lambda _t=\lbrace (p+\alpha (t) x, q(t)+x) , x\in \mathbb{R}\rbrace , \qquad {\rm where}\quad \alpha (t)=\frac{\alpha }{1+\alpha t},\quad \ q(t)=q+tp, \end{equation} \tag{ 92 }$$

which is a line through (p, q(t)) = Φ^t(p, q) with the slope α(t). Note that if α < 0, i.e. if the initial line has a negative slope, then α(t) → ∞ for t → −1/α; this means that the line Λ_t turns vertical, and hence we have a caustic. But if α ⩾ 0, then α(t) ⩽ α for all t ⩾ 0, and no caustics occur. The phase function that generates Λ_t and satisfies the Hamilton–Jacobi equation with the initial condition S₀ is given by

$$\begin{equation} S(t,x)=\frac{1}{2} p^2t+p(x-q(t))+\frac{\alpha (t)}{2}(x-q(t))^2. \end{equation} \tag{ 93 }$$

To find the transport operator, we have to find the map $\phi _{\Lambda }(t,x)=\pi _{\Lambda _t}\Phi ^t\big(\pi ^{-1}_{\Lambda _0}(x)\big)$. To this end, we notice that $\pi ^{-1}_{\Lambda _0}(x)=(\nabla S_0(x),x)=(p+\alpha (x-q),x)$, and hence $\Phi ^t\big(\pi ^{-1}_{\Lambda _0}(x)\big)=(p+\alpha (x-q),x+t(p+\alpha (x-q)),$ and the final projection just takes the second component, therefore,

$$\begin{equation} \phi _{\Lambda }(t,x)=(1+\alpha t)x+t(p-\alpha q). \end{equation} \tag{ 94 }$$

Since ϕ'_Λ(t, x) = (1 + αt), the condition α ⩾ 0 guarantees that ϕ_Λ is non-contracting and we obtain $\phi _{\Lambda }^{-1}(t,x)=\frac{1}{1+\alpha t} (x-t(p-\alpha q))$. Therefore, the action of the transport operator reads

$$\begin{equation} T(t)A(x)=\frac{1}{(1+\alpha t)^{1/2}}A\bigg (\frac{1}{1+\alpha t} (x-t(p-\alpha q))\bigg ), \end{equation} \tag{ 95 }$$

which for α ⩾ 0 is well defined for all t ⩾ 0. In order to compute the metaplectic correction, which is generated in (49), we note that ${\mathcal A}(t,x)=(1+\alpha t)^{-2}$ is independent of x, and hence, we have

$$\begin{equation} \Delta _q(t)=\Delta (t)=\frac{1}{(1+\alpha t)^2}\Delta . \end{equation} \tag{ 96 }$$

Then, the metaplectic correction (55) is

$$\begin{equation} M_q(t)=\mathrm{e}^{\frac{\mathrm{i}\hbar }{2} \frac{t}{1+\alpha t}\Delta }, \end{equation} \tag{ 97 }$$

where the time-dependent factor in the exponent comes from the integration of the time dependence in Δ(t), $\int _0^t\frac{1}{(1+\alpha t^{\prime })^2}\, \mathrm{d}t^{\prime }=\frac{t}{1+\alpha t}$. We see that in particular, as we observed at the end of section 5, that the operator tends to a limit for large t, $M_q(t)\rightarrow M_q(\infty )=\mathrm{e}^{\frac{\mathrm{i}\hbar }{2} \frac{1}{\alpha }\Delta }$.

We computed all the elements in our extended time-dependent WKB propagation scheme for the free particle. Let us apply this to an initial state $\psi _0(x)=A_0(x)\mathrm{e}^{\frac{\mathrm{i}}{\hbar }[p(x-q)+\frac{\alpha }{2}(x-q)^2]}$ with A₀ = L_qa₀ for some $a_0\in S(\mathbb{R})$. Then, M_q(t)A₀ = L_qa_t with a_t = D_q(t)a₀ and since D_q(t) tends to a limit for large t, we have that $a_t\in S(\mathbb{R})$ with bounds uniform in $t\in \mathbb{R}^+$. Hence, we find

$$\begin{equation} \psi (t,x)=A(t,x)\, \mathrm{e}^{\frac{\mathrm{i}}{\hbar }[\frac{1}{2} p^2t+p(x-q(t))+\frac{\alpha (t)}{2}(x-q(t))^2]}, \end{equation} \tag{ 98 }$$

with

$$\begin{equation} \fl A(t,x)=T(t)M_q(t)A_0(x)=\frac{1}{\hbar ^{1/4}(1+\alpha t)^{1/2}}a_t\bigg (\frac{1}{\hbar ^{1/2}(1+\alpha t)} (x-t(p-\alpha q))\bigg ). \end{equation} \tag{ 99 }$$

Since the Hamiltonian is quadratic, this expression is actually exact, and the remainder terms are zero. Note that the behaviour of the amplitude A(t, x) depends on the size of ℏ^1/2(1 + αt), which determines the Ehrenfest time scale (8). If ℏ^1/2(1 + αt) → 0, the amplitude becomes concentrated, but if ℏ^1/2(1 + αt) ∼ 1, then the amplitude A(t, x) is a smooth function and ψ(t) actually becomes a Lagrangian state.

Let us note that for the special case that our initial state is a Gaussian coherent state,

$$\begin{equation} \psi _0(x)=\frac{1}{(\pi \hbar )^{\frac{1}{4}}}\, \mathrm{e}^{\frac{\mathrm{i}}{\hbar }[p(x-q)+\frac{b}{2}(x-q)^2]}, \end{equation} \tag{ 100 }$$

where $b\in \mathbb{C}$ with Im b > 0, our result is identical to the one obtained by the standard propagation of coherent states, since the Hamiltonian is quadratic. The only difference is that we derived it in a different way. But it was still instructive to go through the derivation since it highlights a few important points. We can rewrite the state using the phase function S₀ from (90) as

$$\begin{equation} \psi (x)=[L_qa_0](x)\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S_0(x)} ,\qquad {\rm with}\quad a_0(x)=\mathrm{e}^{\mathrm{i}\frac{b-\alpha }{2}\, x^2}, \end{equation} \tag{ 101 }$$

and then applying (98) and (99) gives us the propagated state. But since the initial state is independent of α, the final result is independent of the choice of α, too, i.e. of the initial Lagrangian submanifold used to propagate the state. But the intermediate steps depend on α, first of all, we needed to choose α ⩾ 0 in order to avoid caustics, and then the choice of α determined how we split the evolution into transport, implemented by the operator T(t), and dispersion, represented by D(t) = M₀(t). The simplest choice would have been α = 0; then, Λ₀ is horizontal, and actually invariant under the classical flow. In this case, T(t) is defined by transport along one central trajectory, and all the dispersion is described by D(t). If α > 0, then Λ₀ is transversal to the flow and the action of T(t) actually takes a family of trajectories into account. Since these trajectories move with different speed, this already accounts for part of the dispersion of the wave packet, and the operator M(t) has only to add the dispersion transversal to Λ_t. This is reflected by the fact that the coefficients of Δ(t) tend to 0 for large t, and fast enough so that the limit

$$\begin{equation} \lim _{t\rightarrow \infty } M_0(t)=\mathrm{e}^{\frac{\mathrm{i}\hbar }{2}\frac{1}{\alpha }\Delta } \end{equation} \tag{ 102 }$$

actually exists if α > 0. This means that for large t we can replace M₀(t) by its limit and still obtain a good approximation of the propagated state. This is why α > 0 is preferable to α = 0, in particular if we extend these constructions to more general systems.

In addition, the fact that for α > 0 the coefficients of Δ(t) tend to 0 for large t reflects the property that the map ϕ_Λ is expanding; this also has consequences for the size of the remainder term if we add higher order terms to the Hamiltonian. To study this effect will be the main focus of the next example.

6.2. Integrable systems

Let us look now at a more general Hamiltonian that is a function of the momentum only, but no longer necessarily quadratic, so H = H(ξ), and we will assume that for ξ > 0 we have H'(ξ) > 0 and H''(ξ) ⩾ 0. This is the form a Hamiltonian of an integrable system takes in action–angle coordinates.

The classical dynamics is given by Φ^t(ξ, x) = (ξ, x + tH'(ξ)), and by the condition H''(ξ) ⩾ 0, the velocity H'(ξ) cannot decrease with increasing ξ. This implies that if we take again an initial phase function of the form S₀(x) = px + αx²/2, then the corresponding Lagrangian submanifold does not develop caustics for t ⩾ 0 if α ⩾ 0, We find $\Lambda _t=\Phi ^t(\Lambda _0)=\lbrace (p+\alpha x, x+tH^{\prime }(p+\alpha x)),x\in \mathbb{R}\rbrace$ and therefore

$$\begin{equation} \fl \phi _{\Lambda }(t,x)=x+tH^{\prime }(p+\alpha x)\qquad {\rm and}\qquad \phi _{\Lambda }^{\prime }(t,x)=1+\alpha tH^{\prime \prime }(p+\alpha x), \end{equation} \tag{ 103 }$$

and so if α > 0 and H'' > 0, then ϕ_Λ(t, x) is actually expanding. Since the flow is now no longer linear, Λ_t is no longer a straight line and we do not attempt to find an explicit expression for S(t, x), but rather take it as defined by the Hamilton–Jacobi equation. But we observe that since Λ_t = {(∇S(t, ϕ_Λ(t, x)), ϕ_Λ(t, x))}, we have ∇S(t, ϕ_Λ(t, x)) = p + αx = ∇S₀(x), which is of course a consequence of momentum conservation.

Having the phase function and the map ϕ_Λ(t), and therefore, the transport operator T(t), we have the ingredients of the standard WKB propagation, and we now want to compute the metaplectic correction for an initial state localized at (p, 0). Here, we have to use the generalized formulation from section 5 since the Hamiltonian is no longer a sum of a quadratic kinetic energy term and a potential energy term. The transport operator is the quantization of the map Φ^t₁(ξ, x) = ([ϕ'_Λ(t, x)†]⁻¹ξ, ϕ_Λ(t, x)) and so the operator $\widehat{\delta H}$ has the symbol δH(t, ξ, x) = δH⁽⁰⁾(t, Φ^t₁(ξ, x)) + O_t(ℏ²), where δH⁽⁰⁾(t, ξ, x) = H(∇S(t, x) + ξ) − H(∇S(t, x)) − ξ · ∇H(∇S(t, x)) and the O_t(ℏ²) term comes from Egorov's theorem. Therefore, using ∇S(t, ϕ_Λ(t, x)) = ∇S₀(x), we find

$$\begin{eqnarray} &&\fl \delta H(t,\xi , x)=H(\nabla _0 S(x)+[\phi _{\Lambda }^{\prime }(t,x)^{\dagger }]^{-1}\xi )\nonumber\\ &&-\,H(\nabla S_0(x))-[\phi _{\Lambda }^{\prime }(t,x)^{\dagger }]^{-1}\xi \cdot \nabla H (\nabla S_0(x))+O_t(\hbar ^2). \end{eqnarray} \tag{ 104 }$$

Since the point (p, 0) corresponds to ξ = 0, x = 0, we find using (103) that the generator of the metaplectic correction is $\widehat{\delta H_2}(t)=-\frac{\hbar ^2}{2}\frac{H^{\prime \prime }(p)}{(1+\alpha t H^{\prime \prime }(p))^2} \Delta$ and hence

$$\begin{equation} M_0(t)=\mathrm{e}^{\frac{\mathrm{i}\hbar }{2}\frac{H^{\prime \prime }(p)t }{1+\alpha t H^{\prime \prime }(p)} \Delta }. \end{equation} \tag{ 105 }$$

We observe the same phenomenon as in the free case, and for α > 0, the large t behaviour has a limit. But now this is also reflected in the behaviour of the remainder term, which is determined by the size of

$$\begin{equation} \fl \delta H(t,\sqrt{\hbar }\xi , \sqrt{\hbar }x)-\delta H_2(t,\sqrt{\hbar }\xi , \sqrt{\hbar }x)= O(\hbar ^{3/2}[\phi _{\Lambda }^{\prime }(t,x)\phi _{\Lambda }^{\prime }(t,x)^{\dagger }]^{-1}), \end{equation} \tag{ 106 }$$

and since the time integral of this quantity is bounded, we obtain for a coherent state ψ₀ concentrated at (p, 0) that

$$\begin{equation} \Vert U(t)\psi _0-U_1(t) M_{0}(t)\psi _0\Vert \le C\sqrt{\hbar }. \end{equation} \tag{ 107 }$$

Here, we made the assumption that we can also control the remainder term O_t(ℏ²) from the application of Egorov's theorem, which seems likely since the transport operator has still a simple form. But we leave a detailed investigation to a future publication.

This is in contrast to the standard coherent state propagation that can only work if $\vert t\vert \sqrt{\hbar }\ll 1$, since $T_E=1/\sqrt{\hbar }$ for an integrable system. The physical reason for this is that for $t\sqrt{\hbar }\sim 1$ the state U₁(t)M(t)ψ₀ is no longer a coherent state but has become a Lagrangian state, i.e. there is a qualitative change in the nature of the state at the Ehrenfest time.

6.3. A parabolic barrier

Let us now look at a case with a hyperbolic trajectory, let $H(\xi ,x)=\frac{1}{2}\xi ^2-\frac{V_0}{2}x^2$; this Hamiltonian describes the motion in a parabolic barrier if V₀ > 0. The classical flow is given by

$$\begin{equation} \Phi ^t(\xi ,x)=(\cosh (\lambda t)\xi +\lambda \sinh (\lambda t)x,\lambda ^{-1}\sinh (\lambda t)\xi +\cosh (\lambda t)x), \end{equation} \tag{ 108 }$$

where $\lambda :=\sqrt{V_0}>0$ is the Liapunov exponent. The origin is a hyperbolic fixed point with stable and unstable manifolds given by

$$\begin{equation} V^{s}=\lbrace (-\lambda x,x) , x\in \mathbb{R}\rbrace ,\qquad V^{u}=\lbrace (\lambda x,x), x\in \mathbb{R}^n\rbrace . \end{equation} \tag{ 109 }$$

Let us first choose an initial state localized on the fixed point, i.e. at (0, 0). An initial phase function S₀(x) = α²x²/2 corresponds to $\Lambda _0=\lbrace (\alpha x,x)\, ,\, x\in \mathbb{R}\rbrace$ and the transported manifold Λ_t = Φ^t(Λ₀) can be written as

$$\begin{equation} \fl \Lambda _t=\lbrace (\alpha (t) x,x), x\in \mathbb{R}\rbrace ,\qquad {\rm with}\quad \alpha (t)=\frac{\cosh (\lambda t)\alpha +\lambda \sinh (\lambda t)}{\lambda ^{-1}\sinh (\lambda t)\alpha +\cosh (\lambda t)}. \end{equation} \tag{ 110 }$$

The corresponding phase function is S(t, x) = α(t)x²/2, and furthermore, we find

$$\begin{equation} \phi _{\Lambda }(t,x)=(\alpha \lambda ^{-1} \sinh (\lambda t)+\cosh (\lambda t))x, \end{equation} \tag{ 111 }$$

and therefore,

$$\begin{equation} \Delta (t)=(\alpha \lambda ^{-1} \sinh (\lambda t)+\cosh (\lambda t))^{-2}\Delta . \end{equation} \tag{ 112 }$$

The manifold Λ_t does not develop caustics if α ⩾ −λ. The case α = −λ is special since then Λ_t = V^s, i.e. Λ_t = Λ₀ is invariant and equal to the stable manifold. In this case, the manifold is contracting, but for all α > −λ the manifold is non-contracting. Since the case α = λ is also special, Λ_t = V^u for all t, and so the initial manifold is the unstable manifold; in all other cases, Λ_t → V^u for large t. These different cases are also reflected in the behaviour of T(t) and M(t). For α = −λ, we have ϕ_Λ(x) = e^−λtx, and hence,

$$\begin{equation} T(t)A(x)=\mathrm{e}^{\lambda t/2}A(\mathrm{e}^{\lambda t}x)\qquad {\rm and}\qquad {M}(t)=\mathrm{e}^{\frac{\mathrm{i}\hbar }{2}\frac{\mathrm{e}^{2\lambda t}-1}{2\lambda }\Delta }. \end{equation} \tag{ 113 }$$

The transport operator T(t) is squeezing at an exponential rate, and the metaplectic correction has to make up for this at an exponential rate, too. In contrast, if α > −λ, we have

$$\begin{equation} \alpha \lambda ^{-1} \sinh (\lambda t)+\cosh (\lambda t)=\frac{\alpha +\lambda }{2\lambda }\,\mathrm{e}^{\lambda t}+O(\mathrm{e}^{-\lambda t}) \end{equation} \tag{ 114 }$$

and so T(t)A(x) = [ϕ'_Λ(t)]^−1/2A(ϕ_Λ⁻¹(t, x)) is stretching at an exponential rate. The metaplectic correction is given by

$$\begin{equation} M(t)=\mathrm{e}^{\frac{\mathrm{i}\hbar }{2}[\frac{2\lambda }{(\alpha +\lambda )^2}+O(\mathrm{e}^{-\lambda t})]\Delta } \end{equation} \tag{ 115 }$$

and tends exponentially fast to a limit. This dichotomy between the cases α = −λ and α > −λ is analogous to the one between α = 0 and α > 0 that we found for the free particle. For α > −λ, the metaplectic correction saturates in time, whereas for α = −λ, its contribution grows exponentially. Since the Hamiltonian is quadratic, both cases give exact solutions, and the different initial manifolds only correspond to different ways to split the time evolution. But as in the integrable case, if we move to a perturbation, the contracting case α = −λ gives remainder terms that blow up, whereas the expanding cases α > −λ give remainder terms that remain bounded by $\sqrt{\hbar }$ independent of time.

So far we have looked at an initial state that is concentrated on top of the barrier; let us now look at a initial state localized at (p, q) with q < 0 and p > 0, hence an incoming state. Then the question of interest is if this state gets transmitted or reflected, and how one can describe the transition between reflection and transmission uniformly if one changes (p, q). We chose an initial phase function S₀(x) = p(x − q) + α(x − q)²/2, which again for α > −λ gives an expanding initial manifold Λ₀. For simplicity, we will choose α = λ, i.e. Λ₀ = (p, q) + V^u; the general case α > −λ can be treated similarly. From computing Φ^t(∇S₀(x), x), we find

$$\begin{equation} \phi _{\Lambda }(t,x)=\mathrm{e}^{\lambda t}(x-q)+q(t), \end{equation} \tag{ 116 }$$

where q(t) = λ⁻¹psinh (λt) + qcosh (λt), and therefore,

$$\begin{equation} T(t)L_qa(x)=\frac{\mathrm{e}^{-\lambda t/2}}{\hbar ^{1/4}}a\bigg (\frac{\mathrm{e}^{-\lambda t}}{\sqrt{\hbar }}[x-q(t)]\bigg ), \qquad M(t) =\mathrm{e}^{\frac{\mathrm{i}}{2}\frac{1-\mathrm{e}^{-2\lambda t}}{2\lambda } \Delta }. \end{equation} \tag{ 117 }$$

For the phase function, we find

$$\begin{eqnarray} &&\fl S(t,x)=\mathcal {L}(t)+p(t)(x-q(t))+\lambda (x-q(t))^2/2,\qquad {\rm where}\quad (p(t),q(t))=\Phi ^t(p,q),\nonumber\\ \end{eqnarray} \tag{ 118 }$$

and $\mathcal {L}(t)=\int _0^t (p\dot{q}-H(p,q))\, \mathrm{d}t$. We see that M(t) again tends with exponential speed to a limit $M^{\infty }= \mathrm{e}^{\frac{\mathrm{i}}{2}\frac{1}{2\lambda } \Delta }$. The behaviour of T(t)L_qa depends crucially on $\mathrm{e}^{-\lambda t}/\sqrt{\hbar }$; for small time, this parameter is large, and hence, the state is localized, but for

$$\begin{equation} \mathrm{e}^{-\lambda t}\approx \sqrt{\hbar }, \end{equation} \tag{ 119 }$$

i.e. on Ehrenfest time scales, the function T(t)L_qa is no longer strongly localized and the transported wave packet is actually a Lagrangian state associated with

$$\begin{equation} \Lambda _t=\Phi ^t(\Lambda _0)=(p(t),q(t))+V^u. \end{equation} \tag{ 120 }$$

Where this state moves depends on q(t), and for long times, we have

$$\begin{equation} q(t)=\frac{p+\lambda q}{2\lambda }\,\mathrm{e}^{\lambda t}+O(\mathrm{e}^{-\lambda t}). \end{equation} \tag{ 121 }$$

So the behaviour of q(t) depends on the value of p + λq. The case p + λq = 0 corresponds to initial conditions (p, q) in the stable manifold V^s of the fixed point at (0, 0), and so the trajectory runs into the fixed point and the Lagrangian state is for long times located on the unstable manifold V^u. The case p + λq > 0 gives a wave packet that is transmitted over the barrier and the case p + λq < 0 corresponds to a wave packet which is reflected. The strength of our approach is that we can describe the transition between these different cases in a uniform way given in (117). We studied a simple example here, but the method works for general potential barriers and allows us to describe the transition between reflection and transmission of time-dependent wave packets in a uniform way.

Such processes are of great importance in the theory of chemical reactions and we expect that our method could be of great use in the theoretical and quantitative description of time resolved chemical reactions, which are experimentally accessible in femto and atto chemistry [35, 36]. A first step will be to analyse the transport of wave packets over barriers, or more generally through bottlenecks in phase space, using the method of quantum normal forms recently developed in [37, 38].

In this section, we want to discuss the case that the initial coherent state is concentrated in a point z = (p, q) on a hyperbolic trajectory z(t) = Φ^t(z). As a test system, we choose the kicked harmonic oscillator

$$\begin{equation} H(t,p,q,t) = \frac{1}{2} ( p^2 + q^2 ) + K \cos q \sum _{n = -\infty }^{\infty } \delta ( t - n ), \end{equation} \tag{ 122 }$$

where K is the chaoticity parameter that we will choose to be K = 2. This system has also been used in [25]. In figure 1(a), we show a portrait of the classical phase space, the system has an unstable fixed point at the origin and we display the unstable manifold. The Liapunov exponent of the unstable fixed point is λ = 0.83... and we use ℏ = 0.0008, so the Ehrenfest time is

$$\begin{equation} T_E=4.29\ldots. \end{equation} \tag{ 123 }$$

Zoom In Zoom Out Reset image size

As an initial state, we choose a Gaussian wave packet centred on the fixed point at the origin:

$$\begin{equation} \psi _0(x)= \frac{1}{ (\pi \hbar )^{1/4} } \, {\rm e}^{ -x^2/2 \hbar }. \end{equation} \tag{ 124 }$$

In figure 1(b), we display the Wigner function of the time-evolved state at t = 4, which can be seen to be stretched along the unstable manifold. As expected at the Ehrenfest time, the state becomes a WKB-type state associated with the unstable manifold. For comparison, the real part of the time-evolved wavefunction ψ(t) is plotted in figure 2(c).

Zoom In Zoom Out Reset image size

The prediction of the theory in section 3 is that, after a suitable metaplectic correction of the initial amplitude, the state can be propagated using the standard time-dependent WKB approximation. This means that if we decompose the initial state as

$$\begin{equation} \psi _0 (x) = (L_0a)(x) \, \mathrm{e}^{i S_0 (x)/\hbar }, \end{equation} \tag{ 125 }$$

then the evolved state is from (19) and (34) given by

$$\begin{equation} \fl \psi (t,x) =[ T(t) D(t) L_0a_0 ](x)\, {\rm e}^{{\rm i} S(t,x)/\hbar } \approx [ T(t) M_0(t) L_0a] (x) \, {\rm e}^{{\rm i} S(t,x)/\hbar } . \end{equation} \tag{ 126 }$$

Here, S(t, x) is the solution of the Hamilton–Jacobi equation satisfying S(0, x) = S₀(x). The operators T(t), D(t) and M₀(t) also depend on the choice of S₀(q). Even if according to section 4.4 this choice is, to some extent, arbitrary, the efficiency of the method may be sensitive to S₀(x). We shall only consider the quadratic case $S_0(x)=\frac{\tan \theta }{2}\, x^2$, corresponding to a linear Lagrangian manifold p = tan (θ) q. We shall see that there is a wide range of initial manifolds that work well.

The metaplectic correction M₀(t) is calculated by considering the linearized dynamics at the fixed point, i.e. we approximate cos (q) by 1 − q² in (122). Within this approximation, the initial manifold remains a straight line through the origin for all times and the transport operator T(t) (14) becomes a scaling transformation:

$$\begin{equation} \left[ T(t) a \right](q)= \sqrt{\kappa } \, a( \kappa q ) , \end{equation} \tag{ 127 }$$

with κ(t) being a time-dependent factor (which can be calculated by following the linear evolution of the initial manifold). Accordingly, the generalized Laplacian (17) reads

$$\begin{equation} \Delta (t)= \kappa ^2(t) \, \Delta, \end{equation} \tag{ 128 }$$

and the metaplectic correction is just

$$\begin{equation} M_0(t)= \exp \left( \frac{{\rm i} \hbar }{2} \int _0^t \kappa ^2(t^{\prime })\, {\rm d}t^{\prime } \, \Delta \right). \end{equation} \tag{ 129 }$$

Note that equation (126) is equivalent to saying that the sequence of operations

$$\begin{equation} T^\ast (t) [ \psi (t) \, {\rm e}^{-{\rm i} S(t)/\hbar } ] =D(t)L_0\, a_0\approx M_0(t)L_0 \, a_0 \, \end{equation} \tag{ 130 }$$

must produce approximately a Gaussian state centred at 0, for a Gaussian initial amplitude a₀(x). This is the first test we shall carry out, with the choice S₀ = 0 corresponding to the initial Lagrangian manifold p = 0. First, we propagate the state (125) exactly during some time, and then propagate it back using time-dependent WKB. Figures 2(a) and (b) show the comparison of the exact result D(t) L₀a₀ with the metaplectic approximation M₀(t) L₀a₀. In the case of the modulus, the agreement between the exact DL₀a₀ and M₀L₀a₀ is perfect within visual resolution. There is a small difference between the phases, but this occurs in a region where the amplitude is very small. We also show a direct comparison between the exact and the metaplectically extended WKB propagated states in figures 2(c) and (d). The agreement is excellent.

We will now investigate how robust the method is with respect to changing the choice of the (linear) initial WKB manifold. We have changed the slope θ in a range of almost 90°, from θ = −0.30π/2 to θ = 0.65π/2. The comparison between exact and metaplectically extended WKB at t = 4 is shown in figure 3. Except for very slight differences (blue triangles), the agreement is still excellent for these 'large' slopes.

Zoom In Zoom Out Reset image size

As we increase the slope, the amplitude concentrates in narrower regions. This happens because, as the slope is large, the map induced in the q-coordinate by the WKB manifold is more expansive and T* is more contracting.

We tested the theory also for different values for ℏ and found good agreement with the expected behaviour (not shown).

In this paper, we derived an extension of the standard time-dependent WKB method that can be applied to highly localized states like coherent states. It allows us to describe in a uniform way the transition in time from a semiclassically highly localized coherent state to a delocalized Lagrangian state, which takes place at the Ehrenfest time.

The main idea on which this extension of time-dependent WKB theory is built is an exact decomposition of the time evolution of an initial state of the form $\psi _0(x)=A_0(x)\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S_0(x)}$, where S₀ is real valued, into several parts

$$\begin{equation} \psi (t)=(T(t)D(t)A_0)\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S(t)}. \end{equation} \tag{ 131 }$$

Here, S(t) is a solution of the Hamilton–Jacobi equation and is hence related to the transport of the Lagrangian manifold Λ₀ generated by S₀ through phase space. The unitary operator T(t) transports functions in position space along the projections of the phase space trajectories emanating from Λ₀. Finally, D(t) is the propagator generated by the time-dependent Hamiltonian −ℏ²T*(t)ΔT(t)/2. S(t) and T(t) are defined purely in terms of transport along classical trajectories and D(t) takes into account the dispersive effect of quantum mechanics.

The standard time-dependent WKB approximation is obtained by approximating D(t) ≈ I; this works fine if the amplitude A₀ is sufficient flat, i.e. has bounded derivatives. But for a coherent state, A₀ is strongly localized around a point x = q; then, it is more natural to freeze the coefficients of the generator of D(t) at x = q, and the resulting operator is a metaplectic operator M_q(t) whose action on functions can be computed quite easily. The main result of this paper is thus

$$\begin{equation} \psi (t)\approx (T(t)M_q(t)A_0)\,\mathrm{e}^{\frac{\mathrm{i}}{\hbar }S(t)}, \end{equation} \tag{ 132 }$$

which is valid for amplitudes A₀ that are strongly localized around x = q. In order to justify the validity of this approximation for long times, particularly on Ehrenfest time scales, we had to analyse the underlying classical dynamics more carefully. We introduced a non-contraction condition on the position space trajectories emanating from a neighbourhood of the initial state with momenta p = ∇S(x), and if this condition holds our propagation scheme is effective. The non-contraction condition excludes caustics, but we have a large freedom in the choice of the initial phase function S(x), which in many cases allows us to avoid caustics, at least until the state becomes delocalized.

For times shorter than the Ehrenfest time, our scheme reproduces the standard coherent state propagation results based on a Taylor expansion of the Hamiltonian around the centre trajectory. But for times of the order of the Ehrenfest time, we find that the state becomes extended and a Lagrangian state, which in the chaotic case, is supported by the unstable manifold of the centre trajectory. From that time onwards, the standard time-dependent WKB theory applies, as has been observed in [25]. In particular, if the system is hyperbolic and mixing, one can apply the results from [29] to conclude that the state becomes equidistributed after the Ehrenfest time.

In order to extend the results to more general Hamiltonians, and to be able to include caustics, we noted that the standard time-dependent WKB approximation can be viewed as the exact quantum time evolution generated by a quantization of a first-order Taylor approximation of the Hamilton function around the Lagrangian manifold associated with the time-evolved WKB state. Based on this insight, we could now choose different first-order approximations, which remained valid for general Hamiltonians and at caustics. The price one has to pay is a more complicated and less explicit formalism. In addition, all the previous results were in principle rigorous, although we refrained from stating them in the form of theorems, but here we have to rely on an assumption that certain special cases of Egorov's theorem remain valid on Ehrenfest time scales. But we get a further benefit from this more general way to look at the time-dependent WKB approximation, we can show that the classical map associated with the metaplectic correction M_q(t) is a shear map on phase space and that it furthermore in many cases converges to a limit for large t, which can be determined from simple geometric considerations. This implies that in many cases

$$\begin{equation} \lim _{t\rightarrow \infty } M_q(t)=M_q^{(\infty )} \end{equation} \tag{ 133 }$$

and we can determine M^(∞)_q up to a phase from simple geometric considerations. If the centre trajectory is hyperbolic, then the limit in (133) is reached exponentially fast.

We illustrated and tested the theory with several examples. For the free particle, and more generally, for integrable systems, the dynamics can be computed quite explicitly. The Ehrenfest time is of order T_E ∼ ℏ^−1/2 and we see a qualitative transition in the nature of the state from a localized coherent state to an extended WKB state at this time scale. Similarly, for a parabolic barrier, the dynamics can be computed explicitly and the formalism developed in this work gives explicit formulae, which describe in a uniform way the transition from reflection to transmission of a wave packet when one varies the energy near the critical energy. We then carried out detailed numerical tests on the kicked harmonic oscillator for an initial state localized on a hyperbolic fixed point. These showed impressive agreement between the metaplectic extension of the WKB method and the exact quantum propagation. These tests also illustrated the simplicity of the method; the metaplectic correction is in fact a very simple operator, and if one has implemented the time-dependent WKB method, then adding the metaplectic correction is easy and at once allows one to propagate a much larger class of states.

There are many areas where the results from this paper should be of interest. In many physical systems, the Ehrenfest time for a realistic set of parameters is actually quite short, e.g., in quantum billiards coherent states spread out after a few bounces [22, 23], and the metaplectic extension of WKB should be able to describe this transition. The scattering of a wave packet of a barrier near the critical energy is another example; the Ehrenfest time is the time when the wave packet reaches the barrier, i.e. the time where the physically interesting processes start to occur. Metaplectically extended WKB allows us to describe this process explicitly and uniformly in ℏ and t. A large class of chemical reactions is described in the framework of transition state theory by the crossing of a barrier on a high-dimensional energy surface and modern experimental methods in atto and femto chemistry allow us to study the dynamics of such chemical reactions with impressive precision, [35, 36]. The metaplectic extension to WKB together with the normal form approach to transition state theory [37, 38] should allow us to give an efficient theoretical description of chemical processes on such time scales.

In this appendix, we collect some material on Weyl quantization, which we use in the main part. Much of it is a development of standard material and only the result on Egorov's theorem seems to be new. For the material we present, see [30, 31].

A.1. Weyl quantization

Let $\hat{W}(\xi , x):=\mathrm{e}^{\frac{\mathrm{i}}{\hbar } (\xi \hat{q}+x\hat{p})}$ be the Weyl operators that represent translations in phase space; then, we can define for a function A(p, q) on phase space an operator, its Weyl quantization, by

$$\begin{equation} \hat{A}=\int\!\!\!\int \mathcal {F}A(\xi ,x)\hat{W}(\xi ,x)\frac{\mathrm{d}x\,\mathrm{d}\xi }{(2\pi \hbar )^{2n}}, \end{equation} \tag{ A.1 }$$

where $\mathcal {F}A(\xi ,x)=\int\!\!\!\int A(p,q)\,\mathrm{e}^{-\frac{\mathrm{i}}{\hbar }[p\xi +qx]}\, \mathrm{d}p\,\mathrm{d}q$ denotes the Fourier transform of A. The function A(p, q) is called the Weyl symbol of $\hat{A}$, and the properties of $\hat{A}$ can often be determined from properties of A, e.g., if $A\in S^{\prime }(\mathbb{R}^n\times \mathbb{R}^n)$, then $\hat{A} :S(\mathbb{R}^n)\rightarrow S^{\prime }(\mathbb{R}^n)$. If $\hat{A}=|\psi \rangle \langle \psi |$ is the projection onto a state ψ, then A(ξ, x) is proportional to the Wigner function of the state ψ.

The product of operators can be expressed in terms of the symbols, and one has $\hat{A}\hat{B}=\widehat{A\sharp B}$, where

$$\begin{eqnarray} &&\fl A\sharp B(p,q)=A(p,q)\,\mathrm{e}^{\frac{\mathrm{i}\hbar }{2}[\overleftarrow{\nabla }_p\cdot \overrightarrow{\nabla }_q-\overleftarrow{\nabla }_q\cdot \overrightarrow{\nabla }_p]}B(p,q)\nonumber\\ &&\hphantom{\fl A\sharp B(p,q)}=A(p,q)B(p,q)+\frac{\mathrm{i}\hbar }{2} \lbrace A,B\rbrace (p,q)\nonumber\\ &&\hphantom{\fl A\sharp B(p,q)=}-\frac{\hbar ^2}{8}A(p,q)[\overleftarrow{\nabla }_p\cdot \overrightarrow{\nabla }_q-\overleftarrow{\nabla }_q\cdot \overrightarrow{\nabla }_p]^2B(p,q)+\cdots . \end{eqnarray} \tag{ A.2 }$$

Here, {A, B}(p, q) denote the Poisson bracket, and the arrows over the derivatives indicate if they act on the function on the left or on the right. Suitable conditions on the functions A and B under which this expansion holds can be found in [30].

A.2. Dynamics generated by first-order operators

Let K(t, x, ξ) = X(t, x) · ξ, where X(t, x) is a time-dependent vector field. We are interested in the time evolution T(t, s) generated by

$$\begin{equation} \hat{K}(t)=-\mathrm{i}\hbar X(t,x)\cdot \nabla -\frac{\mathrm{i}\hbar }{2}\nabla X(t,x). \end{equation} \tag{ A.3 }$$

Let ϕ(t, s, x) be the family of maps generated by the vector field X(t, x):

$$\begin{equation} \partial _t\phi (t,s;x)=X(t,\phi (t,s;x)),\qquad \phi (t,t;x)=I. \end{equation} \tag{ A.4 }$$

Then,

$$\begin{equation} (T(t,s)A)(x)=[\det \phi ^{\prime }(t,s;\phi ^{-1}(t,s;x))]^{-1/2}A(\phi ^{-1}(t,s;x)). \end{equation} \tag{ A.5 }$$

A.3. Egorov's theorem

Egorov's theorem is one way to formulate the correspondence principle; it gives a general relation between classical and quantum dynamics. Let H(t) be a real-valued smooth phase space function, $\hat{H}(t)$ be its Weyl quantization, and Φ^t and U(t) be the classical and quantum time evolution generated by H(t) and $\hat{H}(t)$, respectively. Then, Egorov's theorem states that

$$\begin{equation} U^*(t)\hat{A}U(t)=\hat{A}_t +O_t(\hbar ^2), \end{equation} \tag{ A.6 }$$

where A_t = A○Φ^t, and A and H have to satisfy some conditions on their smoothness and growth at infinity, see [28]. Since A_t is the classical time evolution of A, this means that quantum and classical evolution are close for small ℏ. The remainder term O_t(ℏ²) does depend on time, and the best general estimates are of the form O_t(ℏ²) ≪ ℏ² e^Γt for some constant Γ > 0, which depends on H. We want to discuss two cases in which one has better control over the remainder.

One approach to Egorov's theorem is based on writing Heisenberg's equation of motion for $\hat{A}(t)=U^*(t)\hat{A}U(t)$, i.e. $\mathrm{i}\hbar \partial _t\hat{A}(t)=[\hat{H},\hat{A}]$, in terms of the symbols using (A.2), which gives

$$\begin{eqnarray} &&\fl \partial _t A(p,q)=\frac{2}{\hbar }A(p,q) \sin (\hbar [\overleftarrow{\nabla }_p\cdot \overrightarrow{\nabla }_q-\overleftarrow{\nabla }_q\cdot \overrightarrow{\nabla }_p]/2)H(p,q)\nonumber\\ &&\hphantom{\fl \partial _t A(p,q)}=\lbrace A,H\rbrace (p,q)+\frac{\hbar ^2}{24}A(p,q)\big[\overleftarrow{\nabla }_p\cdot \overrightarrow{\nabla }_q-\overleftarrow{\nabla }_q\cdot \overrightarrow{\nabla }_p\big]^3 H(p,q)+\cdots . \end{eqnarray} \tag{ A.7 }$$

The leading order term is just the Liouville equation and gives A ≈ A○Φ^t, the main problem is then to control the higher order terms.

A.3.1. Quadratic Hamiltonians and metaplectic operators

If H(p, q) is a quadratic function of p and q, with possibly time-dependent coefficients, then the ℏ expansion in (A.7) terminates after the leading term, and hence, the evolution equation for A is just the classical Liouville equation. So in this case

$$\begin{equation} U^*(t)\hat{A}U(t)=\hat{A}_t, \end{equation} \tag{ A.8 }$$

where A_t = A○Φ^t, and one says Egorov's theorem is exact. The corresponding classical maps Φ^t are linear; hence, the operators U(t) for all quadratic H form a quantization of the symplectic group, which turns out to be a double cover of the symplectic group called the metaplectic group. The operators U(t) are often referred to as metaplectic operators [3, 32, 33].

A.3.2. Conjugation by a flow

The second case we need is that U(t) = T(t, 0), i.e. H(t) is linear in p and given by H(t, p, q) = K(t, p, q) = X(t, q) · p. The classical map Φ^t₁(p, q) generated by H is given by the solutions to

$$\begin{equation} \dot{\xi }(t) =-\nabla _q K(t,\xi (t),x(t)),\qquad \dot{x} (t)=\nabla _p K(t,\xi (t),x(t)) \end{equation} \tag{ A.9 }$$

with the initial conditions ξ(0) = p and x(0) = q. We can express Φ^t₁ in terms of the map ϕ(t, 0, x) as

$$\begin{equation} \Phi ^t_1(\xi ,x)=([\phi ^{\prime }(0,t,x)]^{\dagger }\xi ,\phi (t,0,x)), \end{equation} \tag{ A.10 }$$

where [ϕ'(0, t, x)]† is the transpose of the inverse of the matrix ϕ'(t, 0, x).

Let us consider first the case that A is linear in p, i.e. A = b(q) · p for some vector-valued function b(q); then, only the leading order term in (A.7) is non-zero, and hence

$$\begin{equation} \fl T^*(t,0)\hat{A} T(t,0)=\hat{A}_t,\qquad {\rm with}\quad A(t,q,p)=b(\phi (t,0;q))\cdot [\phi ^{\prime }(0,t,x)]^{\dagger }p, \end{equation} \tag{ A.11 }$$

without any remainder terms; hence, Egorov's theorem is exact again. Using this result, we can discuss the case we will need, namely the case that $\hat{A}$ is a second-order differential operator of the form

$$\begin{equation} \hat{A}=({\mathcal B}(q)\nabla ) \cdot {\mathcal B}(q)\nabla \end{equation} \tag{ A.12 }$$

so that $T^*\hat{A}T=T^*{\mathcal B}(q)\nabla T\cdot T^* {\mathcal B}(q)\nabla T,$ and if we denote the rows of ${\mathcal B}$ by b_j(q), we can use the previous result (A.11) and (A.2) to obtain $T^*(t,0)\hat{A} T(t,0)=\hat{A}(t)$ with

$$\begin{eqnarray} &&\fl A(t)=\sum _{i}(b_i(\phi (t,0;q))\cdot [\phi ^{\prime }(0,t,x)]^{\dagger }p) \sharp (b_i(\phi (t,0;q))\cdot [\phi ^{\prime }(0,t,x)]^{\dagger }p)\nonumber\\ &&\hphantom{\fl A(t)}=p \cdot \phi ^{\prime }(0,t,x)({\mathcal B}^{\dagger }{\mathcal B})(\phi (t,0;q))[\phi ^{\prime }(0,t,x)]^{\dagger }p\nonumber\\ &&\hphantom{\fl A(t)=}+\,\frac{\hbar ^2}{8}\sum _{i} (b_i\cdot [\phi ^{\prime }(0,t,q)]^{\dagger }p)(\overleftarrow{\nabla }_p\cdot \overrightarrow{\nabla }_q)(\overleftarrow{\nabla }_q\cdot \overrightarrow{\nabla }_p) (b_i(\phi (t,0;q))\cdot [\phi ^{\prime }(0,t,x)]^{\dagger }p),\nonumber\\ \end{eqnarray} \tag{ A.13 }$$

where we have used (A.2) . This can be further simplified, but we will restrict ourselves to the case ${\mathcal B}=I$; hence, b_i = e_i and then we find

$$\begin{equation} A(t,p,q)= p\cdot \phi ^{\prime }(0,t,x)[\phi ^{\prime }(0,t,x)]^{\dagger }p +\frac{\hbar ^2}{8}\sum _i \mathop {\mathrm{Tr}}\nolimits [\phi _i^{\prime \prime }(0,t,q)]^2, \end{equation} \tag{ A.14 }$$

which is the symbol of T*(t, 0)ΔT(t, 0). The term of order ℏ² contains second derivatives of the inverse of ϕ(t, 0; q) and so if ϕ(t, 0; q) is non-contracting, then these derivatives stay bounded.

The results in this appendix are used at the end of section 5 to show that the classical map associated with the metaplectic correction is uniquely determined by the dynamics of the tangent space to the initial Lagrangian submanifold and the vertical subspace.

Lemma B.1. Let $L_1,L_2\subset \mathbb{R}^n\times \mathbb{R}^n$ be Lagrangian subspaces with L₁∩L₂ = {0} and let L be another Lagrangian subspace with L∩L₁ = {0}. Then, there exists a unique linear symplectic map T with $T|_{L_1}=I$ and T(L₂) = L.

To show this, we will use mainly the following two standard facts from linear symplectic geometry; see, e.g., [31].

• (a)  
  By the Darboux theorem, there exist symplectic coordinates $(v,w)\in \mathbb{R}^n\times \mathbb{R}^n$ such that L₁ = {w = 0} and L₂ = {v = 0}.
• (b)  
  If L₃ is a Lagrangian subspace with L₃∩L₁ = {0}, then there exists a unique symmetric matrix A, such that L₃ = {v = Aw}.

Let us now prove the lemma. The condition $T|_{L_1}=I$ implies that in the coordinates from (a) the matrix representing T is of the form $\mathcal {M}_T=\big(\begin{array}{@{}cc@{}}{\ssty I} & {\ssty A}\\ [-5pt] {\ssty 0} &{\ssty B}\end{array}\big)$, where A and B are n × n matrices. Now this matrix must also be symplectic, i.e. M^t_TΩM_T = Ω with $\Omega =\big(\begin{array}{@{}cc@{}}{\ssty 0} & {\ssty -I}\\ [-5pt] {\ssty I} & {\ssty 0}\end{array}\big)$ and this gives B = I and A = A^t. Therefore,

$$\begin{equation} T(L_2)=\lbrace v=Aw\rbrace , \end{equation} \tag{ B.1 }$$

but from (b) the condition T(L₂) = L then uniquely determines A. So $M_T=\big(\begin{array}{@{}cc@{}} {\ssty I} & {\ssty A}\\ [-4pt] {\ssty 0} &{\ssty I}\end{array}\big)$ is the unique shear relative to L₁, which maps L₂ to L.