Frozen Gaussian approximation for high frequency wave propagation in periodic media

Abstract

Propagation of high-frequency wave in periodic media is a challenging problem due to the existence of multiscale characterized by short wavelength, small lattice constant and large physical domain size. Conventional computational methods lead to extremely expensive costs, especially in high dimensions. In this paper, based on Bloch decomposition and asymptotic analysis in the phase space, we derive the frozen Gaussian approximation for high-frequency wave propagation in periodic media and establish its converge to the true solution. The formulation leads to efficient numerical algorithms, which are presented in a companion paper [SIAM J. Sci. Comput. 38 (2016), A2440–A2463].

Keywords

Frozen Gaussian approximation high-freqency wave propagation lattice potential Bloch decomposition

1. Introduction

We are interested in studying high-frequency wave propagation in periodic media. A typical example is given by the following Schrödinger equation in the semiclassical regime with a superposition of a (highly oscillatory) microscopic periodic potential and a macroscopic smooth potential, $\begin{matrix} (1.1) & i ε \frac{\partial ψ^{ε}}{\partial t} = - \frac{ε^{2}}{2} Δ ψ^{ε} + V (x / ε) ψ^{ε} + U (x) ψ^{ε}, x \in R^{d}, \end{matrix}$ where V and U are smooth potential functions, V is periodic with respect to the lattice $Z^{d}$ : $V (x + e_{i}) = V (x)$ for any $x \in R^{d}$ and ${e_{i}, i = 1, 2, \dots, d}$ is the standard basis of $R^{d}$ . Here $ε ≪ 1$ is the rescaled Planck constant, $ψ^{ε}$ is the wave function, and d is the spatial dimensionality.

The equation (1.1) can be viewed as a model for electron dynamics in a crystal, where V is the effective periodic potential induced by the crystal, and U is some external macroscopic potential. Notice that we have identified the period of $V (x / ε)$ and the “semiclassical parameter” in front of the derivative terms. This parameter choice gives the most interesting case as $ε \to 0$ [2].

The mathematical analysis of this work is motivated by the challenge of numerical simulation of (1.1) when ε is small. In this semiclassical regime, the wave function $ψ^{ε}$ becomes oscillatory with wave length $O (ε)$ . This means a computational domain of order 1 size contains $O (1 / ε)$ wavelengths, and each of them needs to be resolved if conventional numerical methods are applied. For example, even for the simplest case $V = 0$ (no lattice potential), a mesh size of $O (ε)$ is required when using the time-splitting spectral method [1] to compute (1.1) directly; an even worse mesh size of $o (ε)$ is needed if one uses the Crank–Nicolson schemes [27] or the Dufort–Frankel scheme [28]. Besides, the presence of non-zero lattice potential introduces further difficulties which restrict the mesh size to be $o (ε)$ in the standard time-splitting spectral method [1]. Special techniques using Bloch decomposition are needed to relax the mesh size to be of $O (ε)$ [13–15]. Moreover, in these methods, a large domain is demanded in order to avoid the boundary effects. Therefore the total number of grid points is huge, which usually leads to unaffordable computational cost, especially in high ( $d > 1$ ) dimensions.

An alternative efficient approach is to solve (1.1) asymptotically by the Bloch decomposition and modified WKB methods [3,5,8], which lead to eikonal and transport equations in the semi-classical regime. An advantage of this method is that the computational cost is independent of ε. However, the eikonal equation can develop singularities which make the method break down at caustics. The Gaussian beam method (GBM) [32] was then introduced by Popov to overcome this drawback at caustics. The idea is to allow the phase function to be complex and choose the imaginary part properly so that the solution has a Gaussian profile; see [16–21,29,40,41] for recent developments. Similar ideas can be also found in the Hagedorn wave packet method [9,10]. Unlike the geometric optics based method, the Gaussian beam method allows for accurate computation of wave function around caustics [7,34]. But the problem is that the constructed beam must stay near the geometric rays to maintain accuracy. This becomes a drawback when the solution spreads [24,29,33].

The Herman–Kluk propagator [12,22,23] was proposed for Schrödinger equation without the oscillatory periodic background potential. The method was rigorously analyzed in [36,39] and further extended as the frozen Gaussian approximation (FGA) for general high frequency wave propagation in [24–26]. The FGA method uses Gaussian functions with fixed widths, instead of using those that might spread over time, to approximate the wave solution. Despite its superficial similarity with the Gaussian beam method, it is different at a fundamental level. FGA is based on phase plane analysis, while GBM is based on the asymptotic solution to a wave equation with Gaussian initial data. In FGA, the solution to the wave equation is approximated by a superposition of Gaussian functions living in phase space, and each function is not necessarily an asymptotic solution, while GBM uses Gaussian functions (called beams) in physical space, with each individual beam being an asymptotic solution to the wave equation. The main advantage of FGA over GBM is that the problem of beam spreading no longer exists.

In this paper, we extend FGA for computation of high-frequency wave propagation in periodic media. We mainly focus on the derivation of an integral representation formula of FGA in the phase space and establish the rigorous convergence results for FGA. While the FGA works for general strictly hyperbolic equations, we focus in this paper the case of semiclassical Schrödinger equation with periodic media (1.1). The computational algorithm and numerical results will be presented in a separate paper [6]. The rest of the paper is organized as follows. We first recall the Bloch decomposition of periodic media and introduce the windowed Bloch transform in Section 2. In Section 3, we present the formulation of FGA for periodic media and the main convergence result. The proof of the main result is given in Section 4.

Notations.
The absolute value, Euclidean distance, vector norm, induced matrix norm, and sum of components of a multi-index will all be denoted by $| \cdot |$ . We will use the standard notations $S$ , $C^{\infty}$ , and $C_{c}^{\infty}$ for Schwartz class functions, smooth functions, and compactly supported smooth functions, respectively. We will sometimes use subscripts to specify the dependence of a constant on the parameters, for instance, notations like $C_{T}$ to specify the dependence of a constant on a parameter T.

2. Bloch decomposition and windowed Bloch transform

The frozen Gaussian approximation for periodic media relies crucially on the Bloch decomposition to capture the fine scale ( $O (ε)$ spatial scale) oscillation. First we briefly recall the well-known Bloch–Floquet decomposition for Schrödinger operators with a periodic potential.

Consider a Schrödinger operator $\begin{matrix} (2.1) & H = - \frac{1}{2} Δ + V (x), \end{matrix}$ where the potential V is periodic with respect to the lattice $Z^{d}$ . We denote Γ the unit cell of the lattice: ${Γ = [0, 1)}^{d}$ . The unit cell of the reciprocal lattice (known as the first Brillouin zone) is given by ${Γ^{*} = [- π, π)}^{d}$ . It is standard (e.g., [35]) that the spectrum of H is given by energy bands $\begin{matrix} spec (H) = ⋃_{n = 1}^{\infty} ⋃_{ξ \in Γ^{*}} E_{n} (ξ), \end{matrix}$ where for each $ξ \in Γ^{*}$ , ${E_{n} (ξ)}$ are the collection of eigenvalues (in ascending order) of the operator $\begin{matrix} H_{ξ} = \frac{1}{2} {(- i \nabla_{x} + ξ)}^{2} + V (x) \end{matrix}$ with periodic boundary condition on Γ. The Bloch waves are the associated eigenfunctions: For each band n and $ξ \in Γ^{*}$ , it solves $\begin{matrix} (2.2) & H_{ξ} u_{n} (ξ, \cdot) = E_{n} (ξ) u_{n} (ξ, \cdot) . \end{matrix}$ with periodic boundary condition on Γ, where $ξ$ serves as a parameter in the above equation. $u_{n} (ξ, \cdot)$ is normalized that $\begin{matrix} (2.3) & \int_{Γ} {| u_{n} (ξ, x) |}^{2} d x = 1 . \end{matrix}$ We extend $u_{n} (ξ, x)$ periodically with respect to the second variable, so it is defined on $Γ^{*} \times R^{d}$ . We will also write $u_{n, ξ} = u_{n} (ξ, \cdot)$ when the former is more convenient.

These Bloch waves generalize the Fourier modes (complex exponentials) to periodic media (see for example discussions in [8]). In particular, for any function $f \in L^{2} (R^{d})$ , we have the Bloch decomposition $\begin{matrix} (2.4) & f (x) = \frac{1}{{(2 π)}^{d / 2}} \sum_{n = 1}^{\infty} \int_{Γ^{*}} u_{n} (ξ, x) e^{i ξ \cdot x} {(B f)}_{n} (ξ) d ξ, \end{matrix}$ where the Bloch transform $B : L^{2} (R^{d}) \to L^{2} {(Γ^{*})}^{N}$ is given by $\begin{matrix} (2.5) & {(B f)}_{n} (ξ) = \frac{1}{{(2 π)}^{d / 2}} \int_{R^{d}} {\bar{u}}_{n} (ξ, y) e^{- i ξ \cdot y} f (y) d y . \end{matrix}$ As an analog of the Parseval’s identity, we have $\begin{matrix} (2.6) & \int_{R^{d}} {| f (x) |}^{2} d x = \sum_{n = 1}^{\infty} \int_{Γ^{*}} {| {(B f)}_{n} (ξ) |}^{2} d ξ . \end{matrix}$ As suggested by (2.4) and (2.5), we introduce the notation Ω to denote the phase space corresponding to one band ( $Γ^{*}$ is viewed as a torus, i.e., periodic boundary condition is assumed on $Γ^{*}$ ) $\begin{matrix} (2.7) & Ω : = R^{d} \times Γ^{*} = {(x, ξ) ∣ x \in R^{d}, ξ \in Γ^{*}} . \end{matrix}$ Correspondingly, we will use the notation $(q, p)$ for a point in Ω.

For later usage, we define the Berry phase $A_{n}$ for the Bloch waves, $\begin{matrix} (2.8) & A_{n} (ξ) = {⟨ u_{n} (ξ, \cdot), i \nabla_{ξ} u_{n} (ξ, \cdot) ⟩}_{L^{2} (Γ)} . \end{matrix}$ The normalization condition (2.3) implies $A_{n} (ξ)$ is always a real number. Berry phase is an important concept appearing during the adiabatic limit of quantum dynamics, as some external parameters are varied slowly [4,37]. As B. Simon observed in [38], the adiabatic Berry phase has an elegant mathematical interpretation as the holonomy of a certain connection, the Berry connection, in the appropriate fiber bundle. It has been used in the Bloch dynamics to explain various important phenomena in crystals, see for example [42] and references therein.

