Semiclassical wave packet dynamics in Schrodinger equations with periodic potentials

We consider semiclassically scaled Schrodinger equations with an external potential and a highly oscillatory periodic potential. We construct asymptotic solutions in the form of semiclassical wave packets. These solutions are concentrated (both, in space and in frequency) around the effective semiclassical phase-space flow obtained by Peierl's substitution, and involve a slowly varying envelope whose dynamics is governed by a homogenized Schrodinger equation with time-dependent effective mass. The corresponding adiabatic decoupling of the slow and fast degrees of freedom is shown to be valid up to Ehrenfest time scales.


Introduction
1.1. General setting. We consider the following semiclassically scaled Schrödinger equation: ψ ε |t=0 = ψ ε 0 , with d 1, the spatial dimension, and ψ ε = ψ ε (t, x) ∈ C. Here, we already have rescaled all physical parameters such that only one semiclassical parameter ε > 0 (i.e. the scaled Planck's constant) remains. In the following we shall be interested in the asymptotic description of ψ ε (t, x) for ε ≪ 1. To this end, the potential V Γ (y) ∈ R is assumed to be smooth and periodic with respect to some regular lattice Γ ≃ Z d , generated by a given basis {η 1 , . . . , η d }, η ℓ ∈ R d , i.e.
In addition, the slowly-varying potential V is assumed to satisfy the following: Note that this implies that V (x) grows at most quadratically at infinity. Equation 1.1 describes the dynamics of quantum particles in a periodic lattice-potential V Γ under the influence of an external, slowly varying driving force F = −∇V (x). A typical application arises in solid state physics where (1.1) describes the timeevolution of electrons moving in a crystalline lattice (generated by the ionic cores). The asymptotics of (1.1) as ε → 0 + is a natural two-scale problem which is wellstudied in the physics and mathematics literature. Early mathematical results are based on time-dependent WKB type expansions [3,15,37] (see also [7] for a more recent application in the nonlinear case), which, however, suffer from the appearance of caustics and are thus only valid for small times. In order to overcome this problem, other methods based on, e.g., Gaussian beams [11], or Wigner measures [13,14], have been developed. These approaches yield an asymptotic description for time-scales of order O(1) (i.e. beyond caustics). More recently, so-called spaceadiabatic perturbation theory has been used (together with Weyl pseudo-differential calculus) to derive an effective Hamiltonian, governing the dynamics of particles in periodic potentials V Γ under the additional influence of slowly varying perturbations [22,39]. The semi-classical asymptotics of this effective model is then obtained in a second step, invoking an Egorov-type theorem.
On the other hand, it is well known that in the case without periodic potential, semiclassical approximations which are valid up to Ehrenfest time t ∼ O(ln 1/ε) can be constructed in a rather simple way. The corresponding asymptotic method is based on propagating semiclassical wave packets, or coherent states, i.e. approximate solutions of (1.1) which are sufficiently concentrated in space and in frequency around the classical Hamiltonian phase-space flow. More precisely, one considers and the purely time-dependent function S(t) denotes the classical action (see §1.3 below). The right hand side of (1.3) corresponds to a wave function which is equally localized in space and in frequency (at scale √ ε), so the uncertainty principle is optimized. In other words, the three quantities , and have the same order of magnitude, O(1), as ε → 0. The basic idea for this type of asymptotic method can be found in the classical works of [16,27] (see also [4,29] for a broader introduction). It has been developed further in, e.g., [8,9,34,35,38] and in addition also proved to be applicable in the case of nonlinear Schrödinger equations [6] (a situation in which the use of Wigner measures of spaceadiabatic perturbation theory fails). Asymptotic results based on such semiclassical wave packets also have the advantage of giving a rather clear connection between quantum mechanics and classical particle dynamics and are thus frequently used in numerical simulations (see e.g. [12]). Ehrenfest time is the largest time up to which the wave packet approximation is valid, in general. Without any extra geometric assumption, the coherent structure may be lost at some time of order C ln 1/ε, if C is too large. See e.g. [5,10,30,32,33] and references therein.
Interestingly enough, though, it seems that so far this method has not been extended to include also highly oscillatory periodic potentials V Γ x ε , and it will be the main task of this work to do so. To this end, it will be necessary to understand the influence of V Γ x ε on the dispersive properties of the solution ψ ε (t, x). In particular, having in mind the results quoted above, one expects that in this case the usual kinetic energy of a particle E = 1 2 |k| 2 has to be replaced by E m (k), i.e. the energy of the m-th Bloch band associated to V Γ . In physics this is known under the name Peierls substitution. We shall show that under the additional influence of a slowly varying potential V (x), this procedure is in fact asymptotically correct (i.e. for ε ≪ 1) up to Ehrenfest time, provided the initial data ψ ε 0 is sufficiently concentrated around (q 0 , p 0 ) ∈ R 2d . Remark 1.2. Indeed, we could also allow for time-dependent external potentials V (t, x) ∈ R measurable in time, smooth in x, and satisfying . Under this assumptions, it is straightforward to adapt the analysis given below. For the sake of notation, we shall not do so here, but rather leave the details to the reader.
1.2. Bloch and semiclassical wave packets. In order to state our result more precisely, we first recall some well-known results on the spectral theory for periodic Schrödinger operators, cf. [31,40]: Denote by Y ⊂ Γ the (centered) fundamental domain of the lattice Γ, equipped with periodic boundary conditions, i.e. Y ≃ T d . Similarly, we denote by Y * ≃ T d the fundamental domain of the corresponding dual lattice. The latter is usually referred to as the Brillouin zone. Bloch-Floquet theory asserts that H per admits a fiber-decomposition where for k ∈ Y * , we denote It therefore suffices to consider the following spectral problem on Y : where E m (k) ∈ R and χ m (y, k), respectively, denote an eigenvalue/eigenvector pair of H Γ (k), parametrized by k ∈ Y * , the so-called crystal momentum. These eigenvalues can be ordered increasingly, such that where each eigenvalue is repeated according to its multiplicity (which is known to be finite). This implies that where {E m (k); k ∈ Y * } is called the m-th energy band (or Bloch band). The associated eigenfunctions χ m (y, k) are Γ * -periodic w.r.t. k and form a complete orthonormal basis of L 2 (Y ). Moreover, the functions χ m (·, k) ∈ H 2 (Y ) are known to be real-analytic with respect to k on Y * \Ω, where Ω is a set of Lebesgue measure zero (the set of band crossings). Next, we consider for some m ∈ N the corresponding semi-classical band Hamiltonian, obtained by Peierl's substitution, i.e.
and denote the semiclassical phase space trajectories associated to h sc m by (1.6) q(t) = ∇ k E m (p(t)) , q(0) = q 0 , This system is the analog of (1.4) in the presence of an additional periodic potential.
that is, a shift with constant speed ω = ∇E m (p 0 ).
In order to make sure that the system (1.6) is well-defined, we shall from now on impose the following condition on E m (k). Assumption 1.4. We assume that E m (p(t)) is a simple eigenvalue, uniformly for all t ∈ R, i.e. there exists a δ > 0, such that It is known that if E m (k) is simple, it is infinitely differentiable and thus the right hand side of (1.6) is well defined. Under Assumption 1.4, we consequently obtain a smooth semi-classical flow (q 0 , p 0 ) → (q(t), p(t)), for all t ∈ R. In addition, one can choose χ m (y, k) to be Γ-periodic with respect to y and such that (y, t) → χ m (y, p(t)) is bounded together with all its derivatives. Example 1.5. By compactness of Y * , Assumption 1.4 is satisfied in either of the following two cases:

Main result.
With the above definitions at hand, we are now able to state our main mathematical result. To this end, we first define a semiclassical wave packet in the m-th Bloch band (satisfying Assumption 1.4) by with q(t), p(t) given by system (1.6) and u(t, z) ∈ C, a smooth slowly varying envelope which will determined by an envelope equation yet to be derived (see below). In addition, the ε-oscillatory phase is where S m (t) ∈ R is the (purely time-dependent) semi-classical action with L m denoting the Lagrangian associated to the effective Hamiltonian h sc m , i.e.
Remark 1.6. Note that this is nothing but the Legendre transform of the effective Hamiltonian h sc m . As in classical mechanics, one associates to a given Hamiltonian H(p, q) a Lagrangian via L(p, q) = p ·q − H(p, q).
The function ϕ ε given by (1.8) generalizes the usual class of semiclassical wave packets considered in e.g. [16,27]. Note that in contrast to two-scale WKB approximation considered in [3,15,37], it involves an additional scale of the order 1/ √ ε, the scale of concentration of the amplitude u. In addition, (1.8) does not suffer from the appearance of caustics. Nevertheless, in comparison to the highly oscillatory Bloch function χ m , the amplitude is still slowly varying and thus we can expect an adiabatic decoupling between the slow and fast scales to hold on (long) macroscopic time-scales. Indeed, we shall prove the following result: Theorem 1.7. Let V Γ be smooth and V satisfy Assumptions 1.1. In addition, let Assumption 1.4 hold and the initial data be given by with q 0 , p 0 ∈ R d and some given profile u 0 ∈ S(R d ).
Then there exists C > 0 such that the solution of (1.1) can be approximated by Here, ϕ ε is given by (1.8) with where β(t) ∈ iR is the so-called Berry phase term and v ∈ C(R; S(R d )) satisfies the following homogenized Schrödinger equation In particular there exists C 0 > 0 so that Remark 1.8. In fact it is possible to prove the same result under less restrictive regularity assumptions on u 0 and V Γ . Indeed, Proposition 5.1 shows that it is sufficient to require that u 0 belongs to a certain weighted Sobolev space. Concerning the periodic potential, it is possible to lower the regularity considerably, depending on the dimension. For example, in d = 3 it is sufficient to assume V Γ to be infinitesimally bounded with respect to −∆. This implies which, together with several density arguments (to be invoked at different stages of the formal expansion), is enough to justify the analysis given below. Theorem 1.7 provides an approximate description of the solution to (1.1) up to Ehrenfest time and can be seen as the analog of the results given in [16,27,8,9,18,34,35,38] where the case of slowly varying potentials V (x) is considered. The proof does not rely on the use of pseudo-differential calculus or space spaceadiabatic perturbation theory and can thus be considered to be considerably simpler from a mathematical point of view. In fact, our approach is similar to the one given in [18], which derives an analogous result for the so-called Born-Oppenheimer approximation of molecular dynamics. Note however, that we allow for more general initial amplitudes, not necessarily Gaussian. Indeed, in the special case where the initial envelope u 0 is a Gaussian, then its evolution u remains Gaussian, and can be completely characterized; see §4. 3. Also note that in contrast to the closely related method of Gaussian beams presented in, e.g., [11], we do not need to include complex-valued phases and in addition, obtain an approximation valid for longer times.
The Berry phase term is an example for so-called geometric phases in quantum mechanics. It is a well known feature of semiclassical approximation in periodic potentials, see, e.g., [28] for more details and a geometric interpretation. The homogenized Schrödinger equation features a rather unusual dispersive behavior described by a time-dependent effective mass tensor Hessian of E m (k) evaluated at k = p(t). To our knowledge, Theorem 1.7 is the first result in which a Schrödinger type equation with time-dependent effective mass has been rigorously derived (see also the discussion in Remark 3.1).
Remark 1.9. Let us also mention that the same class of initial data has been considered in [1] for a Schrödinger equation with locally periodic potential V Γ (x, y) and corresponding x-dependent Bloch bands E m (k; x). In this work, the authors derive a homogenized Schrödinger equation, provided that ψ ε 0 is concentrated around a stationary point point x 0 , p 0 of the semiclassical phase space flow, i.e.
This implies q(t) = q 0 and p(t) = p 0 , for all t ∈ R, yielding (at least asymptotically) a localization of the wave function. We observe the same phenomenon in our case under the condition V (x) = 0 and ∇ k E m (k) = 0 (see Example 1.3).
This work is now organized as follows: In the next section, we shall formally derive an approximate solution to (1.1) by means of a (formal) multi-scale expansion. This expansion yields a system of three linear equations, which we shall solve in Section 3. In particular, we shall obtain from it the homogenized Schrödinger equation. The corresponding Cauchy problem is then analyzed in Section 4, where we also include a brief discussion on the particularly important case of Gaussian profiles (yielding a direct connection to [16]). A rigorous stability result for our approximation, up to Ehrenfest time, is then given in Section 5.
Remark 1.10. We expect that our results can be generalized to the case of (weakly) nonlinear Schrödinger equations (as considered in [6,7]). This will be the aim of a future work.

Formal derivation of an approximate solution
2.1. Reduction through exact computations. We seek the solution ψ ε of (1.1) in the following form where the phase φ(t, x) is given by (1.9), the function U ε = U ε (t, z, y) is assumed to be smooth, Γ-periodic with respect to y, and admits an asymptotic expansion Note that due to the inclusion of the factor ε −d/4 the L 2 (R d ) norm of the right hand side of (2.1) is in fact uniformly bounded with respect to ε, whereas the L ∞ (R d ) norm in general will grow as ε → 0. The asymptotic expansion 2.2 therefore has to be understood in the L 2 sense. Taking into account that in view of (1.9), ∇ x φ m (t, x) = p(t), we compute: where in all of the above expressions, the various functions have to be understood to be evaluated as follows: Thus, ordering equal powers of ε in equation (1.9) we find that So far, we have neither used the fact that q(t), p(t) are given by the Hamiltonian flow (1.6), nor the explicit dependence of φ m on time. Using these properties, allows Now, recall that in the above lines, U ε is evaluated at the shifted spatial variable z = (x − q(t))/ √ ε. Taking this into account, we notice that the above hierarchy has to be modified, and we find: Next, we perform a Taylor expansion of V around the point q(t): since V is at most quadratic in view of Assumption 1.1. Recalling that h sc m (p, q) = E m (p)+ V (q), the terms involving V (q) cancel out in b ε 0 , the terms involving ∇V (q) cancel out in b ε 1 , and thus, we finally obtain: Lemma 2.1. Let the Assumptions 1.1, 1.4 hold and ψ ε be related to U ε through (2.1). Then it holds and a remainder r ε (t, z, y) satisfying where the constant C > 0 is independent of t, z, y and ε.

2.2.
Introducing the approximate solution. We now expand U ε in powers of ε, according to (2.2). To this end, we introduce the following (time-dependent) linear operators In order to solve (1.1) up to a sufficiently small error term (in L 2 ), we need to cancel the first three terms in our asymptotic expansion. This yields, the following system of equations Assuming for the moment that we can do so, this means that we (formally) solve (1.1) up to errors of order ε 3/2 (in L 2 ), which is expected to generate a small perturbation of the exact solution (in view of the ε in front of the time derivative of ψ ε in (1.1)). We consequently define the approximate solution where the remainder terms r ε 1 , r ε 2 are given by r ε 1 (t, z, y) = L 2 U 1 (t, z, y), r ε 2 (t, z, y) = L 1 U 2 (t, z, y), and r ε satisfies where the constant C > 0 is independent of t, z, y and ε.

Derivation of the homogenized equation
3.1. Some useful algebraic identities. Given the form of L 0 , the equation Before studying the other two equations, we shall recall some algebraic formulas related to the eigenvalues and eigenvectors of H Γ . First, in view of the identity (1.5), we have Taking the scalar product in L 2 (Y ) with χ m , we infer Since H Γ is self-adjoint, the last term is zero, thanks to (1.5). We infer Differentiating (3.2) again, we have, for all j, ℓ ∈ {1, . . . , d}: Taking the scalar product with χ m , we have:

3.2.
Higher order solvability conditions. By Fredholm's alternative, a necessary and sufficient condition to solve the equation L 0 U 1 + L 1 U 0 = 0, is that L 1 U 0 is orthogonal to ker L 0 , that is: Given the expression of L 1 and the formula (3.1), we compute In view of (3.3), we infer that (3.5) is automatically fulfilled. We thus obtain 1 , the part of U 1 which is orthogonal to ker L 0 , is obtained by inverting an elliptic equation: Note that the formula for L 1 U 0 can also be written as ) χ m (y, p(t)) · ∇ z u(t, z), thus taking into account (3.2), we simply have: u ⊥ 1 (t, z, y) = −i∇ k χ m (y, p(t)) · ∇ z u(t, z). At this stage, we shall, for simplicity choose u 1 = 0, in which case U 1 becomes simply a function of u: (3.6) U 1 (t, z, y) = −i∇ k χ m (y, p(t)) · ∇ z u(t, z).
As a next step in the formal analysis, we must solve By the same argument as before, we require With the expression (3.6), we compute and we also have + iuṗ(t) · ∇ k χ m (y, p(t)) .
Recalling thatṗ(t) = −∇V (q(t)), we find: By making the last sum symmetric with respect to j and ℓ, and using (3.4), we finally obtain the homogenized Schrödinger equation with time-dependent effective mass: where we recall that the so-called Berry phase term. From χ m L 2 (Y ) = 1, we infer that In other words, χ m , ∇ k χ m L 2 (Y ) ∈ iR and thus iβ(t) ∈ R, acts like a purely timedependent, real-valued, potential. Thus, invoking the unitary change of variable implies that v(t, z) solves (1.12). Equation (3.8) models a quantum mechanical time-dependent harmonic oscillator, in which the time dependence is present both in the differential operator, and in the potential.
This equation has been derived in [2] (see also [36,20,23] for similar results). Note, however, that in the quoted works the scaling of the original equation (1.1) is different (i.e. not in semiclassical form).

The envelope equation
We examine the Cauchy problem for (3.8), with special emphasis on the large time control of u.

4.1.
The general Cauchy problem. Equation 3.8 can be seen as the quantum mechanical evolutionary problem corresponding to the following time-dependent Hamiltonian, Under Assumptions 1.1 and 1.4, this Hamiltonian is self-adjoint, smooth in time, and quadratic in (z, ζ) (in fact, at most quadratic would be sufficient). Using the result given in [24, p.197] (see also [25]), we directly infer the following existence result: Lemma 4.1 (From [24]). For d 1 and v 0 ∈ L 2 (R d ), consider the equation If the coefficients a jk and b jk are continuous and real-valued, such that the matrices (a jk ) j,k and (b jk ) j,k are symmetric for all time, then (4.2) has a unique solution v ∈ C(R; L 2 (R d )). It satisfies Moreover, if v 0 ∈ Σ k for some k ∈ N, then v ∈ C(R; Σ k ).
In particular, this implies that if u 0 ∈ Σ k , then (1.12) has a unique solution v ∈ C(R; Σ k ). As a consequence, (3.8) has a unique solution u ∈ C(R; Σ k ) such that u |t=0 = u 0 . Remark 4.2. It may happen that the functions a jk are zero on some non-negligible set. In this case, (4.2) ceases to be dispersive. Note that the standard harmonic oscillator is dispersive, locally in time only, since it has eigenvalues. We shall see that this is not a problem in our analysis though.

4.2.
Exponential control of the envelope equation. To prove Theorem 1.7, we need to control the error present in Lemma 2.2 for large time. In general, i.e. without extra geometric assumptions on the wave packet, exponential growth in time must be expected: Then the solution u to (3.8) satisfies u ∈ C(R; Σ k ), and there exists C > 0 such that Proof. The result can be established by induction on k. The constant C must actually be expected to depend on k, as shown by the case of There, the fundamental solution is explicit (generalized Mehler formula, see e.g. [21]), and we check that u(t) Σ k behaves like e kt . For k = 0, the result is obvious, since in view of Lemma 4.1, the L 2 -norm is conserved. The case k = 1 illustrates the general mechanism of the proof, and we shall stick to this case for simplicity. The key remark is that even though the operators z and ∇ z (involved in the definition of Σ 1 ) do not commute with the Hamiltonian (4.1), the commutators yield a closed system of estimates. First, multiplying (3.8) by z, we find and, in view of Assumption 1.1, Summing over the two inequalities and using the conservation of mass, we infer and Gronwall's lemma yields the proposition in the case k = 1. By induction, applying (z, ∇ z ) to (3.8) k times, the defects of commutation always yield the same sort of estimate, and the proposition follows easily.

Gaussian wave packets.
In the case where the initial datum in (3.8) is a Gaussian, we can compute its evolution and show that it remains Gaussian, by following the same strategy as in [16] (see also [17,18]). As a matter of fact, the order in which we have proceeded is different from the one in [16], since we have isolated the envelope equation (3.8) before considering special initial data. As a consequence, we have fewer unknowns. Consider (3.8) with initial datum where the matrices A and B satisfy the following properties: A and B are invertible;   where the matrices A(t) and B(t) evolve according to the differential equations (4.10) In addition, for all time t ∈ R, A(t) and B(t) satisfy (4.5)-(4.8).
Proof. The argument is the same as in [16] (see also [17,18]): One easily checks that if A(t) and B(t) evolve according to (4.10), then u given by (4.9) solves (3.8). On the other hand, it is clear that (4.10) has a global solution. Finally, since ∇ 2 k E m (p(t)) and ∇ 2 x V (q(t)) are symmetric matrices, it follows from [16, Lemma 2.1] that for all time, A(t) and B(t) satisfy (4.5)-(4.8).

Stability of the approximation up to Ehrenfest time
As a final step we need to show that the derived approximation ψ ε app (t) indeed approximates the exact solution ψ ε (t) up to Ehrenfest time.
The a-priori L 2 estimate yields The assertion then follows from Lemma 2.2 (establishing the needed properties for the functions r ε , r ε 1 and r ε 2 ), Proposition 3.2, and Proposition 4.3. With this approach, we need to know that r ε is in L 2 z , so U 0 , U 1 and U 2 have three momenta in L 2 z : in view of Proposition 3.2 and Proposition 4.3, this amounts to demanding u 0 ∈ Σ 5 . This asymptotic stability result directly yields the assertion of Theorem 1.7.
Remark 5.2. The construction of the approximate solution ψ ε app has forced us to introduce non-zero correctors U 1 and U 2 , given by elliptic inversion. Therefore, we had to consider well-prepared initial data for ψ ε app . This aspect is harmless as long as one is interested only in the leading order behavior of ψ ε as ε → 0. As a consequence, our approach would not allow us to construct arbitrary accurate approximations for ψ ε (in terms of powers of ε), unless well-prepared initial data are considered, i.e. data lying in so-called super-adiabatic subspaces, in the terminology of [28] (after [26]). This is due to the spectral analysis implied by the presence of the periodic potential V Γ , and shows a sharp contrast with the case V Γ = 0.
Of course the above given stability result immediately generalizes to situations where, instead of a single ϕ ε , a superposition of finitely many semiclassical wave packets is considered, u n x − q n √ ε χ mn x ε , p n e ipn·(x−qn)/ε .
Since the underlying semiclassical Schrödinger equation (1.1) is linear, each of these initial wave packets will evolve individually from the rest, as in Theorem 1.7. Up to some technical modifications, it should be possible to consider even a continuous superposition of wave packets, yielding a semiclassical approximation known under the name "frozen Gaussians", see [19].