We should be cautious about one subtlety though as the eigenvalue equation (2.2) and the normalization only define $u_{n} (ξ, \cdot)$ up to a unit complex number, in particular, for any function φ periodic in $Γ^{*}$ , $\begin{matrix} (2.9) & v_{n} (ξ, x) = e^{i φ (ξ)} u_{n} (ξ, x), (x, ξ) \in Ω, \end{matrix}$ also provides a set of Bloch waves. This is known as the gauge choice for the Bloch waves. However, different gauge choice gives different values of $A_{n} (ξ)$ and even causes trouble if φ is discontinuous. While for the analysis, it suffices to assume smooth dependence of $u_{n}$ on $ξ$ (which is possible as the nth band is separated from the rest of the spectrum), this gauge freedom makes numerical computation nontrivial. We will further address this by designing a gauge-invariant algorithm in a companion paper [6] on the numerical algorithms.

Differentiating (2.2) with respect to $ξ$ produces $\begin{matrix} (2.10) & H_{ξ} \nabla_{ξ} u_{n} (ξ, x) + (- i \nabla_{x} + ξ) u_{n} (ξ, x) = E_{n} (ξ) \nabla_{ξ} u_{n} (ξ, x) + \nabla_{ξ} E_{n} (ξ) u_{n} (ξ, x) . \end{matrix}$ Taking inner product with $u_{n} (ξ, \cdot)$ yields $\begin{matrix} (2.11) & \nabla_{ξ} E_{n} (ξ) = - i ⟨ u_{n} (ξ, \cdot), \nabla_{x} u_{n} (ξ, \cdot) ⟩ + ξ . \end{matrix}$ Differentiate (2.10) with respect to $ξ$ again gives $\begin{array}{l} H_{ξ} \nabla_{ξ}^{2} u_{n} (ξ, x) + 2 (- i \nabla_{x} + ξ) \nabla_{ξ} u_{n} (ξ, x) + u_{n} (ξ, x) I \\ (2.12) & = E_{n} (ξ) \nabla_{ξ}^{2} u_{n} (ξ, x) + 2 \nabla_{ξ} E_{n} (ξ) \nabla_{ξ} u_{n} (ξ, x) + E_{n} (ξ) \nabla_{ξ}^{2} u_{n} (ξ, x) . \end{array}$ Taking inner product with $u_{n} (ξ, \cdot)$ , one gets $\begin{array}{l} ⟨ u_{n} (ξ, \cdot), - i \nabla_{x} \nabla_{ξ} u_{n} (ξ, x) ⟩ + ξ ⟨ u_{n} (ξ, \cdot), \nabla_{ξ} u_{n} (ξ, \cdot) ⟩ + I / 2 \\ (2.13) & = \nabla_{ξ} E_{n} (ξ) ⟨ u_{n} (ξ, \cdot), \nabla_{ξ} u_{n} (ξ, x) ⟩ + \frac{1}{2} E_{n} (ξ) ⟨ u_{n} (ξ, \cdot), \nabla_{ξ}^{2} u_{n} (ξ, \cdot) ⟩ . \end{array}$ These identities (2.11) and (2.13) will be useful later.

We shall now introduce the windowed Bloch transform. This is an analog of the windowed Fourier transform (also known as the short time Fourier transform) widely used in time-frequency signal analysis.

Definition 2.1.
The windowed Bloch transform $W : L^{2} (R^{d}) \to L^{2} {(Ω)}^{N}$ is defined as $\begin{matrix} (2.14) & {(W f)}_{n} (q, p) = \frac{2^{d / 4}}{{(2 π)}^{3 d / 4}} ⟨ u_{n} (p, \cdot) G_{q, p}, f ⟩ = \frac{2^{d / 4}}{{(2 π)}^{3 d / 4}} \int_{R^{d}} {\bar{u}}_{n} (p, x) {\bar{G}}_{q, p} (x) f (x) d x, \end{matrix}$ where $G_{q, p}$ is a Gaussian centered at $(q, p) \in Ω$ , given by $\begin{matrix} (2.15) & G_{q, p} (x) = exp (- \frac{1}{2} | x - q |^{2} + i p \cdot (x - q)) . \end{matrix}$ The adjoint operator $W^{} : L^{2} {(Ω)}^{N} \to L^{2} (R^{d})$ is then $\begin{matrix} (2.16) & (W^{} g) (x) = \frac{2^{d / 4}}{{(2 π)}^{3 d / 4}} \sum_{n = 1}^{\infty} \iint_{Ω} u_{n} (p, x) G_{q, p} (x) g_{n} (q, p) d q d p . \end{matrix}$
Proposition 2.2.
The windowed Bloch transform and its adjoint satisfies $\begin{matrix} (2.17) & W^{} W = {Id}_{L^{2} (R^{d})} . \end{matrix}$
Remark.
Similar to the windowed Fourier transform, the representation given by the windowed Bloch transform is redundant, so that $W W^{} \neq {Id}_{L^{2} {(Ω)}^{N}}$ . The normalization constant in the definition of $W$ is also due to this redundancy.
Proof.
Fix a $f \in L^{2} (R^{d})$ , by definition, we have $\begin{matrix} (W^{} W f) (x) & = \frac{2^{d / 2}}{{(2 π)}^{3 d / 2}} \sum_{n = 1}^{\infty} \iint_{Ω} u_{n} (p, x) G_{q, p} (x) ⟨ G_{q, p} u_{n} (p, \cdot), f ⟩ d q d p \\ = \frac{2^{d / 2}}{{(2 π)}^{3 d / 2}} \sum_{n = 1}^{\infty} \iint_{Ω} \int_{R^{d}} u_{n} (p, x) G_{q, p} (x) {\bar{G}}_{q, p} (y) {\bar{u}}_{n} (p, y) f (y) d y d q d p . \end{matrix}$ Let us integrate in $q$ first. $\begin{matrix} \int_{R^{d}} G_{q, p} (x) {\bar{G}}_{q, p} (y) d q & = e^{i p \cdot (x - y)} \int_{R^{d}} e^{- | x - q |^{2} / 2 - | y - q |^{2} / 2} d q \\ = e^{i p \cdot (x - y)} e^{- | x - y |^{2} / 4} \int_{R^{d}} exp (- {| q - \frac{x + y}{2} |}^{2}) d q \\ = π^{d / 2} e^{i p \cdot (x - y)} e^{- | x - y |^{2} / 4} . \end{matrix}$ Hence, denoting ${\tilde{f}}_{x} (y) = e^{- | x - y |^{2} / 4} f (y)$ , we have $\begin{matrix} (W^{} W f) (x) & = \frac{1}{{(2 π)}^{d}} \sum_{n = 1}^{\infty} \int_{Γ^{}} \int_{R^{d}} u_{n} (p, x) e^{i p \cdot (x - y)} e^{- | x - y |^{2} / 4} {\bar{u}}_{n} (p, y) f (y) d y d p \\ = \frac{1}{{(2 π)}^{d}} \sum_{n = 1}^{\infty} \int_{Γ^{}} \int_{R^{d}} u_{n} (p, x) e^{i p \cdot (x - y)} {\bar{u}}_{n} (p, y) {\tilde{f}}_{x} (y) d y d p \\ \overset{(2.4)}{=} {\tilde{f}}_{x} (x) = e^{- | x - x |^{2} / 4} f (x) = f (x) . \end{matrix}$ □

The previous proposition motivates us to consider the contribution of each band to the reconstruction formulae (2.17). This gives to the operator $Π_{n}^{W} : L^{2} (R^{d}) \to L^{2} (R^{d})$ for each $n \in N$ $\begin{matrix} (2.18) & (Π_{n}^{W} f) (x) = \frac{2^{d / 4}}{{(2 π)}^{3 d / 4}} \iint_{Ω} u_{n} (p, x) G_{q, p} (x) {(W f)}_{n} (q, p) d q d p . \end{matrix}$ It follows from (2.17) that $\sum_{n} Π_{n}^{W} = {Id}_{L^{2} (R^{d})}$ , while $Π_{n}^{W}$ is not projection due to the redundancy of windowed Bloch transform.
3. Formulation and main results

Let us start with fixing some more notations. We will switch between physical domain and phase space in the FGA formulation. For clarity, we will use $x, y \in R^{d}$ as spatial variables, $(q, p) \in R^{2 d}$ as phase space variables. The capital letters $X$ and $Y$ are shorthand notations for $X = x / ε$ and $Y = y / ε$ .

We define an effective (classical) Hamiltonian corresponding to each energy band by $\begin{matrix} (3.1) & h_{n} (q, p) = E_{n} (p) + U (q) . \end{matrix}$ The associated Hamiltonian flow $κ_{n} (t) = (Q_{n} (t, q, p), P_{n} (t, q, p))$ solves $\begin{matrix} (3.2) & \{\begin{matrix} \frac{d Q_{n}}{d t} = \nabla_{P_{n}} h_{n} (Q_{n}, P_{n}), \\ \frac{d P_{n}}{d t} = - \nabla_{Q_{n}} Q_{n}, P_{n}) \end{matrix} \end{matrix}$ on Ω with initial conditions $Q_{n} (0, q, p) = q$ and $P_{n} (0, q, p) = p$ .

From now on, we will use the short hand notation $(Q_{n}, P_{n})$ for $(Q_{n} (t, q, p), P_{n} (t, q, p))$ . For the long time existence of the Hamiltonian flow (3.2), we will assume that the external potential $U (x)$ is subquadratic as below.

Definition 3.1.
A potential U is called subquadratic, if $‖ \partial_{x}^{α} U (x) ‖_{L^{\infty}}$ is finite for all multi-index $| α | ⩾ 2$ .
Remark.
As a result, since the domain $Γ^{*}$ for $p$ is bounded, the Hamiltonian $h_{n}$ is also subquadratic.

The frozen Gaussian approximation will be formulated by the following Fourier integral operator.
Definition 3.2 (Fourier Integral Operator).

For $u \in S (R^{2 d} \times Ω, C)$ and $φ \in S (R^{d}, C)$ we define the Fourier Integral Operator with symbol u by the oscillatory integral $\begin{matrix} (3.3) & [I^{ε} (u) φ] (x) = \frac{1}{{(2 π ε)}^{3 d / 2}} \iint_{Ω} \int_{R^{d}} e^{\frac{i}{ε} Φ (t, x, y, q, p)} u (x, y, q, p) φ (y) d y d q d p \end{matrix}$ where the complex valued phase function $Φ (t, x, y, q, p)$ is given by $\begin{array}{l} Φ (t, x, y, q, p) \\ (3.4) & = S (t, q, p) - p \cdot (y - q) + P \cdot (x - Q) + \frac{i}{2} | y - q |^{2} + \frac{i}{2} | x - Q |^{2} \end{array}$ and $S (t, q, p)$ is a real-valued action function associated to κ satisfying $\begin{matrix} (3.5) & \begin{matrix} \nabla_{q} S (t, q, p) = - p + \nabla_{q} Q \cdot P, \\ \nabla_{p} S (t, q, p) = \nabla_{p} Q \cdot P . \end{matrix} \end{matrix}$

Note that if $κ (t) = κ_{n} (t)$ , the action $S_{n} (t, q, p)$ can be obtained by solving the evolution equation $\begin{matrix} (3.6) & \frac{d S_{n}}{d t} = P_{n} \cdot \nabla_{P_{n}} h_{n} (Q_{n}, P_{n}) - h_{n} (Q_{n}, P_{n}), \end{matrix}$ with initial condition $S_{n} (0, q, p) = 0$ .

We are now ready to formulate the frozen Gaussian approximation. The FGA approximates the solution of the Schrödinger equation (1.1) on the nth band to the leading order by $\begin{matrix} (3.7) & ψ_{FGA}^{ε} (t, x) = [I^{ε} (a_{n, 0} (t, q, p) u_{n} (P_{n}, \frac{x}{ε}) {\bar{u}}_{n} (p, \frac{y}{ε})) ψ_{0}^{ε}] (x), \end{matrix}$ where $ψ_{0}^{ε}$ is the initial condition. More explicitly, at time t, $ψ_{FGA}^{ε}$ is given by $\begin{array}{rcl} ψ_{FGA}^{ε} (t, x) & = & \frac{1}{{(2 π ε)}^{3 d / 2}} \iint_{Ω} a_{n, 0} (t, q, p) e^{i S_{n} (t, q, p) / ε} G_{Q_{n}, P_{n}}^{ε} (x) u_{n} (P_{n}, x / ε) \\ (3.8) & \cdot ⟨ G_{q, p}^{ε} u_{n} (p, \cdot / ε), ψ_{0}^{ε} ⟩ d q d p . \end{array}$ Here and in the sequel, we use the short-hand notation for Gaussians with semiclassical scaling $\begin{matrix} (3.9) & G_{q, p}^{ε} (x) : = exp (- \frac{| x - q |^{2}}{2 ε} + i \frac{p \cdot (x - q)}{ε}), \end{matrix}$ where the subscripts $(q, p)$ indicate the center of the Gaussian in phase space. Note that the semiclassical Fourier transform of $G_{q, p}^{ε}$ is $\begin{array}{rcl} {\hat{G}}_{q, p}^{ε} (ξ) & = & \frac{1}{{(2 π ε)}^{d / 2}} \int_{R^{d}} G_{q, p}^{ε} (x) e^{- i ξ \cdot x / ε} d x \\ (3.10) & = & exp (- \frac{| ξ - p |^{2}}{2 ε} + i \frac{q \cdot (ξ - p)}{ε}) . \end{array}$ Correspondingly, the semiclassical windowed Bloch transform $W^{ε} : L^{2} (R^{d}) \to L^{2} {(Ω)}^{N}$ is defined as $\begin{array}{rcl} {(W^{ε} f)}_{n} (q, p) & = & \frac{2^{d / 4}}{{(2 π ε)}^{3 d / 4}} ⟨ u_{n} (p, \cdot / ε) G_{q, p}^{ε}, f ⟩ \\ (3.11) & = & \frac{2^{d / 4}}{{(2 π ε)}^{3 d / 4}} \int_{R^{d}} {\bar{u}}_{n} (p, x / ε) {\bar{G}}_{q, p}^{ε} (x) f (x) d x . \end{array}$ Similarly we also have the operator $Π_{n}^{W, ε} : L^{2} (R^{d}) \to L^{2} (R^{d})$ for each $n \in N$ with semiclassical scaling $\begin{matrix} (3.12) & (Π_{n}^{W, ε} f) (y) = \frac{2^{d / 4}}{{(2 π ε)}^{3 d / 4}} \iint_{Ω} u_{n} (ξ, y / ε) G_{x, ξ}^{ε} (y) {(W^{ε} f)}_{n} (x, ξ) d x d ξ . \end{matrix}$ It follows from (2.17) and a change of variable that $\sum_{n} Π_{n}^{W, ε} = {Id}_{L^{2} (R^{d})}$ .

The only term in (3.8) that remains to be specified is the amplitude $a_{n, 0} (t, q, p)$ . It solves the evolution equation $\begin{array}{rcl} \partial_{t} a_{n, 0} & = & - i a_{n, 0} A_{n} (P_{n}) \cdot \nabla U (Q_{n}) + \frac{1}{2} a_{n, 0} tr (\partial_{z} P_{n} \nabla^{2} E (P_{n}) {(Z_{n})}^{- 1}) \\ (3.13) & - \frac{i}{2} a_{n, 0} tr (\partial_{z} Q_{n} \nabla^{2} U (Q_{n}) {(Z_{n})}^{- 1}), \end{array}$ with initial conditions $a_{n, 0} (0, q, p) = 2^{d / 2}$ for each $(q, p)$ and we recall that $A_{n} (ξ) = ⟨ u_{n} (ξ, \cdot), i \nabla_{ξ} u_{n} (ξ, \cdot) ⟩$ is the Berry phase. Here the matrix Z associated with the Hamiltonian flow $κ_{n} (t)$ is defined by $\begin{matrix} (3.14) & Z_{n} (t, q, p) : = \partial_{z} (Q_{n} + i P_{n}), \end{matrix}$ where $\partial_{z} : = \partial_{q} - i \partial_{p}$ .

We now state the main results of this work.

Theorem 3.1.
Assume that the nth Bloch band $E_{n} (ξ)$ does not intersect any other Bloch bands for all $ξ \in Γ^{}$ and the Hamiltonian* $h_{n} (x, ξ)$ is subquadratic. Let $U_{t}^{ε}$ be the propagator of the time-dependent Schrödinger equation ( 1.1 ) with initial condition $ψ_{0}^{ε} \in L^{2} (R^{d})$ . Then for any given T, $0 ⩽ t ⩽ T$ and sufficiently small ε, we have $\begin{matrix} (3.15) & sup_{0 ⩽ t ⩽ T} {‖ U_{t}^{ε} (Π_{n}^{W, ε} ψ_{0}^{ε}) - ψ_{FGA}^{ε} ‖}_{L^{2}} ⩽ C_{T, n} ε {‖ ψ_{0}^{ε} ‖}_{L^{2}} . \end{matrix}$
Remark.
Note that the FGA solution $ψ_{FGA}^{ε}$ approximates the time evolution of $Π_{n}^{W, ε} ψ_{0}^{ε}$ , which is the nth band contribution to the initial condition in the reconstruction formulae (2.17). In particular, if the initial condition is concentrated on the nth band in the sense that $ψ_{0}^{ε} = Π_{n}^{W, ε} ψ_{0}^{ε}$ , the theorem states that the solution to (1.1) is approximated by the FGA solution with $O (ε)$ error.
Remark.
We can also construct higher order approximations by replacing the term $a_{n, 0} u_{n} (P_{n}, x / ε)$ with an ε-expansion of the form $b_{n, 0} + ε b_{n, 1} + ε^{2} b_{n, 2} + \dots + ε^{N - 1} b_{n, N - 1}$ where $b_{n, 0} = a_{n, 0} u_{n} (P_{n}, x / ε)$ . This will give an approximate solution $ψ_{FGA}^{ε, N}$ to $O (ε^{N})$ accuracy. In this paper we shall focus on the first order approximation and omit the formulation and proof for higher orders.
Remark.
Let us also remark that while we take the more explicit approach of using Bloch waves in a modified FGA ansatz for periodic media, as in (3.8). The same approximation can be also derived by first projecting the whole Schrödinger equation using a super-adiabatic projection as developed in [30,31] and then apply the frozen Gaussian approximation to the resulting dynamics. We will not go into the details in this work.

The proof of Theorem 3.1 is given in Section 4. By linearity of (1.1), we have the following more general statement, as an easy corollary from Theorem 3.1. Theorem 3.2.
Assume that the first N Bloch bands $E_{n} (ξ)$ , $n = 1, \dots, N$ do not intersect and are separated from the other bands for all $ξ \in Γ^{}$ ; and assume that the Hamiltonian* $h_{n} (x, ξ)$ is subquadratic. Let $U_{t}^{ε}$ be the propagator of the time-dependent Schrödinger equation ( 1.1 ) with initial condition $ψ_{0}^{ε} \in L^{2} (R^{d})$ . Then for any given T, $0 ⩽ t ⩽ T$ and sufficiently small ε, we have $\begin{array}{l} sup_{0 ⩽ t ⩽ T} {‖ U_{t}^{ε} ψ_{0}^{ε} - \sum_{n = 1}^{N} I^{ε} (a_{n, 0} u_{n} (P_{n}, x / ε) {\bar{u}}_{n} (p, y / ε)) ψ_{0}^{ε} ‖}_{L^{2}} \\ (3.16) & ⩽ C_{T, N} ε {‖ ψ_{0}^{ε} ‖}_{L^{2}} + {‖ ψ_{0}^{ε} - \sum_{n = 1}^{N} Π_{n}^{W, ε} ψ_{0}^{ε} ‖}_{L^{2}} . \end{array}$
Proof.
Taking the short-hand notation $ψ_{0, n}^{ε} = Π_{n}^{W, ε} ψ_{0}^{ε}$ and $V_{t, n}^{ε} = I^{ε} (a_{n, 0} u_{n} (P_{n}, x / ε) {\bar{u}}_{n} (p, y / ε))$ , we have $\begin{array}{rcl} {‖ U_{t}^{ε} ψ_{0}^{ε} - \sum_{n = 1}^{N} V_{t, n}^{ε} ψ_{0}^{ε} ‖}_{L^{2}} & = & {‖ U_{t}^{ε} (\sum_{n = 1}^{\infty} ψ_{0, n}^{ε}) - \sum_{n = 1}^{N} V_{t, n}^{ε} ψ_{0}^{ε} ‖}_{L^{2}} \\ = & {‖ U_{t}^{ε} (\sum_{n = 1}^{N} ψ_{0, n}^{ε}) + U_{t}^{ε} (\sum_{n = N + 1}^{\infty} ψ_{0, n}^{ε}) - \sum_{n = 1}^{N} V_{t, n}^{ε} ψ_{0}^{ε} ‖}_{L^{2}} \\ ⩽ & {‖ U_{t}^{ε} (\sum_{n = 1}^{N} ψ_{0, n}^{ε}) - \sum_{n = 1}^{N} V_{t, n}^{ε} ψ_{0}^{ε} ‖}_{L^{2}} + {‖ U_{t}^{ε} (\sum_{n = N + 1}^{\infty} ψ_{0, n}^{ε}) ‖}_{L^{2}} \\ \overset{(3.15)}{⩽} & \sum_{n = 1}^{N} C_{T, n} ε {‖ ψ_{n, 0}^{ε} ‖}_{L^{2}} + {‖ \sum_{n = N + 1}^{\infty} ψ_{0}^{ε} ‖}_{L^{2}} \\ ⩽ & C_{T, N} ε {‖ ψ_{0}^{ε} ‖}_{L^{2}} + {‖ ψ_{0}^{ε} - \sum_{n = 1}^{N} Π_{n}^{W, ε} ψ_{0}^{ε} ‖}_{L^{2}} . \end{array}$ □

4. Analysis of frozen Gaussian approximation in periodic media

4.1. Initial condition

Let us first study the initial condition for the frozen Gaussian approximation. At time $t = 0$ , observe that by setting $t = 0$ in (3.8) we have $\begin{matrix} ψ_{FGA, n}^{ε} (0, x) = \frac{2^{d / 2}}{{(2 π ε)}^{3 d / 2}} \iint_{Ω} u_{n} (p, x / ε) G_{q, p}^{ε} (x) ⟨ G_{q, p}^{ε} u_{n} (p, \cdot / ε), ψ_{0}^{ε} ⟩ d q d p = Π_{n}^{W, ε} ψ_{0}^{ε} \end{matrix}$ by definition of the operator $Π_{n}^{W, ε}$ . Hence, the FGA solution matches $Π_{n}^{W, ε} ψ_{0}^{ε}$ at $t = 0$ .

4.2. Estimates of the Hamiltonian flows

To control the error for $t > 0$ , we collect here some preliminary results on the estimate of quantities associated with the Hamiltonian flows. We will assume throughout the rest of the paper that the assumptions of Theorem 3.1 hold for a fixed Bloch band n.

The following notation is useful in the proof. For $u \in C^{\infty} (Ω, C)$ , we define for $k \in N$ , $\begin{matrix} (4.1) & M_{k} [u] = max_{| α_{q} | + | α_{p} | ⩽ k} sup_{(q, p) \in Ω} | \partial_{q}^{α_{q}} \partial_{p}^{α_{p}} u (q, p) | \end{matrix}$ where $α_{q}$ and $α_{p}$ are multi-indices corresponding to $q$ and $p$ , respectively.

Definition 4.1 (Canonical Transformation).

Let $κ : R^{2 d} \to R^{2 d}$ be a differentiable map $κ (q, p) = (Q (q, p), P (q, p))$ and denote the Jacobian matrix as $\begin{matrix} (4.2) & (F) = (\begin{matrix} {(\partial_{q} Q)}^{T} (q, p) & {(\partial_{p} Q)}^{T} (q, p) \\ {(\partial_{q} P)}^{T} (q, p) & {(\partial_{p} P)}^{T} (q, p) \end{matrix}) . \end{matrix}$ We say κ is a canonical transformation if F is symplectic for any $(q, p) \in R^{2 d}$ , i.e. $\begin{matrix} (4.3) & {(F)}^{T} (\begin{matrix} 0 & {Id}_{d} \\ - {Id}_{d} & 0 \end{matrix}) F = (\begin{matrix} 0 & {Id}_{d} \\ - {Id}_{d} & 0 \end{matrix}) . \end{matrix}$

It is easy to check by the definition that the map $κ_{n} (t) : R^{2 d} \to R^{2 d}$ defined by $(q, p) \to (Q_{n} (t, q, p), P_{n} (t, q, p))$ solving (3.2) is a canonical transformation.

Proposition 4.2.
We have for all $k ⩾ 0$ $\begin{matrix} (4.4) & sup_{t \in [0, T]} M_{k} [F_{n} (t)] < \infty, sup_{t \in [0, T]} M_{k} [\frac{d}{d t} F_{n} (t)] < \infty . \end{matrix}$
Proof.
Differentiating $F_{n} (t, q, p)$ with respect to t gives $\begin{matrix} (4.5) & \frac{d}{d t} F_{n} (t, q, p) = (\begin{matrix} \partial_{P} \partial_{Q} h_{n} & \partial_{P} \partial_{P} h_{n} \\ - \partial_{Q} \partial_{Q} h_{n} & - \partial_{Q} \partial_{P} h_{n} \end{matrix}) F_{n} (t, q, p) . \end{matrix}$ By our assumption that U is subquadratic on $R^{d}$ and since $E_{n} \in C^{\infty} (Γ^{})$ , there exists a constant C independent of $(q, p)$ such that $\begin{matrix} (4.6) & \frac{d}{d t} | F_{n} (t, q, p) | = |(\begin{matrix} \partial_{P} \partial_{Q} h_{n} & \partial_{P} \partial_{P} h_{n} \\ - \partial_{Q} \partial_{Q} h_{n} & - \partial_{Q} \partial_{P} h_{n} \end{matrix})| | F_{n} (t, q, p) | ⩽ C | F_{n} (t, q, p) | \end{matrix}$ with $| F_{n} (0) | = | {Id}_{2 d} |$ . By an application of Gronwall’s inequality, we obtain $\begin{matrix} (4.7) & | F_{n} (t) | ⩽ e^{C | t |} . \end{matrix}$ Differentiating (4.5) with respect to $(q, p)$ yields $\begin{array}{rcl} \frac{d}{d t} \partial_{q}^{α_{q}} \partial_{p}^{α_{p}} F_{n} (t, q, p) & = & \sum_{β_{q} ⩽ α_{q}, β_{p} ⩽ α_{p}} (\binom{α_{q}}{β_{q}}) (\binom{α_{p}}{β_{p}}) \partial_{q}^{β_{q}} \partial_{p}^{β_{p}} (\begin{matrix} \partial_{P} \partial_{Q} h_{n} & \partial_{P} \partial_{P} h_{n} \\ - \partial_{Q} \partial_{Q} h_{n} & - \partial_{Q} \partial_{P} h_{n} \end{matrix}) \\ (4.8) & \times \partial_{q}^{α_{q} - β_{q}} \partial_{p}^{α_{p} - β_{p}} F_{n} (t, q, p) . \end{array}$ Our estimate now follows by induction. □

Recall that the matrix $Z_{n} (t, q, p)$ is defined by $\begin{matrix} (4.9) & Z_{n} (t, q, p) : = \partial_{z} (Q_{n} (t, q, p) + i P_{n} (t, q, p)) = (\partial_{q} - i \partial_{p}) (Q_{n} (t, q, p) + i P_{n} (t, q, p)) . \end{matrix}$ We have the following. It follows the same proof of [24, Proposition 3.5], which we reproduce here for completeness. Proposition 4.3.
$Z_{n} (t, q, p)$ is invertible for* $(q, p) \in Ω$ . Moreover, for each $k \in N$ , $\begin{matrix} (4.10) & M_{k} [{(Z_{n} (t))}^{- 1}] < \infty . \end{matrix}$
Proof.
$Z_{n} (t, q, p)$ inherits the property that $M_{k} (Z_{n} (t, q, p)) < \infty$ from the same estimate for $F_{n} (t, q, p)$ . Moreover, we have $\begin{array}{rcl} Z_{n} {(Z_{n})}^{} (t, q, p) & = & (\begin{matrix} i {Id}_{d} & {Id}_{d} \end{matrix}) {(F_{n})}^{T} (t, q, p) (\begin{matrix} {Id}_{d} & - i {Id}_{d} \\ i {Id}_{d} & {Id}_{d} \end{matrix}) F_{n} (t, q, p) (\begin{matrix} - i {Id}_{d} \\ {Id}_{d} \end{matrix}) \\ = & (\begin{matrix} i {Id}_{d} & {Id}_{d} \end{matrix}) ({(F_{n})}^{T} (F_{n})) (t, q, p) (\begin{matrix} - i {Id}_{d} \\ {Id}_{d} \end{matrix}) \\ + (\begin{matrix} i {Id}_{d} & {Id}_{d} \end{matrix}) {(F_{n})}^{T} (t, q, p) (\begin{matrix} 0 & - i {Id}_{d} \\ i {Id}_{d} & 0 \end{matrix}) F_{n} (t, q, p) (\begin{matrix} - i {Id}_{d} \\ {Id}_{d} \end{matrix}) \\ (4.11) & = & (\begin{matrix} i {Id}_{d} & {Id}_{d} \end{matrix}) ({(F_{n})}^{T} F_{n}) (t, q, p) (\begin{matrix} - i {Id}_{d} \\ {Id}_{d} \end{matrix}) + 2 {Id}_{d} . \end{array}$

This calculation shows that, since ${(F_{n} (t))}^{T} F_{n} (t)$ is semi-positive definite, for any $v \in C^{2 d}$ , $\begin{matrix} (4.12) & v^{} Z_{n} (t) {(Z_{n} (t))}^{} v ⩾ 2 | v |^{2} . \end{matrix}$ Therefore $Z_{n} (t, q, p)$ is invertible and $det (Z_{n} (t))$ is uniformly bounded away from 0 for all $q$ and $p$ , so by representing ${(Z_{n})}^{- 1} (t, q, p)$ by minors, $M_{k} ({(Z_{n})}^{- 1} (t, q, p)) < \infty$ , as $M_{k} (Z_{n} (t, q, p))$ is. □
Proposition 4.4.
For each* $k \in N$ , $\begin{matrix} (4.13) & sup_{t \in [0, T]} M_{k} [u_{n} (P_{n}, x)] < \infty . \end{matrix}$
Proof.
$u_{n} (P_{n}, x)$ is smooth on the compact set $Γ^{} \times Γ$ since the nth band is separated from the rest of the spectrum (see e.g., [35, Sec XIII.16]). Thus $u_{n} (P_{n}, x)$ is uniformly bounded on $Γ^{} \times Γ$ and hence $Γ^{*} \times R^{d}$ due to periodicity. We also see from Proposition 4.2 that the derivatives of $u_{n} (P_{n}, x)$ are also bounded. Thus, $M_{k} [u_{n} (P_{n}, x)] < \infty$ for any finite time t. □

4.3. Higher order asymptotic solution

To prove the theorem, we will need to construct a solution to the Schrödinger equation that is accurate up to $O (ε^{2})$ . The construction is based on matched asymptotic expansion. Let us fix a band $n \in N$ and consider the ansatz $\begin{matrix} (4.14) & ψ_{FGA, \infty}^{ε} = \frac{1}{{(2 π ε)}^{3 d / 2}} \iint_{Ω} b^{ε} (t, X, q, p) G_{Q_{n}, P_{n}}^{ε} e^{i S_{n} / ε} ⟨ G_{q, p}^{ε} u_{n} (p, \cdot / ε), ψ_{0} ⟩ d q d p, \end{matrix}$ where the coefficient b assumes the asymptotic expansion $\begin{array}{rcl} b^{ε} (t, X, q, p) & : = & \sum_{j = 0}^{\infty} ε^{j} b_{j} (t, X, q, p) \\ = & a_{n, 0} (t, q, p) u_{n} (P_{n}, X) \\ + ε (a_{n, 1} (t, q, p) u_{n} (P_{n}, X) + b_{n, 1}^{⊥} (t, X, q, p)) \\ + ε^{2} (a_{n, 2} (t, q, p) u_{n} (P_{n}, X) + b_{n, 2}^{⊥} (t, X, q, p)) \\ (4.15) & + \sum_{j = 3}^{\infty} ε^{j} b_{j} (t, X, q, p) \end{array}$

To determine the terms in the expansion, we will make use of the following Lemma.

Definition 4.5.
For $f = f (t, x, y, q, p)$ and $g = g (t, x, y, q, p)$ such that for any t and $x$ , $\begin{matrix} f (t, x, \cdot, \cdot, \cdot), g (t, x, \cdot, \cdot, \cdot) \in L^{\infty} (R^{d}; S (R^{d} \times Γ^{})), \end{matrix}$ we say that f and g are equivalent for the nth Bloch band, denoted as $f \sim_{n} g$ if for any $t ⩾ 0$ and $Ψ_{0} \in L^{2} (R^{d})$ $\begin{matrix} (4.16) & \iint_{Ω} \int_{R^{d}} (f - g) (t, x, y, q, p) G_{Q_{n}, P_{n}}^{ε} e^{i S_{n} (t, q, p) / ε} {\bar{G}}_{q, p}^{ε} (y) Ψ_{0} (y) d y d q d p = 0 . \end{matrix}$
Lemma 4.6.
For any d-vector function* $v (y, q, p)$ such that each component is in $L^{\infty} (R^{d}; S (R^{d} \times Γ^{}))$ $\begin{matrix} (4.17) & v (y, q, p) \cdot (x - Q_{n}) \sim_{n} - ε \partial_{z} \cdot (v Z_{n}^{- 1}), \end{matrix}$ and for any* $d \times d$ matrix function $M (y, q, p)$ such that each component is in $L^{\infty} (R^{d}; S (R^{d} \times Γ^{}))$ $\begin{array}{l} tr (M (y, q, p) {(x - Q_{n})}^{2}) \\ \sim_{n} ε tr (\partial_{z} Q_{n} M Z_{n}^{- 1}) - ε tr (\partial_{z} M (x - Q_{n}) Z_{n}^{- 1} + M (x - Q_{n}) \partial_{z} Z_{n}^{- 1}) \\ (4.18) & = ε tr (\partial_{z} Q_{n} M Z_{n}^{- 1}) + ε^{2} tr (\partial_{z} (\partial_{z} M {(Z_{n}^{- 1})}^{2}) + \partial_{z} (M \partial_{z} Z_{n}^{- 1}) Z_{n}^{- 1}) . \end{array}$ Higher order terms can be obtained recursively. In general we have for any multi-index α that* $| α | ⩾ 3$ , $\begin{matrix} (4.19) & {(x - Q_{n})}^{α} \sim_{n} O (ε^{⌊ \frac{| α | + 1}{2} ⌋}) . \end{matrix}$
Proof.
The proof of Lemma 4.6 is essentially the same as in Lemma 3 of [39] and thus is omitted here. □

We now substitute (4.14) into the Schrödinger equation. For this we first compute the time and space derivatives on $ψ_{FGA, \infty}^{ε}$ : $\begin{array}{l} i ε \partial_{t} ψ_{FGA, \infty}^{ε} \\ = \frac{1}{{(2 π ε)}^{3 d / 2}} \iint_{Ω} {i ε \partial_{t} b^{ε} - (\partial_{t} S_{n} - P_{n} \cdot \partial_{t} Q_{n} + (\partial_{t} P_{n} - i \partial_{t} Q_{n}) \cdot (x - Q_{n})) b^{ε}} \\ (4.20) & \times G_{Q_{n}, P_{n}}^{ε} e^{i S_{n} / ε} ⟨ G_{q, p}^{ε} u_{n} (p, \cdot / ε), ψ_{0} ⟩ d q d p . \\ \frac{1}{2} ε^{2} Δ ψ_{FGA, \infty}^{ε} \\ = \frac{1}{{(2 π ε)}^{3 d / 2}} \iint_{Ω} [- \frac{1}{2} {(- i \nabla_{X} + P_{n})}^{2} b^{ε} - (\nabla_{X} b^{ε} + i b^{ε} P_{n}) \cdot (x - Q_{n}) \\ + \frac{1}{2} b^{ε} | x - Q_{n} |^{2} - \frac{1}{2} ε b^{ε} d] \\ (4.21) & \times G_{Q_{n}, P_{n}}^{ε} e^{i S_{n} / ε} ⟨ G_{q, p}^{ε} u_{n} (p, \cdot / ε), ψ_{0} ⟩ d q d p . \end{array}$ Hence, after rearranging terms, we arrive at $\begin{array}{l} (i ε \partial_{t} + \frac{1}{2} ε^{2} Δ - V (X) - U (x)) ψ_{FGA, \infty}^{ε} \\ = \frac{1}{{(2 π ε)}^{3 d / 2}} \iint_{Ω} {[- \frac{1}{2} {(- i \nabla_{X} + P_{n})}^{2} - V (X) - U (x) - \partial_{t} S_{n}] b^{ε} \\ + ε (i \partial_{t} b^{ε} - \frac{1}{2} b^{ε} d) - [(\nabla_{X} b^{ε} + i b^{ε} P_{n}) + (\partial_{t} P_{n} - i \partial_{t} Q_{n}) b^{ε}] \cdot (x - Q_{n}) \\ (4.22) & + \frac{1}{2} | x - Q_{n} |^{2} b^{ε} + P_{n} \cdot \partial_{t} Q_{n} b^{ε}} G_{Q_{n}, P_{n}}^{ε} e^{i S_{n} / ε} ⟨ G_{q, p}^{ε} u_{n} (p, \cdot / ε), ψ_{0} ⟩ d q d p . \end{array}$ Define $\begin{array}{rcl} f (t, x, y, q, p) & = & {[- \frac{1}{2} {(- i \nabla_{X} + P_{n})}^{2} - V (X) - U (x) - \partial_{t} S_{n}] b^{ε} \\ + ε (i \partial_{t} b^{ε} - \frac{1}{2} b^{ε} d) - [(\nabla_{X} b^{ε} + i b^{ε} P_{n}) + (\partial_{t} P_{n} - i \partial_{t} Q_{n}) b^{ε}] \cdot (x - Q_{n}) \\ (4.23) & + \frac{1}{2} | x - Q_{n} |^{2} b^{ε} + P_{n} \cdot \partial_{t} Q_{n} b^{ε}} {\bar{u}}_{n} (p, Y), \end{array}$ then we can write $\begin{array}{l} (i ε \partial_{t} + \frac{1}{2} ε^{2} Δ - V (X) - U (x)) ψ_{FGA, \infty}^{ε} \\ (4.24) & = \frac{1}{{(2 π ε)}^{3 d / 2}} \iint_{Ω} \int_{R^{d}} f (t, x, y, p, q) G_{Q_{n}, P_{n}}^{ε} e^{i S_{n} / ε} {\bar{G}}_{q, p}^{ε} (y) ψ_{0} (y) d y d q d p . \end{array}$ Applying Lemma 4.6 and adding and subtracting $U (Q_{n})$ , we get $\begin{array}{rcl} f & \sim_{n} & (- \frac{1}{2} {(- i \nabla_{X} + P_{n})}^{2} - V (X) - (U (x) - U (Q_{n})) - \partial_{t} S_{n}) b^{ε} {\bar{u}}_{n} (p, Y) \\ + ε (i \partial_{t} b^{ε} - \frac{1}{2} b^{ε} d) {\bar{u}}_{n} (p, Y) \\ + ε \partial_{z} ([(\nabla_{X} b^{ε} + i b^{ε} P_{n}) + (\partial_{t} P_{n} - i \partial_{t} Q_{n}) b^{ε}] {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) \\ + ε \frac{1}{2} b^{ε} tr [\partial_{z} Q_{n} Z_{n}^{- 1}] {\bar{u}}_{n} (p, Y) + ε^{2} \frac{1}{2} tr [\partial_{z} (\partial_{z} (b^{ε} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) Z_{n}^{- 1})] \\ (4.25) & + P_{n} \cdot \partial_{t} Q_{n} b^{ε} {\bar{u}}_{n} (p, Y) - U (Q_{n}) b^{ε} {\bar{u}}_{n} (p, Y) \end{array}$ Use the Taylor expansion of $U (x)$ about $Q_{n}$ $\begin{array}{l} (U (x) - U (Q_{n})) \\ = \nabla U (Q_{n}) (x - Q_{n}) + \frac{1}{2!} \nabla^{2} U (Q_{n}) {(x - Q_{n})}^{2} \\ + \frac{1}{3!} \nabla^{3} U (Q_{n}) {(x - Q_{n})}^{3} + \frac{1}{4!} \nabla^{4} U (Q_{n}) {(x - Q_{n})}^{4} \\ (4.26) & + \sum_{| α | = 5} R_{α} (x) {(x - Q_{n})}^{α} \end{array}$ with $\begin{matrix} (4.27) & R_{α} (x) = \frac{| α |}{5!} \int_{0}^{1} {(1 - τ)}^{| α | - 1} \partial_{Q_{n}}^{α} U (Q_{n} + τ (x - Q_{n})) d τ . \end{matrix}$ From now on, let us denote the remainder term in (4.26) by $R (x, q, p)$ .

Applying Lemma 4.6 again to (4.25) together with (4.26), we obtain $\begin{array}{rcl} f & \sim_{n} & (- \frac{1}{2} {(- i \nabla_{X} + P_{n})}^{2} - V (X) - \partial_{t} S_{n}) b^{ε} {\bar{u}}_{n} (p, Y) \\ + P_{n} \cdot \partial_{t} Q_{n} b^{ε} {\bar{u}}_{n} (p, Y) - U (Q_{n}) b^{ε} {\bar{u}}_{n} (p, Y) \\ + ε (i \partial_{t} b^{ε} - \frac{1}{2} b^{ε} d) {\bar{u}}_{n} (p, Y) + ε \partial_{z} (\nabla U (Q_{n}) b^{ε} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) \\ + ε \partial_{z} ([(\nabla_{X} b^{ε} + i b^{ε} P_{n}) + (\partial_{t} P_{n} - i \partial_{t} Q_{n}) b^{ε}] {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) \\ + ε \frac{1}{2!} tr [\partial_{z} Q_{n} (I - \nabla^{2} U (Q_{n})) b^{ε} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}] \\ + ε^{2} \frac{1}{2!} tr [\partial_{z} (\partial_{z} ((I - \nabla^{2} U (Q_{n})) b^{ε} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) Z_{n}^{- 1})] \\ + ε^{2} \frac{2}{3!} tr [\partial_{z} (\partial_{z} Q_{n} \nabla^{3} U (Q_{n}) b^{ε} {\bar{u}}_{n} (p, Y) {(Z_{n}^{- 1})}^{2})] \\ + ε^{2} \frac{1}{3!} tr [\partial_{z} Q_{n} \partial_{z} (\nabla^{3} U (Q_{n}) b^{ε} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) Z_{n}^{- 1}] \\ (4.28) & - ε^{2} \frac{3}{4!} tr [{(\partial_{z} Q_{n})}^{2} \nabla^{4} U (Q_{n}) b^{ε} {\bar{u}}_{n} (p, Y) {(Z_{n}^{- 1})}^{2}] + R (x, q, p) b^{ε} {\bar{u}}_{n, p} (Y) . \end{array}$ Let us define three operators $L_{0}^{n}$ , $L_{1}^{n}$ , and $L_{2}^{n}$ acting on $Φ = Φ (t, x, y, q, p)$ by $\begin{array}{l} L_{0}^{n} (Φ) : = (- \frac{1}{2} {(- i \nabla_{X} + P_{n})}^{2} - V (X) - \partial_{t} S_{n}) Φ \\ + P_{n} \cdot \partial_{t} Q_{n} Φ - U (Q_{n}) Φ \\ (4.29) & = (- H_{P_{n}} + E_{n} (P_{n})) Φ, \\ L_{1}^{n} (Φ) : = (i \partial_{t} Φ - \frac{1}{2} Φ d) + \partial_{z} (\nabla U (Q_{n}) Φ Z_{n}^{- 1}) \\ + \partial_{z} ([(\nabla_{X} Φ + i Φ P_{n}) + (\partial_{t} P_{n} - i \partial_{t} Q_{n}) Φ] Z_{n}^{- 1}) \\ (4.30) & + \frac{1}{2!} tr [\partial_{z} Q_{n} (I - \nabla^{2} U (Q_{n})) Φ Z_{n}^{- 1}], \end{array}$ and $\begin{array}{l} L_{2}^{n} (Φ) : = & \frac{1}{2!} tr [\partial_{z} (\partial_{z} ((I - \nabla^{2} U (Q_{n})) Φ Z_{n}^{- 1}) Z_{n}^{- 1})] \\ + \frac{2}{3!} tr [\partial_{z} (\partial_{z} Q_{n} \nabla^{3} U (Q_{n}) Φ {(Z_{n}^{- 1})}^{2})] \\ + \frac{1}{3!} tr [\partial_{z} Q_{n} \partial_{z} (\nabla^{3} U (Q_{n}) Φ Z_{n}^{- 1}) Z_{n}^{- 1}] \\ (4.31) & - \frac{3}{4!} tr [{(\partial_{z} Q_{n})}^{2} \nabla^{4} U (Q_{n}) Φ {(Z_{n}^{- 1})}^{2}] . \end{array}$ We thus arrive at $\begin{array}{l} (i ε \partial_{t} + \frac{1}{2} ε^{2} Δ - V (X) - U (x)) ψ_{FGA, \infty}^{ε} \\ = \frac{1}{{(2 π ε)}^{3 d / 2}} \iint_{Ω} \int_{R^{d}} {L_{0}^{n} (b^{ε} {\bar{u}}_{n} (p, Y)) + ε L_{1}^{n} (b^{ε} {\bar{u}}_{n} (p, Y)) \\ (4.32) & + ε^{2} L_{2}^{n} (b^{ε} {\bar{u}}_{n} (p, Y)) + R (x, q, p) b^{ε} {\bar{u}}_{n} (p, Y)} G_{Q_{n}, P_{n}}^{ε} e^{i S_{n} / ε} {\bar{G}}_{q, p}^{ε} (y) ψ_{0} (y) d y d q d p . \end{array}$

Note that by the choice $b_{n, 0} = a_{n, 0} u_{n} (P_{n}, X)$ , the $O (1)$ term in the integrand on the right hand side of (4.32) vanishes as $\begin{matrix} (4.33) & L_{0}^{n} (a_{n, 0} (t, q, p) u_{n} (p, X)) = a_{n, 0} (t, q, p) (- H_{P_{n}} + E_{n} (P_{n})) u_{n} (P_{n}, X) = 0 \end{matrix}$ for any $a_{n, 0}$ .
4.3.1. Leading order term $b_{n, 0}$

To determine $a_{n, 0}$ , we set the order $O (ε)$ term on the right hand side of (4.32) to zero and get $\begin{matrix} (4.34) & L_{0}^{n} (b_{n, 1} {\bar{u}}_{n} (p, Y)) = - L_{1}^{n} (b_{n, 0} {\bar{u}}_{n} (p, Y)) . \end{matrix}$ We multiply the equation by ${\bar{u}}_{n} (P_{n}, X)$ and integrate over Γ; this gives $\begin{matrix} (4.35) & \partial_{t} a_{n, 0} = \frac{1}{2} a_{n, 0} tr (\partial_{z} P_{n} (\nabla_{P_{n}}^{2} E_{n}) Z_{n}^{- 1}) - i a_{n, 0} A (P_{n}) \cdot \nabla_{Q_{n}} U - \frac{i}{2} a_{n, 0} tr (\partial_{z} Q_{n} (\nabla_{Q_{n}}^{2} U) Z_{n}^{- 1}) . \end{matrix}$ Indeed, by integration, we get (index n is suppressed) $\begin{array}{l} \int_{Γ} {\bar{u}}_{n} (P_{n}, X) (- \frac{1}{2} {(- i \nabla_{X} + P_{n})}^{2} - V (X) - \partial_{t} S_{n}) b_{1} {\bar{u}}_{n} (p, Y) d X \\ + \int_{Γ} {{\bar{u}}_{n} (P_{n}, X) (i \partial_{t} b_{0} - \frac{1}{2} b_{0} d) {\bar{u}}_{n} (p, Y) + {\bar{u}}_{n} (P_{n}, X) \partial_{z} (\nabla U (Q_{n}) b_{0} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) \\ + {\bar{u}}_{n} (P_{n}, X) \partial_{z} ([(\nabla_{X} b_{0} + i b_{0} P_{n}) + (\partial_{t} P_{n} - i \partial_{t} Q_{n}) b_{0}] {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) \\ (4.36) & + {\bar{u}}_{n} (P_{n}, X) \frac{1}{2!} tr [\partial_{z} Q_{n} (I - \nabla^{2} U (Q_{n})) b_{0} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}]} d X = 0 . \end{array}$ The perpendicular terms in the $b_{j}$ ’s will now drop out and we can symplify this equation to $\begin{array}{l} - ⟨ u_{n} (P_{n}, X), \partial_{z} ([i u_{n} (P_{n}, X) \nabla_{P_{n}} E_{n} - \nabla_{X} u_{n} (P_{n}, X) - i u_{n} (P_{n}, X) P_{n}] a_{0} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) ⟩ \\ + (i \partial_{t} a_{0} - a_{0} A (P_{n}) \cdot \nabla_{Q_{n}} U - \frac{d}{2} a_{0}) {\bar{u}}_{n} (p, Y) + \frac{1}{2} a_{0} tr (\partial_{z} Q_{n} (I - \nabla_{Q_{n}}^{2} U) Z_{n}^{- 1}) {\bar{u}}_{n} (p, Y) \\ (4.37) & = 0 . \end{array}$ Using (2.11), we observe that $\begin{array}{l} ⟨ u_{n} (P_{n}, X), [i u_{n} (P_{n}, X) \nabla_{P_{n}} E_{n} - \nabla_{X} u_{n} (P_{n}, X) - i u_{n} (P_{n}, X) P_{n}] \cdot \partial_{z} (a_{0} {\bar{u}}_{n} (P_{n}, Y) Z_{n}^{- 1}) ⟩ \\ (4.38) & = 0 . \end{array}$ Hence, we arrive at $\begin{array}{l} a_{0} tr (⟨ u_{n} (P_{n}, X), \partial_{z} \cdot [i u_{n} (P_{n}, X) \nabla_{P_{n}} E_{n} - \nabla_{X} u_{n} (P_{n}, X) - i u_{n} (P_{n}, X) P_{n}] ⟩ Z_{n}^{- 1}) \\ + (i \partial_{t} a_{0} - a_{0} A (P_{n}) \cdot \nabla_{Q_{n}} U - \frac{d}{2} a_{0}) + \frac{1}{2} a_{0} tr (\partial_{z} Q_{n} (I - \nabla_{Q_{n}}^{2} U) Z_{n}^{- 1}) \\ (4.39) & = 0 . \end{array}$ To further simplify the equation, observe that $\begin{array}{l} ⟨ u_{n} (P_{n}, X), \partial_{z} \cdot [i u_{n} (P_{n}, X) \nabla_{P_{n}} E_{n} - \nabla_{X} u_{n} (P_{n}, X) - i u_{n} (P_{n}, X) P_{n}] ⟩ \\ = i ⟨ u_{n} (P_{n}, X), \partial_{z} u_{n} (P_{n}, X) ⟩ (\nabla_{P_{n}} E_{n} - P_{n}) - ⟨ u_{n} (P_{n}, X), \partial_{z} \nabla_{X} u_{n} (P_{n}, X) ⟩ \\ + i (\partial_{z} \nabla_{P_{n}} E_{n} - \partial_{z} P_{n}) \\ = i \partial_{z} P_{n} ⟨ u_{n} (P_{n}, X), \partial_{P_{n}} u_{n} (P_{n}, X) ⟩ (\nabla_{P_{n}} E_{n} - P_{n}) \\ - \partial_{z} P_{n} {⟨ u_{n} (P_{n}, X), \nabla_{p} \nabla_{X} u_{n} (P_{n}, X) ⟩}_{Γ} + i \partial_{z} P_{n} (\nabla_{P_{n}}^{2} E_{n} - I) \\ (4.40) & \overset{(2.13)}{=} \frac{1}{2} i \partial_{z} P_{n} (\nabla_{P_{n}}^{2} E_{n} - I) . \end{array}$ Putting this into (4.39), we have $\begin{array}{l} \frac{1}{2} i a_{0} tr (\partial_{z} P_{n} (I - \nabla_{P_{n}}^{2} E_{n}) Z^{- 1}) + (i \partial_{t} a_{0} - a_{0} A (P_{n}) \cdot \nabla_{Q_{n}} U - \frac{d}{2} a_{0}) \\ (4.41) & + \frac{1}{2} a_{0} tr (\partial_{z} Q_{n} (I - \nabla_{Q_{n}}^{2} U) Z_{n}^{- 1})] = 0 . \end{array}$ We arrive at (4.35) finally by noting that $\begin{matrix} (4.42) & \frac{1}{2} a tr [\partial_{z} Q_{n} Z_{n}^{- 1}] + \frac{i}{2} a tr [\partial_{z} P_{n} Z_{n}^{- 1}] = \frac{1}{2} a tr [Z_{n} Z_{n}^{- 1}] = \frac{d}{2} a . \end{matrix}$

4.3.2. Next order term $b_{n, 1}$

To characterize $b_{n, 1}$ , we set the order $O (ε^{2})$ term in (4.32) to zero, we have $\begin{matrix} (4.43) & \int_{Γ} {\bar{u}}_{n} (P_{n}, X) (L_{0}^{n} (b_{n, 2} {\bar{u}}_{n} (p, Y)) + L_{1}^{n} (b_{n, 1} {\bar{u}}_{n} (p, Y)) + L_{2}^{n} (b_{n, 0} {\bar{u}}_{n} (p, Y))) d X = 0 . \end{matrix}$ Let us first derive the equation for $a_{1}$ . We start with (4.43) written in expanded form $\begin{array}{l} \int_{Γ} {\bar{u}}_{n} (P_{n}, X) {\frac{1}{2!} tr [\partial_{z} (\partial_{z} ((I - \nabla^{2} U (Q_{n})) b_{0} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) Z_{n}^{- 1})] \\ + \frac{2}{3!} tr [\partial_{z} (\partial_{z} Q_{n} \nabla^{3} U (Q_{n}) b_{0} {\bar{u}}_{n} (p, Y) {(Z_{n}^{- 1})}^{2})] \\ + \frac{1}{3!} tr [\partial_{z} Q_{n} \partial_{z} (\nabla^{3} U (Q_{n}) b_{0} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) Z_{n}^{- 1}] \\ - \frac{3}{4!} tr [{(\partial_{z} Q_{n})}^{2} \nabla^{4} U (Q_{n}) b_{0} {\bar{u}}_{n} (p, Y) {(Z_{n}^{- 1})}^{2}] \\ + (i \partial_{t} b_{1} - \frac{1}{2} b_{1} d) {\bar{u}}_{n} (p, Y) + \partial_{z} (\nabla U (Q_{n}) b_{1} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) \\ + \partial_{z} ([(\nabla_{X} b_{1} + i b_{1} P_{n}) + (\partial_{t} P_{n} - i \partial_{t} Q_{n}) b_{1}] {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}) \\ + \frac{1}{2!} tr [\partial_{z} Q_{n} (I - \nabla^{2} U (Q_{n})) b_{1} {\bar{u}}_{n} (p, Y) Z_{n}^{- 1}] \\ + (- H_{P_{n}} + E (P_{n})) b_{2} {\bar{u}}_{n} (p, Y)} d X \\ (4.44) & = 0 . \end{array}$ Making use of the Hamiltonian flow (3.2) and the identity (4.38), we arrive at $\begin{array}{l} - tr (⟨ u_{P_{n}}, \partial_{z} \cdot [u (i \nabla E_{n} (P_{n})) - \nabla_{X} u - i u P_{n}] (a_{1}) ⟩ Z_{n}^{- 1}) \\ + (i \partial_{t} a_{1} - a_{1} A (P_{n}) \cdot \nabla_{Q_{n}} U - \frac{d}{2} a_{1}) + a_{0} \frac{1}{2} tr (\partial_{z} (\partial_{z} [(I - \nabla_{Q_{n}}^{2} U) Z_{n}^{- 1}] Z_{n}^{- 1})) \\ + a_{1} \frac{1}{2} tr (\partial_{z} Q_{n} (I - \nabla_{Q_{n}}^{2} U) Z_{n}^{- 1}) + \frac{2}{3!} a_{0} tr (\partial_{z} (\partial_{z} Q_{n} \nabla_{Q_{n}}^{3} U {(Z_{n}^{- 1})}^{2})) \\ (4.45) & + \frac{1}{3!} a_{0} tr (\partial_{z} Q_{n} \partial_{z} (\nabla_{Q_{n}}^{3} U Z_{n}^{- 1}) Z_{n}^{- 1}) - \frac{3}{4!} a_{0} tr ({(\partial_{z} Q_{n})}^{2} \nabla_{Q_{n}}^{4} U {(Z_{n}^{- 1})}^{2}) = 0 . \end{array}$ Then using (4.40) and (4.42), upon simplification we obtain the equation for $a_{n, 1}$ $\begin{array}{rcl} \partial_{t} a_{n, 1} & = & - i a_{n, 1} A (P_{n}) \cdot \nabla_{Q_{n}} U + \frac{1}{2} a_{n, 1} tr (\partial_{z} P_{n} (\nabla_{P_{n}}^{2} E_{n}) Z_{n}^{- 1}) - \frac{i}{2} a_{n, 1} tr (\partial_{z} Q_{n} (\nabla_{Q_{n}}^{2} U) Z_{n}^{- 1}) \\ + \frac{i}{2} a_{n, 0} tr (\partial_{z} (\partial_{z} [(I - \nabla_{Q_{n}}^{2} U) Z_{n}^{- 1}] Z_{n}^{- 1})) + \frac{2 i}{3!} a_{n, 0} tr (\partial_{z} (\partial_{z} Q_{n} \nabla_{Q_{n}}^{3} U {(Z_{n}^{- 1})}^{2})) \\ (4.46) & + \frac{i}{3!} a_{n, 0} tr (\partial_{z} Q_{n} \partial_{z} (\nabla_{Q_{n}}^{3} U Z_{n}^{- 1}) Z_{n}^{- 1}) - \frac{3 i}{4!} a_{n, 0} tr ({(\partial_{z} Q_{n})}^{2} \nabla_{Q_{n}}^{4} U {(Z_{n}^{- 1})}^{2}) . \end{array}$

Define the operator $Q = Id - Π_{n}$ where $Π_{n}$ is the projection operator onto the nth Bloch wave. $b_{n, 1}^{⊥}$ satisfies $Π_{n} b_{n, 1}^{⊥} = ⟨ u_{n} (P_{n}, x), b_{n, 1}^{⊥} ⟩ = 0$ , and is hence determined by applying $Q$ to $L_{0}^{n} (b_{n, 1} {\bar{u}}_{n} (p, Y)) = - L_{1}^{n} (b_{n, 0} {\bar{u}}_{n} (p, Y))$ . We obtain $\begin{matrix} (4.47) & b_{n, 1}^{⊥} {\bar{u}}_{n} (p, Y) = - {(L_{0}^{n})}^{- 1} Q (L_{1}^{n} (b_{n, 0} {\bar{u}}_{n} (p, Y))) . \end{matrix}$ Note that the inverse of the operator $L_{0}^{n}$ can be defined on its range.

Thus, we have obtained the equations for $a_{n, 0}$ (4.35), $a_{n, 1}$ (4.46), and $b_{n, 1}^{⊥}$ (4.47). This can be continued to higher orders. Let us summarize the estimate of these terms in the following propositions.

Proposition 4.7.
For each $k \in N$ , the amplitudes $a_{n, 0}$ and $a_{n, 1}$ , given by ( 4.35 ) and ( 4.46 ) satisfy $\begin{matrix} (4.48) & sup_{t \in [0, T]} M_{k} [a_{n, 0}] < \infty and sup_{t \in [0, T]} M_{k} [a_{n, 1}] < \infty . \end{matrix}$
Proof.
By (4.2), (4.3) and (4.4), we see that the right hand side of (4.35) and (4.46) are bounded by some constants independent of $q$ and $p$ times $a_{n, 0}$ and $a_{n, 1}$ , respectively. An application of Gronwall’s inequality yields the result. □
Proposition 4.8.
For each $k \in N$ we have that $\begin{matrix} (4.49) & sup_{t \in [0, T]} M_{k} [b_{n, 1}^{⊥} {\bar{u}}_{n} (P_{n}, Y)] < \infty . \end{matrix}$
Proof.
The equation for $b_{n, 1}^{⊥}$ is given by equation (4.47). We thus obtain a bound by using the spectrum of $L_{0}^{n}$ . We can write $\begin{matrix} (4.50) & {(L_{0}^{n})}^{- 1} (Φ) = \sum_{m \neq n} \frac{{⟨ u_{m} (P_{m}, \cdot), Φ (\cdot, Y, q, p) ⟩}_{L^{2} (Γ)} u_{m} (P_{m}, X)}{E_{n} (P_{n}) - E_{m} (P_{m})} . \end{matrix}$ Let $g = {min}_{ξ \in {[- π, π]}^{d}} {| E_{n} (ξ) - E_{n - 1} (ξ) |, | E_{n} (ξ) - E_{n + 1} (ξ) |}$ . Then for each $k \in N$ , we obtain $\begin{array}{l} M_{k} [b_{n, 1}^{⊥} {\bar{u}}_{n} (P_{n}, Y)] \\ ⩽ M_{k} [\frac{1}{g} \sum_{m \neq n} {⟨ u_{m} (P_{m}, \cdot), b_{n, 0} (t, \cdot, q, p) {\bar{u}}_{n} (p, Y) ⟩}_{L^{2} (Γ)} u_{n} (P_{n}, X)] \\ (4.51) & = M_{k} [\frac{{\bar{u}}_{n} (p, Y)}{g} \sum_{m \neq n} a_{n, 0} (t, q, p) {⟨ u_{m} (P_{m}, \cdot), u_{n} (P_{n}, \cdot) ⟩}_{L^{2} (Γ)} u_{n} (P_{n}, X)] . \end{array}$ Hence, by Propositions 4.4 and 4.7, it suffices to control $\begin{matrix} (4.52) & M_{k} [\sum_{m \neq n} {⟨ u_{m} (P_{m}, X), u_{n} (P_{n}, X) ⟩}_{L^{2} (Γ)}] . \end{matrix}$ Since $\int_{Γ} | u_{n} (ξ, x) |^{2} d x = 1$ , Bessel’s inequality implies that the above is finite. □

4.4. Proof of Theorem 3.1

We will need the following estimate, which is proved in [11, Lemma 2.8].

Lemma 4.9.

Suppose $H (ε)$ is a family of self-adjoint operators for $ε > 0$ . Suppose $ψ (t, ε)$ belongs to the domain of $H (ε)$ , is continuously differentiable in t and approximately solves the Schrodinger equation in the sense that $\begin{matrix} (4.53) & i ε \frac{\partial ψ}{\partial t} (t, ε) = H (ε) ψ (t, ε) + ζ (t, ε), \end{matrix}$ where $ζ (t, ε)$ satisfies $\begin{matrix} (4.54) & ‖ ζ (t, ε) ‖ ⩽ μ (t, ε) . \end{matrix}$ Then, $\begin{matrix} (4.55) & e^{- i t H (ε) / ε} ψ (0, ε) - ψ (t, ε) ⩽ ε^{- 1} \int_{0}^{t} μ (s, ε) d s . \end{matrix}$

Moreover, for the Fourier integral operator, we have

Lemma 4.10.

If, for fixed $x, y \in R^{d}$ , $u (x, y, q, p) \in L^{\infty} (Ω; C)$ , for each $n \in N$ and any t, $I^{ε} (u)$ can be extended to a linear bounded operator on $L^{2} (R^{d}, C)$ , and we have $\begin{matrix} (4.56) & {‖ I^{ε} (u) ‖}_{L (L^{2} (R^{d}; C))} ⩽ 2^{- d / 2} ‖ u ‖_{L^{\infty} (R^{2 d}; C)} . \end{matrix}$

Proof.

The proof of lemma 4.10 is essentially the same as Proposition 3.8 in [25] and thus is omitted here. □

We are now ready to prove Theorem 3.1.

Proof of Theorem 3.1.

Computing $i ε \frac{\partial}{\partial t} + \frac{1}{2} ε^{2} \nabla^{2} - V (X) - U (x)$ applied to $I^{ε} (b_{n}^{ε, 1} (t, \frac{x}{ε}, q, p) {\bar{u}}_{n} (p, \frac{y}{ε}))$ , we obtain $\begin{matrix} (4.57) & (i ε \frac{d}{d t} + \frac{1}{2} ε^{2} \nabla^{2} - V (X) - U (x)) I^{ε} (b_{n}^{ε, 1} {\bar{u}}_{n} (p, Y)) = I^{ε} (\sum_{j = 0}^{1} ε^{j} v_{n, j}) + ε^{2} I^{ε} (v_{n, 2}^{ε}) . \end{matrix}$ The expressions for $v_{n, 0}$ , $v_{n, 1}$ , and $v_{n, 2}$ follows from (4.32) by expanding $b^{ε}$ and the linearity of $L_{0}^{n}$ , $L_{1}^{n}$ , and $L_{2}^{n}$ . By equations (4.33) and (4.34), $v_{n, 0}$ and $v_{n, 1}$ vanish. The remaining term $\begin{matrix} (4.58) & v_{n, 2}^{ε} = L_{2}^{n} (b_{n}^{ε, 1} {\bar{u}}_{n} (p, Y)) + R (x, q, p) b_{n}^{ε, 1} {\bar{u}}_{n} (p, Y) \end{matrix}$ satisfies $M_{k} [v_{n, 2}^{ε}] < \infty$ by Propositions 4.3, 4.7, and 4.8. Finally, applying Lemma 4.10 and Lemma 4.9 we obtain the inequality in Theorem 3.1. □

Footnotes

Acknowledgements

R.D. and X.Y. were partially supported by the NSF grants DMS-1418936 and DMS-1107291: NSF Research Network in Mathematical Sciences “Kinetic description of emerging challenges in multiscale problems of natural science”. The work of J.L. was supported in part by the Alfred P. Sloan foundation and the National Science Foundation under award DMS-1312659. X.Y. was also partially supported by the Regents Junior Faculty Fellowship and Hellman Family Foundation Faculty Fellowship of University of California, Santa Barbara.

References

Bao ,

Jin and

P.A.

Markowich , On time-splitting spectral approximations for the Schrödinger equation in the semiclassical regime, J. Comput. Phys. 175 (2002), 487–524. doi:10.1006/jcph.2001.6956.

Bechouche ,

N.J.

Mauser and

Poupaud , Semiclassical limit for the Schrödinger–Poisson equation in a crystal, Comm. Pure Appl. Math. 54 (2001), 851–890. doi:10.1002/cpa.3004.

Bensoussan ,

J.L.

Lions and

Papanicolaou , Asymptotic Analysis for Periodic Structures, North-Holland Publishing Co., Amsterdam-New York, 1978.

M.V.

Berry , Quantal phase factors accompanying adiabatic changes, Proc. Roy. Soc. Lond. A 392 (1984), 45–57. doi:10.1098/rspa.1984.0023.

Carles and

Sparber , Semiclassical wave packet dynamics in Schrödinger equations with periodic potentials, Discrete Contin. Dyn. Syst. Ser. B 12 (2012), 759–774. doi:10.3934/dcdsb.2012.17.759.

Delgadillo ,

Lu and

Yang , Gauge-invariant frozen Gaussian approximation method for the Schrödinger equation with periodic potentials, SIAM J. Sci. Comput. 38 (2016), A2440–A2463. doi:10.1137/15M1040384.

Dimassi ,

J.C.

Guillot and

Ralston , Gaussian beam construction for adiabatic pertubations, Math. Phys. Anal. Geom. 9 (2006), 187–201. doi:10.1007/s11040-006-9009-9.

E ,

Lu and

Yang , Asymptotic analysis of the quantum dynamics in crystals: The Bloch–Wigner transform, Bloch dynamics and Berry phase, Acta Math. Appl. Sin. Engl. Ser. 29 (2013), 465–476. doi:10.1007/s10255-011-0095-5.

Faou ,

Gradinaru and

Lubich , Computing semiclassical quantum dynamics with Hagedorn wavepackets, Commun. Math. Phys. 31 (2009), 3027–3041.

10.

G.A.

Hagedorn , Semiclassical quantum mechanics I:

ℏ \to 0

limit for coherent states, Commun. Math. Phys. 71 (1980), 77–93. doi:10.1007/BF01230088.

11.

G.A.

Hagedorn , Raising and lowering operators for semiclassical wave packets, Annals of Physics 269 (2002), 77–104. doi:10.1006/aphy.1998.5843.

12.

M.F.

Herman and

Kluk , A semiclassical justification for the use of non-spreading wavepackets in dynamics calculations, Chem. Phys. 91 (1984), 27–34. doi:10.1016/0301-0104(84)80039-7.

13.

Huang ,

Jin ,

P.A.

Markowich and

Sparber , A Bloch decomposition-based split-step pseudospectral method for quantum dynamics with periodic potentials, SIAM J. Sci. Comput. 29 (2007), 515–538. doi:10.1137/060652026.

14.

Huang ,

Jin ,

P.A.

Markowich and

Sparber , Numerical simulation of the nonlinear Schrödinger equation with multi-dimensional periodic potentials, Multiscale Model. Simul. 7 (2008), 539–564. doi:10.1137/070699433.

15.

Huang ,

Jin ,

P.A.

Markowich and

Sparber , On the Bloch decomposition based spectral method for wave propagation in periodic media, Wave Motion 46 (2009), 15–28. doi:10.1016/j.wavemoti.2008.07.001.

16.

Jefferis and

Jin , A Gaussian beam method for high frequency solution of symmetric hyperbolic systems with polarized waves, 2014, preprint.

17.

Jin ,

P.A.

Markowich and

Sparber , Mathematical and computational methods for semiclassical Schrödinger equations, Acta Numer. 20 (2011), 211–289. doi:10.1017/S0962492911000031.

18.

Jin ,

Wu and

Yang , Gaussian beam methods for the Schrödinger equation in the semi-classical regime: Lagrangian and Eulerian formulations, Commun. Math. Sci. 6 (2008), 995–1020. doi:10.4310/CMS.2008.v6.n4.a10.

19.

Jin ,

Wu and

Yang , A numerical study of the Gaussian beam methods for one-dimensional Schrödinger–Poisson equations, J. Comput. Math. 28 (2010), 261–272.

20.

Jin ,

Wu and

Yang , Semi-Eulerian and high order Gaussian beam methods for the Schrödinger equation in the semiclassical regime, Commun. Comput. Phys. 9 (2011), 668–687. doi:10.4208/cicp.091009.160310s.

21.

Jin ,

Wu ,

Yang and

Huang , Bloch decomposition-based Gaussian beam method for the Schrödinger equation with periodic potentials, J. Comput. Phys. 229 (2010), 4869–4883. doi:10.1016/j.jcp.2010.01.025.

22.

Kay , Integral expressions for the semi-classical time-dependent propagator, J. Chem. Phys. 100 (1994), 4377–4392. doi:10.1063/1.466320.

23.

Kay , The Herman–Kluk approximation: Derivation and semiclassical corrections, Chem. Phys. 322 (2006), 3–12. doi:10.1016/j.chemphys.2005.06.019.

24.

Lu and

Yang , Frozen Gaussian approximation for high frequency wave propagation, Commun. Math. Sci. 9 (2011), 663–683. doi:10.4310/CMS.2011.v9.n3.a2.

25.

Lu and

Yang , Convergence of frozen Gaussian approximation for high frequency wave propagation, Comm. Pure Appl. Math. 65 (2012), 759–789. doi:10.1002/cpa.21384.

26.

Lu and

Yang , Frozen Gaussian approximation for general linear strictly hyperbolic system: Formulation and Eulerian methods, Multiscale Model. Simul. 10 (2012), 451–472. doi:10.1137/10081068X.

27.

Markowich ,

Pietra and

Pohl , Numerical approximation of quadratic observable of Schrödinger equation-type equations in the semiclassical limit, Numer. Math. 81 (1999), 595–630. doi:10.1007/s002110050406.

28.

Markowich ,

Pietra ,

Pohl and

Stimming , A Wigner-measure analysis of the Dufort–Frankel scheme for the Schrödinger equation, SIAM J. Numer. Anal. 40 (2000), 1281–1310. doi:10.1137/S0036142900381734.

29.

Motamed and

Runborg , Taylor expansion and discretization errors in Gaussian beam superposition, Wave Motion 47 (2010), 421–439.

30.

Panati ,

Spohn and

Teufel , Effective dynamics for Bloch electrons: Peierls substitution and beyond, Comm. Math. Phys. 242 (2003), 547–578. doi:10.1007/s00220-003-0950-1.

31.

Panati ,

Spohn and

Teufel , Motions of Electrons in Adiabatically Perturbed Periodic Structures, Springer, Berlin, 2006, pp. 595–617.

32.

M.M.

Popov , A new method of computation of wave fields using Gaussian beams, Wave Motion 4 (1982), 85–97. doi:10.1016/0165-2125(82)90016-6.

33.

Qian and

Ying , Fast multiscale Gaussian wavepacket transforms and multiscale Gaussian beams for the wave equation, Multiscale Model. Simul. 8 (2010), 1803–1837. doi:10.1137/100787313.

34.

Ralston , Gaussian beams and the propagation of singularities, Studies in PDEs, MAA Stud. Math. 23 (1982), 206–248.

35.

Reed and

Simon , Methods of Modern Mathematical Physics, Vol IV, Academic Press, New York, 1980.

36.

Robert , On the Herman–Kluk semiclassical approximation, Rev. Math. Phys. 22 (2010), 1123–1145. doi:10.1142/S0129055X1000417X.

37.

Shapere and

Wilczek , Geometric Phases in Physics, World Scientific, Singapore, 1989. doi:10.1142/0613.

38.

Simon , Holonomy, the quantum adiabatic theorem, and Berry’s phase, Phys. Rev. Lett. 51 (1983), 2167. doi:10.1103/PhysRevLett.51.2167.

39.

Swart and

Rousse , A mathematical justification of the Herman–Kluk propagator, Commun. Math. Phys. 286 (2009), 725–750. doi:10.1007/s00220-008-0681-4.

40.

N.M.

Tanushev , Superpositions and higher order Gaussian beams, Commun. Math. Sci. 6 (2008), 449–475. doi:10.4310/CMS.2008.v6.n2.a9.

41.

N.M.

Tanushev ,

Qian and

Ralston , Mountain waves and Gaussian beams, Multiscale Model. Simul. 6 (2007), 688–709. doi:10.1137/060673667.

42.

Xiao ,

M.-C.

Chang and

Niu , Berry phase effects on electronic properties, Rev. Mod. Phys. 82 (2010), 1959. doi:10.1103/RevModPhys.82.1959.

Frozen Gaussian approximation for high frequency wave propagation in periodic media

Abstract

Keywords

1. Introduction

4.1. Initial condition

4.2. Estimates of the Hamiltonian flows

Definition 4.1 (Canonical Transformation).

4.3.2. Next order term b n , 1

Footnotes

Acknowledgements

References

4.3.2. Next order term $b_{n, 1}$