Mathematical and computational methods for semiclassical Schrödinger equations*

We consider time-dependent (linear and nonlinear) Schrödinger equations in a semiclassical scaling. These equations form a canonical class of (nonlinear) dispersive models whose solutions exhibit high-frequency oscillations. The design of efficient numerical methods which produce an accurate approximation of the solutions, or at least of the associated physical observables, is a formidable mathematical challenge. In this article we shall review the basic analytical methods for dealing with such equations, including WKB asymptotics, Wigner measure techniques and Gaussian beams. Moreover, we shall give an overview of the current state of the art of numerical methods (most of which are based on the described analytical techniques) for the Schrödinger equation in the semiclassical regime.

We consider time-dependent (linear and nonlinear) Schrödinger equations in a semiclassical scaling. These equations form a canonical class of (nonlinear) dispersive models whose solutions exhibit high-frequency oscillations. The design of efficient numerical methods which produce an accurate approximation of the solutions, or at least of the associated physical observables, is a formidable mathematical challenge. In this article we shall review the basic analytical methods for dealing with such equations, including WKB asymptotics, Wigner measure techniques and Gaussian beams. Moreover, we shall give an overview of the current state of the art of numerical methods (most of which are based on the described analytical techniques) for the Schrödinger equation in the semiclassical regime.

Introduction
This goal of this article is to give an overview of the currently available numerical methods used in the study of highly oscillatory partial differential equations (PDEs) of Schrödinger type. This type of equation forms a canonical class of (nonlinear) dispersive PDEs, i.e., equations in which waves of different frequency travel with different speed. The accurate and efficient numerical computation of such equations usually requires a lot of analytical insight, and this applies in particular to the regime of high frequencies.
The following equation can be seen as a paradigm for the PDEs under consideration: for (t, x) ∈ R × R d , with d ∈ N denoting the space dimension. In addition, ε ∈ (0, 1] denotes the small semiclassical parameter (the scaled Planck's constant), describing the microscopic/macroscopic scale ratio. Here, we have already rescaled all physical parameters, such that only one dimensionless parameter ε 1 remains. The unknown u ε = u ε (t, x) ∈ C is the quantum mechanical wave function whose dynamics is governed by a static potential function V = V (x) ∈ R (time-dependent potentials V (t, x) can usually also be taken into account without requiring too much extra work, but for the sake of simplicity we shall not do so here). In this article, several different classes of potentials, e.g., smooth, discontinuous, periodic, random, will be discussed, each of which requires a different numerical strategy. In addition, possible nonlinear effects can be taken into account (as we shall do in Section 15) by considering nonlinear potentials V = f (|u ε | 2 ).
In the absence of V (x) a particular solution to the Schrödinger equation is given by a single plane wave, for any given wave vector ξ ∈ R d . We see that u ε features oscillations with frequency 1/ε in space and time, which inhibit strong convergence of the wave function in the classical limit ε → 0 + . In addition, these oscillations pose a huge challenge for numerical computations of (2.1); in particular, they strain computational resources when run-of-the-mill numerical techniques are applied in order to numerically solve (1.1) in the semiclassical regime ε 1. For the linear Schrödinger equation, classical numerical analysis methods (such as the stability-consistency concept) are sufficient to derive meshing strategies for discretizations (say, of finite difference, finite element or even time-splitting spectral type) which guarantee (locally) strong convergence of the discrete wave functions when the semiclassical parameter ε is fixed: see Chan, Lee and Shen (1986), Chan and Shen (1987), Wu (1996) and Dörfler (1998); extensions to nonlinear Schrödinger equations can be found in, e.g., Delfour, Fortin and Payre (1981), Taha and Ablowitz (1984) and Pathria and Morris (1990). However, the classical numerical analysis strategies cannot be employed to investigate uniform in ε properties of discretization schemes in the semiclassical limit regime. As we shall detail in Section 4, even seemingly reasonable, i.e., stable and consistent, discretization schemes, which are heavily used in many practical application areas of Schrödinger-type equations, require huge computational resources in order to give accurate physical observables for ε 1. The situation gets even worse when an accurate resolution of u ε itself is required. To this end, we remark that time-splitting spectral methods tend to behave better than finite difference/finite element methods, as we shall see in more detail in Section 5.
In summary, there is clearly a big risk in using classical discretization techniques for Schrödinger calculations in the semiclassical regime. Certain schemes produce completely incorrect observables under seemingly reasonable meshing strategies, i.e., an asymptotic resolution of the oscillation is not always enough. Even worse, in these cases there is no warning from the scheme (such as destabilization) that something has gone wrong in the computation (since local error control is computationally not feasible in the semiclassical regime). The only safety anchor here lies in asymptotic mathematical analysis, such as WKB analysis, and/or a physical insight into the problem. They typically yield a (rigorous) asymptotic description of u ε for small ε 1, which can consequently be implemented on a numerical level, providing an asymptotic numerical scheme for the problem at hand. In this work, we shall discuss several asymptotic schemes, depending on the particular type of potentials V considered.
While one can not expect to be able to pass to the classical limit directly in the solution u ε of (1.1), one should note that densities of physical observables, which are the quantities most interesting in practical applications, are typically better behaved as ε → 0, since they are quadratic in the wave function (see Section 2.1 below). However, weak convergence of u ε as ε → 0 is not sufficient for passing to the limit in the observable densities (since weak convergence does not commute with nonlinear operations). This makes the analysis of the semiclassical limit a mathematically highly complex issue. Recently, much progress has been made in this area, particularly by using tools from micro-local analysis, such as H-measures (Tartar 1990) and Wigner measures (Lions and Paul 1993, Markowich and Mauser 1993, Gérard, Markowich, Mauser and Poupaud 1997. These techniques go far beyond classical WKB methods, since the latter suffers from the appearance of caustics: see, e.g., Sparber, Markowich and Mauser (2003) for a recent comparison of the two methods. In contrast, Wigner measure techniques reveal a kinetic equation in phase space, whose solution, the so-called Wigner measure associated to u ε , does not exhibit caustics (see Section 3 for more details).
A word of caution is in order. First, a reconstruction of the asymptotic description for u ε itself (for ε 1) is in general not straightforward, since, typically, some phase information is lost when passing to the Wigner picture. Second, phase space techniques have proved to be very powerful in the linear case and in certain weakly nonlinear regimes, but they have not yet shown much strength when applied to nonlinear Schrödinger equations in the regime of supercritical geometric optics (see Section 15.2). There, classical WKB analysis (and in some special cases techniques for fully integrable systems) still prevails. The main mathematical reason for this is that the initial value problem for the linear Schrödinger equation propagates only one ε-scale of oscillations, provided the initial datum in itself is ε-oscillatory (as is always assumed in WKB analysis). New (spatial) frequencies ξ may be generated during the time evolution (typically at caustics) but no new scales of oscillations will arise in the linear case. For nonlinear Schrödinger problems this is different, as new oscillation scales may be generated through the nonlinear interaction of the solution with itself. Further, one should note that this important analytical distinction, i.e., no generation of new scales but possible generation of new frequencies (in the linear case), may not be relevant on the numerical level, since, say, 100ε is analytically just a new frequency, but numerically a new scale.
Aside from semiclassical situations, modern research in the numerical solution of Schrödinger-type equations goes in a variety of directions, of which the most important are as follows.
(i) Stationary problems stemming from, e.g., material sciences. We mention band diagram computations (to be touched upon below in Sections 12 and 13) and density functional theory for approximating the full microscopic Hamiltonian (not to be discussed in this paper). The main difference between stationary and time-dependent semiclassical problems is given by the fact that in the former situation the spatial frequency is fixed, whereas in the latter (as mentioned earlier) new frequencies may arise over the course of time. (ii) Large space dimensions d 1, arising, for example, when the number of particles N 1, since the quantum mechanical Hilbert space for N indistinguishable particles (without spin) is given by L 2 (R 3N ). This is extremely important in quantum chemistry simulations of atomistic/molecular applications. Totally different analytical and numerical techniques need to be used and we shall not elaborate on these issues in this paper. We only remark that if some of the particles are very heavy and can thus be treated classically (invoking the so-called Born-Oppenheimer approximation: see Section 11), a combination of numerical methods for both d 1 and ε 1 has to be used.

Basic existence results and physical observables
We recall the basic existence theory for linear Schrödinger equations of the form For the sake of simplicity we assume the (real-valued) potential V = V (x) to be continuous and bounded, i.e., The Kato-Rellich theorem then ensures that the Hamiltonian operator is essentially self-adjoint on D(−∆) = C ∞ 0 ⊂ L 2 (R d ; C) and bounded from below by −K: see, e.g., Reed and Simon (1975). Its unique self-adjoint extension (to be denoted by the same symbol) therefore generates a strongly continuous semi-group U ε (t) = e − itH ε /ε on L 2 (R d ), which ensures the global existence of a unique (mild) solution u ε (t) = U ε (t)u in of the Schrödinger equation (2.1). Moreover, since U ε (t) is unitary, we have u ε (t, ·) 2 L 2 = u ε in 2 L 2 , ∀ t ∈ R. In quantum mechanics this is interpreted as conservation of mass. In addition, we also have conservation of the total energy which is the sum of the kinetic and the potential energies. In general, expectation values of physical observables are computed via quadratic functionals of u ε . To this end, denote by a W (x, εD x ) the operator corresponding to a classical (phase space) observable a ∈ C ∞ b (R d × R d ), obtained via Weyl quantization, x + y 2 , εξ f (y) e i(x−y)·ξ dξ dy, (2.4) where εD x := −iε∂ x . Then the expectation value of a in the state u ε at time t ∈ R is given by where ·, · L 2 denotes the usual scalar product on L 2 (R d ; C).
Remark 2.1. The convenience of the Weyl calculus lies in the fact that an (essentially) self-adjoint Weyl operator a W (x, εD x ) has a real-valued symbol a(x, ξ): see Hörmander (1985).
The quantum mechanical wave function u ε can therefore be considered only an auxiliary quantity, whereas (real-valued) quadratic quantities of u ε yield probability densities for the respective physical observables. The most basic quadratic quantities are the particle density (2.6) and the current density j ε (t, x) := ε Im u ε (t, x)∇u ε (t, x) . (2.7) It is easily seen that if u ε solves (2.1), then the following conservation law holds: ∂ t ρ ε + div j ε = 0. (2.8) In view of (2.3) we can also define the energy density e ε (t, x) := 1 2 |ε∇u ε (t, x)| 2 + V (x)ρ ε (t, x). (2.9) As will be seen (see Section 5), computing these observable densities numerically is usually less cumbersome than computing the actual wave function u ε accurately. From the analytical point of view, however, we face the problem that the classical limit ε → 0 can only be regarded as a weak limit (in a suitable topology), due to the oscillatory nature of u ε . Quadratic operations defining densities of physical observables do not, in general, commute with weak limits, and hence it remains a challenging task to identify the (weak) limits of certain physical observables, or densities, respectively.

Asymptotic description of high frequencies
In order to gain a better understanding of the oscillatory structure of u ε we invoke the following WKB approximation (see Carles (2008) and the references given therein): with real-valued phase S and (possibly) complex-valued amplitude a ε , satisfying the asymptotic expansion a ε ε→0 ∼ a + εa 1 + ε 2 a 2 + · · · . (2.11) Plugging the ansatz (2.10) into (2.1), one can determine an approximate solution to (2.1) by subsequently solving the equations obtained in each order of ε.
To leading order, i.e., terms of order O(1), one obtains a Hamilton-Jacobi equation for the phase function S: This equation can be solved by the method of characteristics, provided V (x) is sufficiently smooth, say V ∈ C 2 (R d ). The characteristic flow is given by the following Hamiltonian system of ordinary differential equations: (2.13) Remark 2.2. The characteristic trajectories y → x(t, y) obtained via (2.13) are usually interpreted as the rays of geometric optics. The WKB approximation considered here is therefore also regarded as the geometric optics limit of the wave field u ε . By the Cauchy-Lipschitz theorem, this system of ordinary differential equations can be solved at least locally in time, and consequently yields the phase function where y(τ, x) denotes the inversion of the characteristic flow X t : y → x(t, y). This yields a smooth phase function S ∈ C ∞ ([−T, T ] × R d ) up to some time T > 0 but possibly very small. The latter is due to the fact that, in general, characteristics will cross at some finite time |T | < ∞, in which case the flow map X t : R d → R d is no longer one-to-one. The set of points at which X t ceases to be a diffeomorphism is usually called a caustic set. See Figure 2.1 (taken from Gosse, Jin and Li (2003)) for examples of caustic formulation. Ignoring the problem of caustics for a moment, one can proceed with our asymptotic expansion and obtain at order O(ε) the following transport equation for the leading-order amplitude: (2.14) In terms of the leading-order particle density ρ := |a| 2 , this reads which is reminiscent of the conservation law (2.8).
The transport equation (2.14) is again solved by the methods of characteristics (as long as S is smooth, i.e., before caustics) and yields where J t (y) := det ∇ y x(t, u) denotes the Jacobi determinant of the Hamiltonian flow. All higher-order amplitudes a n are then found to be solutions of inhomogeneous transport equations of the form ∂ t a n + ∇S · ∇a n + a n 2 ∆S = ∆a n−1 . (2.17) These equations are consequently solved by the method of characteristics. At least locally in time (before caustics), this yields an approximate solution of WKB type, x) + · · · , including amplitudes (a n ) N n=1 up to some order N ∈ N. It is then straightforward to prove the following stability result. Theorem 2.3. Assume that the initial data of (2.1) are given in WKB form, , and let a in ∈ S(R d ), i.e., smooth and rapidly decaying. Then, for any closed time interval I ⊂ T , before caustic onset, there exists a C > 0, independent of ε ∈ (0, 1], such that The first rigorous result of this type goes back to Lax (1957). Its main drawback is the fact that the WKB solution breaks down at caustics, where S develops singularities. In addition, the leading-order amplitude a blows up in L ∞ (R d ), in view of (2.16) and the fact that lim t→T J t (y) = 0. Of course, these problems are not present in the exact solution u ε but are merely an artifact of the WKB ansatz (2.10). Caustics therefore indicate the appearance of new ε-oscillatory scales within u ε , which are not captured by the simple oscillatory ansatz (2.10).

Beyond caustics
At least locally away from caustics, though, the solution can always be described by a superposition of WKB waves. This can be seen rather easily in the case of free dynamics where V (x) = 0. The corresponding solution of the Schrödinger equation (2.1) with WKB initial data is then explicitly given by (2.20) The representation formula (2.19) comprises an oscillatory integral, whose main contributions stem from stationary phase points at which ∂ y,ξ ϕ(x, t) = 0. In view of (2.20) this yields The corresponding map y → x(t, y) is the characteristic flow of the free Hamilton-Jacobi equation Inverting the relation y → x(t, y) yields the required stationary phase points {y j (t, x)} j∈N ∈ R d for the integral (2.19). Assuming for simplicity that there are only finitely many such points, then with constant phase shifts m j = sgn D 2 S(y j (t, x))) ∈ N (usually referred to as the Keller-Maslov index). The right-hand side of this expression is usually referred to as multiphase WKB approximation. The latter can be interpreted as an asymptotic description of interfering wave trains in u ε .
Remark 2.4. The case of non-vanishing V (x), although similar in spirit, is much more involved in general. In order to determine asymptotic description of u ε beyond caustics, one needs to invoke the theory of Fourier integral operators: see, e.g., Duistermaat (1996). In particular, it is in general very hard to determine the precise form and number of caustics appearing throughout the time evolution of S(t, x). Numerical schemes for capturing caustics have been developed in, e.g., Benamou and Solliec (2000), or Benamou, Lafitte, Sentis and Solliec (2003) and the references therein.

The Wigner-transformed picture of quantum mechanics
Whereas WKB-type methods aim for approximate solutions of u ε , the goal of this section is to directly identify the weak limits of physical observable densities as ε → 0. To this end, one defines the so-called Wigner transform of u ε , as given in Wigner (1932): Plancherel's theorem together with a simple change of variables yields The real-valued Wigner transform w ε ∈ L 2 (R d x × R d ξ ) can be interpreted as a phase space description of the quantum state u ε . In contrast to classical phase space distributions, w ε in general also takes negative values (except for Gaussian wave functions).
Applying this transformation to the Schrödinger equation (2.1), the time- is a pseudo-differential operator, taking into account the influence of V (x). Explicitly, it is given by where the symbol δV ε reads Note that in the free case where V (x) = 0, the Wigner equation becomes the free transport equation of classical kinetic theory. Moreover, if V ∈ C 1 (R d ) we obviously have that in which case the ε → 0 limit of (3.2) formally becomes the classical Liouville equation in phase space: see (3.7) below. The most important feature of the Wigner transform is that it allows for a simple computation of quantum mechanical expectation values of physical observables. Namely, where a(x, ξ) is the classical symbol of the operator a W (x, εD x ). In addition, at least formally (since w ε ∈ L 1 (R d × R d ) in general), the particle density (2.6) can be computed via and the current density (2.7) is given by S. Jin and P. Markowich and C. Sparber Similarly, the energy density (2.9) is where the classical (phase space) Hamiltonian function is denoted by Remark 3.1. It can be proved that the Fourier transform of w ε with respect to ξ satisfies w ε ∈ C 0 (R d y ; L 1 (R d x )) and likewise for the Fourier transformation of w ε with respect to x ∈ R d . This allows us to define the integrals of w ε via a limiting process after convolving w ε with Gaussians; see Lions and Paul (1993) for more details.

Classical limit of Wigner transforms
The main point in the formulae given above is that the right-hand side of (3.4) involves only linear operations of w ε , which is compatible with weak limits. To this end, we recall the main result proved in Lions and Paul (1993) and Gérard et al. (1997).
Then, the set of Wigner functions is weak- * compact and thus, up to extraction of subsequences, is called the Wigner measure. If, in addition ∀ t : (ε∇u(t)) ∈ L 2 (R d ) uniformly with respect to ε, then we also have Note that although w ε (t) in general also takes negative values, its weak limit w(t) is indeed a non-negative measure in phase space.
are also often referred to as semiclassical measures and are closely related to the so-called H-measures used in homogenization theory (Tartar 1990). The fact that their weak limits are non-negative can be seen by considering the corresponding Husimi transformation, i.e., The Husimi transform w ε H is non-negative a.e. and has the same limit points as the Wigner function w ε ; see Markowich and Mauser (1993).
This result allows us to exchange limit and integration on the limit on the right-hand side of (3.4) to obtain The Wigner transformation and its associated Wigner measure are therefore highly useful tools for computing the classical limit of the expectation values of physical observables. In addition, it is proved in Lions and Paul (1993) and Gérard et al. (1997) that w(t, x, ξ) is the push-forward under the flow corresponding to the classical Hamiltonian H(x, ξ), i.e., where w in is the initial Wigner measure and F t : R 2d → R 2d is the phase space flow given byẋ In other words w(t, x, ξ) is a distributional solution of the classical Liouville equation in phase space, i.e., denotes the Poisson bracket. Note that in the case where H(x, ξ) is given by (3.5), this yields with characteristic equations given by the Newton trajectories: Strictly speaking we require V ∈ C 1 b (R d ), in order to define w as a distributional solution of (3.8). Note, however, in contrast to WKB techniques, the equation for the limiting Wigner measure (3.8) is globally well-posed, i.e., one does not experience problems of caustics. This is due to the fact that the Wigner measure w(t, x, ξ) lives in phase space. 134 S. Jin and P. Markowich and C. Sparber

Connection between Wigner measures and WKB analysis
A particularly interesting situation occurs for u ε in given in WKB form (2.18). The corresponding Wigner measure is found to be i.e., a mono-kinetic measure concentrated on the initial velocity v = ∇S in . In this case the phase space flow F t is projected onto physical space R d , yielding X t , the characteristic flow of the Hamilton-Jacobi equation (2.12). More precisely, the following result has been proved in Sparber et al. (2003).
Theorem 3.4. Let w(t, x, ξ) be the Wigner measure of the exact solution u ε to (2.1) with WKB initial data. Then if and only if ρ = |a| 2 and v = ∇S are smooth solutions of the leading-order WKB system given by (2.12) and (2.15).
This theorem links the theory of Wigner measures with the WKB approximation before caustics. After caustics, the Wigner measure is in general no longer mono-kinetic. However, it can be shown (Sparber et al. 2003, Jin and) that for generic initial data u ε in and locally away from caustics, the Wigner measure can be decomposed as (3.10) which is consistent with the multiphase WKB approximation given in (2.21).

Basic setting
A basic numerical scheme for solving linear partial differential equations is the well-known finite difference method (FD), to be discussed in this section; see, e.g., Strikwerda (1989) for a general introduction. In the following, we shall be mainly interested in its performance as ε → 0. To this end, we shall allow for more general Schrödinger-type PDEs in the form (Markowich and Poupaud 1999) where H W denotes the Weyl quantization of a classical real-valued phase space Hamiltonian H(x, ξ) ∈ C ∞ (R d × R d ), which is assumed to grow at most quadratically in x and ξ. For the following we assume that the symbol is a polynomial of order K ∈ N in ξ with C ∞ -coefficients H k (x), i.e., where k = (k 1 , . . . , k d ) ∈ N d denotes a multi-index with |k| := k 1 + · · · + k d . The differential operator H(x, εD x ) W can now be written as (4.2) In addition, assume that and, for simplicity, Under these conditions H(x, εD x ) W can be shown (Kitada 1980) to generate a unitary (strongly continuous) semi-group of operators U ε (t) = e −itH W /ε , which provides a unique global-in-time solution u ε = u ε (t) ∈ L 2 (R d ). Next, let Γ := γ = 1 r 1 + · · · + m r m : j ∈ Z for 1 ≤ j ≤ d ⊆ R d be the lattice generated by the linearly independent vectors r 1 , . . . , r d ∈ R d . For a multi-index k ∈ N d we construct a discretization of order N of the operator ∂ k x as follows: Here ∆x = h ∈ (0, h 0 ] is the mesh size, Γ k ⊆ Γ is the finite set of discretization points and a γ,k ∈ R are coefficients satisfying where δ ,k = 1 if = k, and zero otherwise. It is an easy exercise to show that the local discretization error of (4.3) is O(h N ) for all smooth functions if (B1) holds. For a detailed discussion of the linear problem (B1) (i.e., possible choices of the coefficients a γ,k ) we refer to Markowich and Poupaud (1999).

Spatial discretization
We shall now define the corresponding finite difference discretization of H(x, εD x ) W by applying (4.3) directly to (4.2). To this end, we let with = ε h being the ratio between the small semiclassical parameter ε and the mesh size h. Then we obtain the finite difference discretization of (4.2) in the form In view of (4.4), the discretization H h,ε (x, εD x ) W is seen to be a bounded operator on L 2 (R d ) and self-adjoint if We shall now collect several properties of such finite difference approximations, proved in Markowich, Pietra and Pohl (1999). We start with a spatial consistency result.
For a given ε > 0, choosing h such that = ε h → ∞ corresponds to asymptotically resolving the oscillations of wavelength O(ε) in the solution u ε (t, x) to the Schrödinger-type equation (4.1). In the case = const., i.e., putting a fixed number of grid points per oscillation, the symbol H h,ε (x, ξ) ≡ H (x, ξ) is independent of h and ε, i.e., (4.6) In the case ε,h→0 −→ 0, which corresponds to a scheme ignoring the εoscillations, we find and hence H h,ε (x, εD x ) W does not approximate H(x, εD x ) W . We thus cannot expect reasonable numerical results in this case (which will not be investigated further).

Temporal discretization and violation of gauge invariance
For the temporal discretization one can employ the Crank-Nicolson scheme with time step ∆t > 0. This is a widely used time discretization scheme for the Schrödinger equation, featuring some desirable properties (see below). We shall comment on the temporal discretizations below. The scheme reads subject to initial data u σ in = u ε in (x), where from now on, we shall denote the vector of small parameters by σ = (ε, h, ∆t). Note that the self-adjointness of H h,ε (x, εD x ) W implies that the operator Therefore the scheme (4.7) gives welldefined approximations u σ n for n = 1, 2, . . . if u ε in ∈ L 2 (R d ). Moreover, we remark that it is sufficient to evaluate (4.7) at x ∈ hΓ in order to obtain discrete equations for u σ n (hγ), γ ∈ Γ.
Remark 4.2. For practical computations, one needs to impose artificial 'far-field' boundary conditions. Their impact, however, will not be taken into account in the subsequent analysis.
By taking the L 2 scalar product of (4.7) with 1 2 u σ n+1 + 1 2 u σ n , one can directly infer the following stability result.
Lemma 4.3. The solution of (4.7) satisfies u σ n L 2 = u ε in L 2 , n = 0, 1, 2, . . . . In other words, the physically important property of mass-conservation also holds on the discrete level.
On the other hand, the scheme can be seen to violate the gauge invariance of (4.1). More precisely, one should note that expectation values of physical observables, as defined in (2.5), are invariant under the substitution (gauge transformation) In other words, the average value of the observable in the state u ε is equal to its average value in the state v ε . Remark 4.4. Note that in view of (3.1), the Wigner function is seen to also be invariant under this substitution, i.e., On the other hand, using this gauge transformation the Schrödinger equation (4.1) transforms to which implies that the zeroth-order term H 0 (x) in (4.2) is replaced by H 0 (x) + ω, while the other coefficients H k (x), k = 0, remain unchanged. In physical terms, H 0 (x) corresponds to a scalar (static) potential V (x). The corresponding force field obtained via F (x) = ∇H 0 (x) = ∇(H 0 (x) + ω) is unchanged by the gauge transformation and thus (4.8) can be considered (physically) equivalent to (4.1). The described situation, however, is completely different for the difference scheme outlined above. Indeed, a simple calculation shows that the discrete gauge transformation v σ n = u σ n e iωtn/ε does not commute with the discretization (4.7), up to adding a real constant to the potential. Thus, the discrete approximations of average values of observables depend on the gauging of the potential. In other words, the discretization method is not time-transverse-invariant.

Stability-consistency analysis for FD in the semiclassical limit
The consistency-stability concept of classical numerical analysis provides a framework for the convergence analysis of finite difference discretizations of linear partial differential equations. Thus, for ε > 0 fixed it is easy to prove that the scheme (4.7) is convergent of order N in space and order 2 in time if the exact solution u ε (t, x) is sufficiently smooth. Therefore, again for fixed ε > 0, we conclude convergence of the same order for average values of physical observables provided a(x, ξ) is smooth.
However, due to the oscillatory nature solutions to (4.1) the local discretization error of the finite difference schemes and, consequently, also the global discretization error, in general tend to infinity as ε → 0. Thus, the classical consistency-stability theory does not provide uniform results in the classical limit. Indeed, under the reasonable assumption that, for all multi-indices j 1 and j 2 ∈ N d , locally uniformly in t ∈ R, the classical stability-consistency analysis gives the following bound for the global L 2 -discretization error: The situation is further complicated by the fact that for any fixed t ∈ R, the solution u ε (t, ·) of (4.1) and its discrete counterpart u σ n in general converge only weakly in L 2 (R d ) as ε → 0, respectively, σ → 0. Thus, the limit processes ε → 0, σ → 0 do not commute with the quadratically nonlinear operation (2.5), needed to compute the expectation value of physical observables a[u ε (t)].
In practice, one is therefore interested in finding conditions on the mesh size h and the time step ∆t, depending on ε in such a way that the expectation values of physical observables in discrete form approximate a[u ε (t)] uniformly as ε → 0. To this end, let t n = n∆t, n ∈ N, and denote a σ (t n ) := a(·, εD x ) W u σ n , u σ n . The function a σ (t), t ∈ R, is consequently defined by piecewise linear interpolation of the values a σ (t n ). We seek conditions on h, k such that, for all and locally uniformly in t ∈ R. A rigorous study of this problem will be given by using the theory of Wigner measures applied in a discrete setting. Denoting the Wigner transformation (on the scale ε) of the finite difference solution u σ n by w σ (t n ) := w ε [u σ n ], and defining, as before, w σ (t), for any t ∈ R, by the piecewise linear interpolation of w σ (t n ), we conclude that (4.9) is equivalent to proving, locally uniformly in t, where w ε (t) is the Wigner transform of the solution u ε (t) of (4.1). We shall now compute the accumulation points of the sequence {w σ (t)} σ as σ → 0. We shall see that for any given subsequence {σ n } n∈N , the set of Wigner measures of the difference schemes depends decisively on the relative sizes of ε, h and ∆t. Clearly, in those cases in which µ = w, where w denotes the Wigner measure of the exact solution u ε (t), the desired property (4.10) follows. On the other hand (4.10) does not hold if the measures µ and w are different. Such a Wigner-measure-based study of finite difference schemes has been conducted in  and Markowich and Poupaud (1999). The main result given there is as follows.
Theorem 4.5. Fix a scale ε > 0 and denote by µ the Wigner measure of the discretization (4.7) as σ → 0. Then we have the following results.
The proof of this result proceeds similarly to the derivation of the phase space Liouville equation (3.7), in the continuous setting. Note that Theorem 4.5 implies that, as ε → 0, expectation values for physical observables in the state u ε (t), computed via the Crank-Nicolson finite difference scheme, are asymptotically correct only if both spatial and temporal oscillations of wavelength ε are accurately resolved.
Remark 4.6. Time-irreversible finite difference schemes, such as the explicit (or implicit) Euler scheme, behave even less well, as they require ∆t = o(ε 2 ) in order to guarantee asymptotically correct numerically computed observables; see .

Basic setting, first-and second-order splittings
As has been discussed before, finite difference methods do not perform well in computing the solution to semiclassical Schrödinger equations. An alternative is given by time-splitting trigonometric spectral methods, which will be discussed in this subsection; see also McLachlan and Quispel (2002) for a broad introduction to splitting methods. For the sake of notation, we shall introduce the method only in the case of one space dimension, d = 1.
Generalizations to d > 1 are straightforward for tensor product grids and the results remain valid without modifications. In the following, we shall therefore study the one-dimensional version of equation (2.1), i.e., We choose the spatial mesh size ∆x = h > 0 with h = (b − a)/M for some M ∈ 2N, and an ε-independent time step ∆t ≡ k > 0. The spatio-temporal grid points are then given by In the following, let u ε,n j be the numerical approximation of u ε (x j , t n ), for j = 1, . . . , M.
First-order time-splitting spectral method (SP1 ) The Schrödinger equation (5.1) is solved by a splitting method, based on the following two steps.
Step 1. From time t = t n to time t = t n+1 , first solve the free Schrödinger equation Step 2. On the same time interval, i.e., t ∈ [t n , t n+1 ], solve the ordinary differential equation (ODE) with the solution obtained from Step 1 as initial data for Step 2. Equation In Step 1, the linear equation (5.2) will be discretized in space by a (pseudo-) spectral method (see, e.g., Fornberg (1996) for a general introduction) and consequently integrated in time exactly. More precisely, one obtains at time where γ = 2πl b−a and u ε,n denote the Fourier coefficients of u ε,n , i.e., Note that the only time discretization error of this method is the splitting error, which is first-order in k = ∆t, for any fixed ε > 0.
Strang splitting (SP2 ) In order to obtain a scheme which is second-order in time (for fixed ε > 0), one can use the Strang splitting method , i.e., on the time interval [t n , t n+1 ] we compute with u ε, * denoting the Fourier coefficients of u ε, * given by u ε, * j = e iV (x j )k/2ε u ε,n j . Again, the overall time discretization error comes solely from the splitting, which is now (formally) second-order in ∆t = k for fixed ε > 0.
Remark 5.1. Extensions to higher-order (in time) splitting schemes can be found in the literature: see, e.g., Bao and Shen (2005). For rigorous investigations of the long time error estimates of such splitting schemes we refer to Faou (2007a, 2007b) and the references given therein.
In comparison to finite difference methods, the main advantage of such splitting schemes is that they are gauge-invariant: see the discussion in Section 4 above. Concerning the stability of the time-splitting spectral approximations with variable potential V = V (x), one can prove the following lemma (see Bao, Jin and Markowich (2002)), in which we denote U = (u 1 , . . . , u M ) and · l 2 the usual discrete l 2 -norm on the interval [a, b], i.e., Lemma 5.2. The time-splitting spectral schemes SP1 and SP2 are unconditionally stable, i.e., for any mesh size h and any time step k, we have where u ε,n int denotes the trigonometric polynomial interpolating In other words, time-splitting spectral methods satisfy mass-conservation on a fully discrete level.

Error estimate of SP1 in the semiclassical limit
To get a better understanding of the stability of spectral methods in the classical limit ε → 0, we shall establish the error estimates for SP1. Assume for some constant C m > 0. Under these assumptions it can be shown that the solution u ε = u ε (t, x) of (5.1) is (b−a)-periodic and smooth. In addition, we assume for all m, m 1 , m 2 ∈ N ∪ {0}. Thus, we assume that the solution oscillates in space and time with wavelength ε, but no smaller.
Remark 5.3. The latter is known to be satisfied if the initial data u ε in only invoke oscillations of wavelength ε (but no smaller).
Theorem 5.4. Let V (x) satisfy assumption (A) and let u ε (t, x) be a solution of (5.1) satisfying (B). Denote by u ε,n int the interpolation of the discrete approximation obtained via SP1. Then, if where C > 0 is independent of ε and m and G m > 0 is independent of ε, ∆x, ∆t. The proof of this theorem is given in Bao et al. (2002), where a similar result is also shown for SP2. Now let δ > 0 be a desired error bound such that holds uniformly in ε. Then Theorem 5.4 suggests the following meshing strategy on O(1) time and space intervals: where m ≥ 1 is an arbitrary integer, assuming that G m does not increase too fast as m → ∞. This meshing is already more efficient than what is needed for finite differences. In addition, as will be seen below, the conditions (5.5) can be strongly relaxed if, instead of resolving the solution u ε (t, x), one is only interested in the accurate numerical computation of quadratic observable densities (and thus asymptotically correct expectation values).
Example 5.5. This is an example from Bao et al. (2002). The Schrödinger equation (2.1) is solved with V (x) = 10 and the initial data The computational domain is restricted to [0, 1] equipped with periodic boundary conditions. Figure 5.1 shows the solution of the limiting position density ρ and current density j obtained by taking moments of w, satisfying the Liouville equation (3.8). This has to be compared with the oscillatory ρ ε and j ε , obtained by solving the Schrödinger equation (2.1) using SP2. As one can see, these oscillations are averaged out in the weak limits ρ, j.

Accurate computation of quadratic observable densities using time-splitting
We shall again invoke the theory of Wigner functions and Wigner measures.
To this end, let u ε (t, x) be the solution of (5.1) and let w ε (t, x, ξ) be the corresponding Wigner transform. Keeping in mind the results of Section 3, we see that the first-order splitting scheme SP1 corresponds to the following time-splitting scheme for the Wigner equation (3.2).
Step 1. For t ∈ [t n , t n+1 ], first solve the linear transport equation Step 2. On the same time interval, solve the non-local (in space) ordinary differential equation with initial data obtained from Step 1 above.
In (5.6), the only possible ε-dependence stems from the initial data. In addition, in (5.7) the limit ε → 0 can be easily carried out (assuming sufficient regularity of the potential V (x)) with k = ∆t fixed. In doing so, one consequently obtains a time-splitting scheme of the limiting Liouville equation (3.8) as follows.
Step 1. For t ∈ [t n , t n+1 ], solve Step 2. Using the outcome of Step 1 as initial data, solve, on the same time interval, Note that in this scheme no error is introduced other than the splitting error, since the time integrations are performed exactly.
These considerations, which can easily be made rigorous (for smooth potentials), show that a uniform time-stepping (i.e., an ε-independent k = ∆t) of the form combined with the spectral mesh size control given in (5.5) yields the following error: uniformly in ∆t as ε → 0. Essentially this implies that a fixed number of grid points in every spatial oscillation of wavelength ε combined with ε-independent time-stepping is sufficient to guarantee the accurate computation of (expectation values of) physical observables in the classical limit. This strategy is therefore clearly superior to finite difference schemes, which require k/ε → 0 and h/ε → 0, even if one only is interested in computing physical observables.
Remark 5.6. Time-splitting methods have been proved particularly successful in nonlinear situations: see the references given in Section 15.4 below.

Moment closure methods
We have seen before that a direct numerical calculation of u ε is numerically very expensive, particularly in higher dimensions, due to the mesh and time step constraint (5.5). In order to circumvent this problem, the asymptotic analysis presented in Sections 2 and 3 can be invoked in order to design asymptotic numerical methods which allow for an efficient numerical simulation in the limit ε → 0. The initial value problem (3.8)-(3.9) is the starting point of the numerical methods to be described below. Most recent computational methods are derived from, or closely related to, this equation. The main advantage is that (3.8)-(3.9) correctly describes the limit of quadratic densities of u ε (which in itself exhibits oscillations of wavelength O(ε)), and thus allows a numerical mesh size independent of ε. However, we face the following major difficulties in the numerical approximation.
(1) High-dimensionality. The Liouville equation (3.8) is defined in phase space, and thus the memory requirement exceeds the current computational capability in d ≥ 3 space dimensions.
(2) Measure-valued initial data. The initial data (3.9) form a delta measure and the solution at later time remains one (for a single-valued solution), or a summation of several delta functions (for a multivalued solution (3.10)).
In the past few years, several new numerical methods have been introduced to overcome these difficulties. In the following, we shall briefly describe the basic ideas of these methods.

The concept of multivalued solutions
In order to overcome the problem of high-dimensionality one aims to approximate w(t, x, p) by using averaged quantities depending only on t, x. This is a well-known technique in classical kinetic theory, usually referred to as moment closure. A basic example is provided by the result of Theorem 3.4, which tells us that, as long as the WKB analysis of Section 2.2 is valid (i.e., before the appearance of the first caustic), the Wigner measure is given by a mono-kinetic distribution in phase space, i.e., where one identifies ρ = |a| 2 and v = ∇S. The latter solve the pressure-less Euler system which, for smooth solutions, is equivalent to the system of transport equation (2.14) coupled with the Hamilton-Jacobi equation (2.12), obtained through the WKB approximation. Thus, instead of solving the Liouville equation in phase space, one can as well solve the system (6.1), which is posed on physical space R t × R d x . Of course, this can only be done until the appearance of the first caustic, or, equivalently, the emergence of shocks in (6.1).
In order to go beyond that, one might be tempted to use numerical methods based on the unique viscosity solution (see Crandall and Lions (1983)) for (6.1). However, the latter does not provide the correct asymptotic description -the multivalued solution -of the wave function u ε (t, x) beyond caustics. Instead, one has to pass to so-called multivalued solutions, based on higher-order moment closure methods. This fact is illustrated in Figure 6.1, which shows the difference between viscosity solutions and multivalued solutions. Figures 6.1(a) and 6.1(b) are the two different solutions for the following eikonal equation (in fact, the zero-level set of S): This equation, corresponding to H(ξ) = |ξ|, arises in the geometric optics limit of the wave equation and models two circular fronts moving outward in the normal direction with speed 1; see Osher and Sethian (1988). As one can see, the main difference occurs when the two fronts merge. Similarly, Figures 6.1(c) and 6.1(d) show the difference between the viscosity and the multivalued solutions to the Burgers equation This is simply the second equation in the system (6.1) for V (x) = 0, written in divergence form. The solution begins as a sinusoidal function and then forms a shock. Clearly, the solutions are different after the shock formation.

Moment closure
The moment closure idea was first introduced by Brenier and Corrias (1998) in order to define multivalued solutions to the Burgers equation, and seems to be the natural choice in view of the multiphase WKB expansion given in (2.21). The method was then used numerically in Engquist and Runborg (1996) (see also Engquist and Runborg (2003) for a broad review) and Gosse (2002) to study multivalued solutions in the geometrical optics regime of hyperbolic wave equations. A closely related method is given in Benamou (1999), where a direct computation of multivalued solutions to Hamilton-Jacobi equations is performed. For the semiclassical limit of the Schrödinger equation, this was done in  and then Gosse et al. (2003).
In order to describe the basic idea, let d = 1 and define i.e., the th moment (in velocity) of the Wigner measure. By multiplying the Liouville equation (3.8) by ξ and integrating over R ξ , one obtains the following moment system: Note that this system is not closed , since the equation determining the th moment involves the ( + 1)th moment.
The δ-closure As already mentioned in (3.10), locally away from caustics the Wigner measure of u ε as ε → 0 can be written as where the number of velocity branches J can in principle be determined a priori from ∇S in (x). For example, in d = 1, it is the total number of inflection points of v(0, x): see Gosse et al. (2003). Using this particular form (6.5) of w with L = 2J provides a closure condition for the moment system above. More precisely, one can express the last moment m 2J as a function of all of the lower-order moments , i.e., This consequently yields a system of 2J × 2J equations (posed in physical space), which effectively provides a solution of the Liouville equation, before the generation of a new phase, yielding a new velocity v j , j > J. It was shown in  that this system is only weakly hyperbolic, in the sense that the Jacobian matrix of the flux is a Jordan block, with only J distinct eigenvalues v 1 , v 2 , . . . , v J . This system is equivalent to J pressureless gas equations (6.1) for (ρ j , v j ) respectively. In  the explicit flux function g in (6.6) was given for J ≤ 5. For larger J a numerical procedure was proposed for evaluating g. Since the moment system is only weakly hyperbolic, with phase jumps which are under-compressive shocks (Gosse et al. 2003), standard shock-capturing schemes such as the Lax-Friedrichs scheme and the Godunov scheme face severe numerical difficulties as in the computation of the pressure-less gas dynamics: see Bouchut, Jin and Li (2003), Engquist and Runborg (1996) or Jiang and Tadmor (1998). Following the ideas of Bouchut et al. (2003) for the pressure-less gas system, a kinetic scheme derived from the Liouville equation (3.8) with the closure condition (6.6) was used in  for this moment system.
The Heaviside closure Another type of closure was introduced by Brenier and Corrias (1998) using the following ansatz, called H-closure, to obtain the J-branch velocities v j , with j = 1, . . . , J: This type of closure condition for (3.8) arises from an entropy-maximization principle: see Levermore (1996). Using (6.7), one arrives at (6.6) with L = J. The explicit form of the corresponding function g(m 0 , . . . , m 2J−1 ) for J < 5 is available analytically in Runborg (2000). Note that this method decouples the computation of velocities v j from the densities ρ j . In fact, to obtain the latter, Gosse (2002) has proposed solving the following linear conservation law (see also Gosse and James (2002) and Gosse et al. (2003)): The numerical approximation to this linear transport with variable or even discontinuous flux is not straightforward. Gosse et al. (2003) used a semi-Lagrangian method that uses the method of characteristics, requiring the time step to be sufficiently small for the case of non-zero potentials. The corresponding method is usually referred to as H-closure. Note that in d = 1 the H-closure system is a non-strictly rich hyperbolic system, whereas the δ-closure system described before is only weakly hyperbolic. Thus one expects a better numerical resolution from the H-closure approach, which, however, is much harder to implement in the higher dimension. In d = 1, the mathematical equivalence of the two moment systems was proved in Gosse et al. (2003).
Remark 6.1. Multivalued solutions also arise in the high-frequency approximation of nonlinear waves, for example, in the modelling of electron transport in vacuum electronic devices: see, e.g., Granastein, Parker and Armstrong (1999). There the underlying equations are the Euler-Poisson equations, a nonlinearly coupled hyperbolic-elliptic system. The multivalued solution of the Euler-Poisson system also arises for electron sheet initial data, and can be characterized by a weak solution of the Vlasov-Poisson equation: see Majda, Majda and Zheng (1994). Similarly, the work of Li, Wöhlbier, Jin and Booske (2004) uses the moment closure ansatz (6.6) for the Vlasov-Poisson system; see also Wöhlbier, Jin and Sengele (2005). For multivalued (or multiphase) solution of the semiclassical limit of nonlinear dispersive waves using the closely related method of Whitham's modulation theory, we refer to Whitham (1974) and Flaschka, Forest and McLaughlin (1980). Finally, we mention that multivalued solutions also arise in supply chain modelling: see, e.g., Armbruster, Marthaler and Ringhofer (2003).
In summary, the moment closure approach yields an Eulerian method defined in the physical space which offers a greater efficiency compared to the computation in phase space. However, when the number of phases J ∈ N becomes very large and/or in dimensions d > 1, the moment systems become very complex and thus difficult to solve. In addition, in high space dimensions, it is very difficult to estimate a priori the total number of phases needed to construct the moment system. Thus it remains an interesting and challenging open problem to develop more efficient and general physicalspace-based numerical methods for the multivalued solutions.

Eulerian approach
Level set methods have been recently introduced for computing multivalued solutions in the context of geometric optics and semiclassical analysis. These methods are rather general, and applicable to any (scalar) multi-dimensional quasilinear hyperbolic system or Hamilton-Jacobi equation (see below). We shall now review the basic ideas, following the lines of Jin and Osher (2003). The original mathematical formulation is classical; see for example Courant and Hilbert (1962).

Computation of the multivalued phase
Consider a general d-dimensional Hamilton-Jacobi equation of the form For example, in present context of semiclassical analysis for Schrödinger equations, while for applications in geometrical optics (i.e., the high-frequency limit of the wave equation), with c(x) denoting the local sound (or wave) speed. Introducing, as before, a velocity v = ∇S and taking the gradient of (7.1), one gets an equivalent equation (at least for smooth solutions) in the form (Jin and Xin 1998): In other words, the (intersection of the) zero-level sets of all {φ j } d j=1 yield the graph of the multivalued solution v j (t, x) of (7.2). Using (7.2) it is easy to see that φ j solves the following initial value problem: which is simply the phase space Liouville equation. Note that in contrast to (7.2), this equation is linear and thus can be solved globally in time. In doing so, one obtains, for all t ∈ R, the multivalued solution to (7.2) needed in the asymptotic description of physical observables. See also Cheng, Liu and Osher (2003).

Computation of the particle density
It remains to compute the classical limit of the particle density ρ(t, x). To do so, a simple idea was introduced in Jin, Liu, Osher and Tsai (2005a). This method is equivalent to a decomposition of the measure-valued initial data (3.9) for the Liouville equation. More precisely, a simple argument based on the method of characteristics (see Jin, Liu, Osher and Tsai (2005b)) shows that the solution to (3.8)-(3.9) can be written as where φ j (t, x, ξ) ∈ R n , j = 1, . . . , d, solves (7.3) and the auxiliary function ψ(t, x, ξ) again satisfies the Liouville equation (3.8), subject to initial data: The first two moments of w with respect to ξ (corresponding to the particle ρ and current density J = ρu) can then be recovered through Thus the only time one has to deal with the delta measure is at the numerical output, while during the time evolution one simply solves for φ j and ψ, both of which are smooth L ∞ -functions. This avoids the singularity problem mentioned earlier, and gives numerical methods with much better resolution than solving (3.8)-(3.9) directly, e.g., by approximating the initial delta function numerically. An additional advantage of this level set approach is that one only needs to care about the zero-level sets of φ j . Thus the technique of local level set methods developed in Adalsteinsson and Sethian (1995) and Peng, Merriman, Osher, Zhao and Kang (1999) can be used.
One thereby restricts the computational domain to a narrow band around the zero-level set, in order to reduce the computational cost to O (N ln N ), for N computational points in the physical space. This is a nice alternative for dimension reduction of the Liouville equation. When solutions for many initial data need to be computed, fast algorithms can be used: see Fomel and Sethian (2002) or Ying and Candès (2006).
, and WKB initial data corresponding to Figure 7.1 shows the time evolution of the velocity and the corresponding density computed by the level set method described above. The velocity eventually develops some small oscillations with higher frequency, which require a finer grid to resolve.
Remark 7.2. The outlined ideas have been extended to general linear symmetric hyperbolic systems in Jin et al. (2005a). So far, however, level set methods have not been formulated for nonlinear equations, except for the one-dimensional Euler-Poisson equations (Liu and Wang 2007), where a three-dimensional Liouville equation has to be used in order to calculate the corresponding one-dimensional multivalued solutions.

The Lagrangian phase flow method
While the Eulerian level set method is based on solving the Liouville equation (3.8) on a fixed mesh, the Lagrangian (or particle) method , is based on solving the Hamiltonian system (3.6), which is simply the characteristic flow of the Liouville equation (3.7). In geometric optics this idea is referred to as ray tracing (Cervený 2001), and the curves x(t, y, ζ), ξ(t, y, ζ) ∈ R d , obtained by solving (3.6), are usually called bi-characteristics. Remark 7.3. Note that finding an efficient way to numerically solve Hamiltonian ODEs, such as (3.6), is a problem of great (numerical) interest in its own right: see, e.g., Leimkuhler and Reich (2004). Here we shall briefly describe a fast algorithm, called the phase flow method in Ying and Candès (2006), which is very efficient if multiple initial data, as is often the case in practical applications, are to be propagated by the Hamiltonian flow (3.6). Let F t : R 2d → R 2d be the phase flow defined by F t (y, ζ) = (x(t, y, ζ), ξ(t, y, ζ)), t ∈ R.
For the autonomous ODEs, such as (3.6), a key property of the phase map is the one-parameter group structure, F t • F s = F t+s .
Instead of integrating (3.6) for each individual initial condition (y, ζ), up to, say, time T , the phase flow method constructs the complete phase map F T . To this end, one first constructs the F t for small times using standard ODE integrators and then builds up the phase map for larger times via a local interpolation scheme together with the group property of the phase flow. Specifically, fix a small time τ > 0 and suppose that T = 2 n τ .
Step 1. Begin with a uniform or quasi-uniform grid on M.
Step 2. Compute an approximation of the phase map F τ at time τ . The value of F τ at each grid point is computed by applying a standard ODE or Hamiltonian integrator with a single time step of length τ . The value of F τ at any other point is defined via a local interpolation.
When the algorithm terminates, one obtains an approximation of the whole phase map at time T = 2 n τ . This method is clearly much faster than solving each for initial condition independently.

Gaussian beam methods: Lagrangian approach
A common numerical problem with all numerical approaches based on the Liouville equation with mono-kinetic initial data (3.8)-(3.9) is that the particle density ρ(t, x) blows up at caustics. Another problem is the loss of phase information when passing through a caustic point, i.e., the loss of the Keller-Maslov index (Maslov 1981). To this end, we recall that the Wigner measure only sees the gradient of the phase: see (3.10). The latter can be fixed by incorporating this index into a level set method as was done in Jin and Yang (2008)). Nevertheless, one still faces the problem that any numerical method based on the Liouville equation is unable to handle wave interference effects. The Gaussian beam method (or Gaussian wave packet approach, as it is called in quantum chemistry: see Heller (2006)), is an efficient approximate method that allows an accurate computation of the wave amplitude around caustics, and in addition captures the desired phase information. This now classical method was developed in Popov (1982), Ralston (1982) and Hill (1990), and has seen increasing activity in recent years. In the following, we shall describe the basic ideas, starting with its classical Lagrangian formulation.

Lagrangian dynamics of Gaussian beams
Similar to the WKB method, the approximate Gaussian beam solution is given in the form where the variable y = y(t, y 0 ) will be determined below and the phase T (t, x, y) is given by This is reminiscent of the Taylor expansion of the phase S around the point y, upon identifying p = ∇S ∈ R d , M = ∇ 2 S, the Hessian matrix. The idea is now to allow the phase T to be complex-valued (in contrast to WKB analysis), and choose the imaginary part of M ∈ C n×n positive definite so that (8.1) indeed has a Gaussian profile. Plugging the ansatz (8.1) into the Schrödinger equation (2.1), and ignoring the higher-order terms in both ε and (y − x), one obtains the following system of ODEs: where p, V, M, S and A have to be understood as functions of (t, y(t, y 0 )). The latter defines the centre of a Gaussian beam. Equations (8.2)-(8.4) can be considered as the Lagrangian formulation of the Gaussian beam method, with (8.2) furnishing a classical ray-tracing algorithm. We further note that (8.3) is a Riccati equation for M . We state the main properties of (8.3), (8.4) in the following theorem, the proof of which can be found in Ralston (1982) (see also Jin, Wu and Yang (2008b)).
Theorem 8.1. Let P (t, y(t, y 0 )) and R(t, y(t, y 0 )) be the (global) solutions of the equations dP dt = R, dR dt = −(∇ 2 y V )P, (8.5) with initial conditions where Id is the identity matrix and Im(M (0, y 0 )) is positive definite. Assume that M (0, y 0 ) is symmetric. Then, for each initial position y 0 , we have the following properties.
(4) The Hamiltonian H = 1 2 |p| 2 + V is conserved along the y-trajectory, as is (A 2 det P ), i.e., A(t, y(t, y 0 )) can be computed via where the square root is taken as the principal value.
In particular, since (A 2 det P ) is a conserved quantity, we infer that A does not blow up along the time evolution (provided it is initially bounded).

Lagrangian Gaussian beam summation
It should be noted that a single Gaussian beam given by (8.1) is not an asymptotic solution of (2.1), since its L 2 (R 2d )-norm goes to zero, in the classical limit ε → 0. Rather, one needs to sum over several Gaussian beams, the number of which is O(ε −1/2 ). This is referred to as the Gaussian beam summation; see for example Hill (1990). In other words, one first needs to approximate given initial data through Gaussian beam profiles. For WKB initial data (2.18), a possible way to do so is given by the next theorem, proved by Tanushev (2008).
Theorem 8.2. Let the initial data be given by , r θ ≥ 0 is a truncation function with r θ ≡ 1 in a ball of radius θ > 0 around the origin, and C is a constant related to θ.
In view of Theorem 8.2, one can specify the initial data for (8.2)-(8.4) as Then, the Gaussian beam solution approximating the exact solution of (2.1) is given by x, y(t, y 0 )) dy 0 .
In discretized form this reads where the y j 0 are equidistant mesh points, and N y 0 is the number of the beams initially centred at y j 0 .
Remark 8.3. Note that the cut-off error introduced via r θ becomes large when the truncation parameter θ is taken too small. On the other hand, a big θ for wide beams makes the error in the Taylor expansion of T large. As far as we know, it is still an open mathematical problem to determine an optimal size of θ when beams spread. However, for narrow beams one can take a fairly large θ, which makes the cut-off error almost zero. For example, a one-dimensional constant solution can be approximated by in which r θ ≡ 1.

Higher-order Gaussian beams
The above Gaussian beam method can be extended to higher order in ε: see Tanushev (2008),  and Liu and Ralston (2010). For notational convenience we shall only consider the case d = 1. Consider the Schrödinger equation (2.1) with initial data Let a ray y(t, y 0 ) start at a point y 0 ∈ R. Expand S in (x) in a Taylor series around y 0 : Then, a single kth-order Gaussian beam takes the form where the phase is given by and the amplitude reads Here the bi-characteristic curves (y(t, y 0 ), p(t, y(t, y 0 ))) satisfy the Hamiltonian system (8.2) with initial data y(0, y 0 ) = y 0 , p(0, y 0 ) = ∂ y S 0 (y 0 ).
In d = 1, the kth-order Gaussian beam superposition is thus formed by where, as before, r θ ∈ C ∞ 0 (R; R) is some cut-off function. For this type of approximation, the following theorem was proved in Liu, Runborg and Tanushev (2011).
Theorem 8.4. If u ε (t, x) denotes the exact solution to the Schrödinger equation (2.1) and u ε G,k is the kth-order Gaussian beam superposition, then for any T > 0.

Eulerian dynamics of Gaussian beams
The Gaussian beam method can be reformulated in an Eulerian framework. To this end, let us first define the linear Liouville operator as In addition, we shall denote where φ j is the level set function defined in (7.1) and (7.3). Using this, Jin et al. (2005b) and Jin and Osher (2003) showed that one can obtain from the original Lagrangian formulation (8.2)-(8.4) the following (inhomogeneous) Liouville equations for velocity, phase and amplitude, respectively: In addition, if one introduces the quantity (Jin et al. 2005b) then f (t, y, ξ) again satisfies the Liouville equation, i.e., Two more inhomogeneous Liouville equations, which are the Eulerian version of (8.5) for P and R, were introduced in Leung, Qian and Burridge (2007) to construct the Hessian matrix. More precisely, one finds Note that the equations (9.1)-(9.4) are real , while (9.5) and (9.6) are complex and consist of 2n 2 equations.

Gaussian beam dynamics using complex level set functions
In  the following observation was made. Taking the gradient of the equation (9.1) with respect to y and ξ separately, we have Comparing (9.5)-(9.6) with (9.7)-(9.8), one observes that −∇ y Φ and ∇ ξ Φ satisfy the same equations as R and P , respectively. Since the Liouville operator is linear, one can allow Φ to be complex-valued and impose for −∇ y Φ, ∇ ξ Φ the same initial conditions as for R and P , respectively. By doing so, R = −∇ x Φ, P = ∇ ξ Φ hold true for any time t ∈ R. In view of (8.6) and (8.10), this suggests the following initial condition for Φ: With this observation, one can now solve (9.1) for complex Φ, subject to initial data (9.9). Then the matrix M can be constructed by where the velocity v = ∇ y S is given by the intersection of the zero-level contours of the real part of Φ, i.e., for each component φ j , Re(φ j (t, y, ξ)) = 0, at ξ = v(t, y) = ∇ y S. (9.11) Note that in order to compute v, S and M , one only needs to solve d complex-valued homogeneous Liouville equations (9.1). The Eulerian level set method of  (in complex phase space) can then be summarized as follows.
Step 1. Solve (9.1) for Φ complex , with initial condition (9.9), and obtain the velocity v by the intersection of the zero-level sets of Re φ j , j = 1, . . . , n.
Step 2. Use −∇ y Φ and ∇ ξ Φ to construct M by (9.10) (note that these quantities are already available from the first step after discretizing the Liouville equation for Φ).
Step 3. Integrate the velocity v along the zero-level sets (Gosse 2002, Jin andYang 2008) to get the phase S. To do so, one performs a numerical integration following each branch of the velocity. The integration constants are obtained by the boundary condition and the fact that the multivalued phase is continuous when passing from one branch to the other. For example, if one considers a bounded domain [a, b] in space dimension d = 1, the phase function is given by For more details on this and its extension to higher dimensions, see Jin and Yang (2008).
Step 4. Solve (9.4) with the initial condition f 0 (y, ξ) = A 2 0 (y, ξ). Then amplitude A is given by where the square root has to be understood as the principal value. (We also refer to Jin, Wu and Yang (2011) for a more elaborate computation of A.) Note that all functions appearing in Steps 2-4 only need to be solved locally around the zero-level sets of Re φ j , j = 1, . . . , n. Thus, the entire algorithm can be implemented using the local level set methods of Osher, Cheng, Kang, Shim and Tsai (2002) and Peng et al. (1999). For a given mesh size ∆y, the computational cost is therefore O((∆y) −d ln(∆y) −1 ), about the same as the local level set methods for geometrical optics computation: see Jin et al. (2005b).
Remark 9.1. If one is only interested in computing the classical limit of (the expectation values of) physical observables, one observes that the only term in (9.12) which affects a quadratic observable density for fixed time t is x a v(t, y) dy. Thus, as long as one is only interested in physical observables, one can simply take S(t, x) = x a v(t, y) dy (9.14) in the numerical simulations.
That M and A are indeed well-defined via (9.10) and (9.13) is justified by the following theorem (which can be seen as the Eulerian version of Theorem 8.1), proved in .

Eulerian Gaussian beam summation
As before, we face the problem of Gaussian beam summation, i.e., in order to reconstruct the full solution u ε a superposition of single Gaussian beams has to be considered. To this end, we define a single Gaussian beam, obtained through the Eulerian approach, bỹ where A is solved via (9.13), and Then, the wave function is constructed via the following Eulerian Gaussian beam summation formula (Leung et al. 2007): which is consistent with the Lagrangian summation formula (8.2). Indeed, the above given double integral for u ε G can be evaluated as a single integral in y as follows. We again denote by v j , j = 1, . . . , J the jth branch of the multivalued velocity and writẽ However, since det Re(∇ ξ Φ) = 0 at caustics, a direct numerical integration of (9.16) loses accuracy around caustics. To get better accuracy, one can split (9.16) into two parts, where Ω 1 := y : det(Re(∇ p φ)(t, y, p j )) ≥ τ , with τ being a small parameter. The latter is chosen sufficiently small to minimize the cost of computing (9.18), yet large enough to make I 1 a regular integral.
The regular integral I 1 can then be approximated by a standard quadrature rule, such as the trapezoid quadrature rule, while the singular integral I 2 is evaluated by the semi-Lagrangian method introduced in Leung et al. (2007). Remark 9.3. When the velocity contours are complicated due to large numbers of caustics, implementation of the local semi-Lagrangian method is hard. In such situations one can use a discretized δ-function method for numerically computing (9.18), as was done in Wen (2010). In this method one needs to numerically solve (9.2) in order to obtain the phase function, since all values of φ j near the support of δ(Re(φ j )) are needed to evaluate (9.18).
Example 9.4. This is an example from . It considers the free motion of particles in d = 1 with V (x) = 0. The initial conditions for the Schrödinger equation (2.1) are induced by (5x). Figure 9.1 shows the l ∞ -errors between the square modulus of u ε , the exact solution of the Schrödinger equation (2.1), and the approximate solution constructed by: (i) the level set method described in Section 7, (ii) the level set method with Keller-Maslov index built in (Jin and Yang 2008), and (iii) the Eulerian Gaussian beam method described above. As one can see, method (ii) improves the geometric optics solution (i) away from caustics, while the Gaussian beam method offers a uniformly small error even near the caustics.
Compared to the Lagrangian formulation based on solving the ODE system (8.2)-(8.4), the Eulerian Gaussian beam method has the advantage of maintaining good numerical accuracy since it is based on solving PDEs on fixed grids. Moreover, higher-order (in ε) Eulerian Gaussian beam methods have been constructed: see Liu and Ralston (2010) and .

Frozen Gaussian approximations
The construction of the Gaussian beam approximation is based on the truncation of the Taylor expansion of the phase T (t, x, y) (8.1) around the beam centre y. Hence it loses accuracy when the width of the beam becomes large, i.e., when the imaginary part of M (t, y) in (8.1) becomes small so that the Gaussian function is no longer localized . This happens, for example, when the solution of the Schrödinger equation spreads (namely, the ray determined by (8.2) diverges), which can be seen as the time-reversed situation of caustic formation. The corresponding loss in the numerical computation can be overcome by re-initialization every once in a while: see Tanushev et al. (2009), Ariel et al. (2011, Qian and Ying (2010) and Yin and Zheng (2011). However, this approach increases the computational complexity in particular, when beams spread quickly.
The frozen Gaussian approximation (as it is referred to in quantum chemistry), first proposed in Heller (1981), uses Gaussian functions with fixed widths to approximate the exact solution u ε . More precisely, instead of using Gaussian beams only in the physical space, the frozen Gaussian approximation uses a superposition of Gaussian functions in phase space. That is why the method is also known by the name coherent state approximation. To this end, we first decompose the initial data into several Gaussian functions in phase space, and then propagate the centre of each function (y(t), p(t)) along the Hamiltonian flow (8.2), subject to initial data at (y 0 , p 0 ). The frozen Gaussian beam solution takes the form a(t, y 0 , p 0 )ψ(y 0 , p 0 ) e ( ip(t)·(x−y(t))− 1 2 |x−y(t)| 2 )/ε dp 0 dy 0 , where the complex-valued amplitude a(t, y 0 , p 0 ), is the so-called Herman-Kluk pre-factor : see Herman and Kluk (1984). Since the width of the Gaussians is fixed, one does not encounter the problem of beam spreading here. However, since this method is based in phase space, the computational cost is considerably higher than the standard Gaussian beam methods. For subsequent developments in this direction, see Herman and Kluk (1984), Kay (1994Kay ( , 2006, Robert (2010), Swart and Rousse (2009) and Lu and Yang (2011).

Asymptotic methods for discontinuous potentials
Whenever a medium is heterogeneous, the potential V can be discontinuous, creating a sharp potential barrier or interface where waves can be partially reflected and transmitted (as in the Snell-Descartes Law of Refraction). This gives rise to new mathematical and numerical challenges not present in the smooth potential case. Clearly, the semiclassical limit (3.8)-(3.9) does not hold at the barrier. Whenever V is discontinuous, the Liouville equation (3.8) contains characteristics which are discontinuous or even measure-valued. To this end, we recall that the characteristic curves (x(t), ξ(t)) are determined by the Hamiltonian system (3.6). The latter is a nonlinear system of ODEs whose right-hand side is not Lipschitz due to the singularity in ∇ x H(x, ξ), and thus the classical well-posedness theory for the Cauchy problem of the ODEs fails. Even worse, the coefficients in the Liouville equation are in general not even BV (i.e., of bounded variation), for which almost everywhere solutions were introduced by DiPerna and Lions (1989) and Ambrosio (2004). Analytical studies of semiclassical limits in situations with interface were carried out by, e.g., Bal, Keller, Papanicolaou and Ryzhik (1999b), Miller (2000), Nier (1995Nier ( , 1996 and Benedetto, Esposito and Pulvirenti (2004) using more elaborate Wigner transformation techniques, such as two-scale Wigner measures.

The interface condition
In order to allow for discontinuous potentials, one first needs to introduce a notion of solutions of the underlying singular Hamiltonian system (3.6). This can be done by providing interface conditions for reflection and transmission, based on Snell's law. The solution so constructed will give the correct transmission and reflection of waves through the barrier, obeying the laws of geometrical optics.  . However, in this section we will not discuss the diffusive interface, which was treated analytically in Bal et al. (1999b) and numerically in . Jin and Wen (2006a) provide an interface condition connecting the Liouville equations at both sides of a given (sharp) interface. Let us focus here only on the case of space dimension d = 1 and consider a particle moving with velocity ξ > 0 towards the barrier. The interface condition at a given fixed time t is given by

An Eulerian interface condition
Here the superscripts '±' represent the right and left limits of the quantities, α T ∈ [0, 1] and α R ∈ [0, 1] are the transmission and reflection coefficients respectively, satisfying α R + α T = 1. For a sharp interface x + = x − , however, ξ + and ξ − are connected by the Hamiltonian preservation condition The latter is motivated as follows. In classical mechanics, the Hamiltonian H = 1 2 ξ 2 + V (x) is conserved along the particle trajectory, even across the barrier. In this case, α T , α R = 0 or 1, i.e., a particle can be either completely transmitted or completely reflected. In geometric optics (corresponding to H(x, ξ) = c(x)|ξ|), condition (10.2) is equivalent to Snell's Law of Refraction for a flat interface (Jin and Wen 2006b), i.e., waves can be partially transmitted or reflected.
Remark 10.1. In practical terms, the coefficients α T and α R are determined from the original Schrödinger equation (2.1) before the semiclassical limit is taken. Usually one invokes stationary scattering theory to do so. Thus (10.1) represents a multiscale coupling between the (macroscopic) Liouville equation and the (microscopic) Schrödinger (or wave) equation. Furthermore, by incorporating the diffraction coefficients, determined from the geometrical theory of diffraction developed in Keller and Lewis (1995), into the interface condition, one could even simulate diffraction phenomena near boundaries, interfaces or singular geometries (Jin and Yin 2008a. The well-posedness of the initial value problem of the singular Liouville equation with the interface condition (10.1) was established in Jin and Wen (2006a), using the method of characteristics. To determine a solution at (x, p, t) one traces back along the characteristics determined by the Hamiltonian system (3.6) until hitting the interface. At the interface, the solution bifurcates according to the interface condition (10.1). One branch of the solution thereby corresponds to the transmission of waves and the other to the reflection. This process continues until one arrives at t = 0. The interface condition (10.1) thus provides a generalization of the method of characteristics.

A Lagrangian Monte Carlo particle method for the interface
A notion of the solution of the (discontinuous) Hamiltonian system (3.6) was introduced in Jin (2009) (see also Jin and Novak (2006)) using a probability interpretation. One thereby solves the system (3.6) using a standard ODE or Hamiltonian solver, but at the interface, the following Monte Carlo solution can be constructed (we shall only give a solution in the case of ξ − > 0; the other case is similar).
(1) With probability α R , the particle (wave) is reflected with (2) With probability α T , the particle (wave) is transmitted, with x → x, ξ + obtained from ξ − using (10.2). (10.4) Although the original problem is deterministic, this probabilistic solution allows one to go beyond the interface with the new value of (x, ξ) defined in (10.3)-(10.4). This is the Lagrangian formulation of the solution determined by using the interface condition (10.1). This is the basis of a (Monte-Carlobased) particle method for thin quantum barriers introduced in Jin and Novak (2007).

Modification of the numerical flux at the interface
A typical one-dimensional semi-discrete finite difference method for the Liouville equation (3.7) is Here w ij is the cell average or pointwise value of w(t, x i , ξ j ) at fixed t. The numerical fluxes w i+ 1 2 ,j , w i,j+ 1 2 are typically defined by a (first-, or higherorder) upwind scheme, and DV i is some numerical approximation of ∂ When V is discontinuous, such schemes face difficulties when the Hamiltonian is discontinuous, since ignoring the discontinuity of V in the actual numerical computation will result in solutions which are inconsistent with the notion of the physically relevant solution, defined in the preceding subsection. Even with a smoothed Hamiltonian, it is usually impossible (at least in the case of partial transmission and reflection) to obtain transmission and reflection with the correct transmission and reflection coefficients. A smoothed V will also give a severe time step constraint like ∆t ∼ O(∆x∆ξ): see, e.g., Cheng, Kang, Osher, Shim and Tsai (2004). This is a parabolic-type CFL condition, despite the fact that we are solving a hyperbolic problem.
A simple method to solve this problem was introduced by Jin andWen (2005, 2006a). The basic idea is to build the interface condition (10.1) into the numerical flux , as follows. Assume V is discontinuous at x i+1/2 . First one should avoid discretizing V across the interface at x i+1/2 . One possible discretization is where, for example, The numerical flux in the ξ direction, W i,j±1/2 , can be the usual numerical flux (for example, the upwind scheme or its higher-order extension). To define the numerical flux w ± i+1/2,j , without loss of generality, consider the case ξ j > 0. Using the upwind scheme, w − while ξ − is obtained from (10.2) with ξ + j = ξ j . Since ξ − may not be a grid point, one has to define it approximately. A simple approach is to locate the two cell centres that bound it, then use a linear interpolation to approximate the needed numerical flux at ξ − j . The case of ξ j < 0 is treated similarly. The detailed algorithm to generate the numerical flux is given in Jin andWen (2005, 2006a).
This numerical scheme overcomes the aforementioned analytic and numerical difficulties. In particular, it possesses the following properties.
(1) It produces the correct physical solution crossing the interface (as defined in the previous subsection). In particular, in the case of geometric optics, this solution is consistent with the Snell-Descartes Law of Refraction at the interface. (2) It allows a hyperbolic CFL condition ∆t = O(∆x, ∆ξ).
The idea outlined above has its origin in so-called well-balanced kinetic schemes for shallow water equations with bottom topography: see Perthame and Simeoni (2001). It has been applied to computation of the semiclassical limit of the linear Schrödinger equation with potential barriers in Jin and Wen (2005), and the geometrical optics limit with complete transmission/reflection in Jin and Wen (2006b), for thick interfaces, and Jin and Wen (2006a), for sharp interfaces. Positivity and both l 1 and l ∞ stabilities were also established, under a hyperbolic CFL condition. For piecewise constant Hamiltonians, an l 1 -error estimate of the first-order finite difference of this type was established in Wen (2009), following Wen and Jin (2009).
Remark 10.2. Let us remark that this approach has also been extended to high-frequency elastic waves (Jin and Liao 2006), and high-frequency waves in random media with diffusive interfaces . When the initial data are measure-valued, such as (3.9), the level set method introduced in Jin et al. (2005b) becomes difficult for interfaces where waves undergo partial transmissions and reflections, since one needs to increase the number of level set functions each time a wave hits the interface. A novel method to get around this problem has been introduced in Wei, Jin, Tsai and Yang (2010). It involves two main ingredients.
(1) The solutions involving partial transmissions and partial reflections are decomposed into a finite sum of solutions, obtained by solving problems involving only complete transmissions or complete reflections. For the latter class of problems, the method of Jin et al. (2005b) applies.
(2) A re-initialization technique is introduced such that waves coming from multiple transmissions and reflections can be combined seamlessly as new initial value problems. This is implemented by rewriting the sum of several delta functions as one delta measure with a suitable weight, which can easily be implemented numerically.

Semiclassical computation of quantum barriers
Correct modelling of electron transport in nano-structures, such as resonant tunnelling diodes, superlattices or quantum dots, requires the treatment of quantum phenomena in highly localized regions within the devices (so-called quantum wells), while the rest of the device can be dealt with by classical mechanics. However, solving the Schrödinger equation in the entire physical domain is usually too expensive, and thus it is attractive to use a multiscale approach as given in Ben Abdallah, Degond and Gamba (2002): that is, solve the quantum mechanics only in the quantum well, and couple the solution to classical mechanics outside the well. To this end, the following semiclassical approach for thin quantum barriers was proposed in Jin and Novak (2006).
Step 1. Solve the time-independent Schrödinger equation (either analytically or numerically) for the local barrier well to determine the scattering data, i.e., the transmission and reflection coefficients α T , α R .
Step 2. Solve the classical Liouville equation elsewhere, using the scattering data generated in Step 1 and the interface condition (10.1) given above.
The results for d = 1 and d = 2 given in Jin andNovak (2006, 2007) demonstrate the validity of this approach whenever the well is either very thin (i.e., of the order of only a few ε) or well separated. In higher dimensions, the interface condition (10.1) needs to be implemented in the direction normal to interface, and the interface condition may be non-local for diffusive transmissions or reflections (Jin and Novak 2007). This method correctly captures both the transmitted and the reflected quantum waves, and the results agree (in the sense of weak convergence) with the solution obtained by directly solving the Schrödinger equation with small ε. Since one obtains the quantum scattering information only in a preprocessing step (i.e., Step 1), the rest of the computation (Step 2) is classical, and thus the overall computational cost is the same as for computing classical mechanics. Nevertheless, purely quantum mechanical effects, such as tunnelling, can be captured.
If the interference needs to be accounted for, then such Liouville-based approaches are not appropriate. One attempt was made in Jin and Novak (2010) for one-dimensional problems, where a complex Liouville equation is used together with an interface condition using (complex-valued ) quantum scattering data obtained by solving the stationary Schrödinger equation. Its extension to multi-dimensional problems remains to be done. A more general approach could use Gaussian beam methods, which do capture the phase information. This is an interesting line of research yet to be pursued.

Schrödinger equations with matrix-valued potentials and surface hopping
Problems closely related to those mentioned in Section 10 arise in the study of semiclassical Schrödinger equations with matrix-valued potentials. This type of potential can be seen as a highly simplified model of the full manybody quantum dynamics of molecular dynamics. Using the celebrated Born-Oppenheimer approximation to decouple the dynamics of the electrons from the one for the much heavier nuclei (see, e.g., Spohn and Teufel (2001)), one finds that the electrons are subject to external forces which can be modelled by a system of Schrödinger equations for the nuclei along the electronic energy surfaces. The nucleonic Schrödinger system has matrixvalued potentials, which will be treated in this section.
To this end, we consider the following typical situation, namely, a timedependent Schrödinger equation with R 2×2 -matrix-valued potential (see, e.g., Spohn andTeufel (2001), or Teufel (2003)): with v 1 (x), v 2 (x) ∈ R. The matrix V then has two eigenvalues, Remark 11.1. In the Born-Oppenheimer approximation, the dimensionless semiclassical parameter ε > 0 is given by ε = m M , where m and M are the masses of an electron and a nucleus respectively (Spohn and Teufel 2001). Then, all oscillations are roughly characterized by the frequency 1/ε, which typically ranges between one hundred and one thousand.

Wigner matrices and the classical limit for matrix-valued potentials
In this section, we shall discuss the influence of matrix-valued potentials on the semiclassical limit of the Schrödinger equations (11.1). Introduce the Wigner matrix as defined in Gérard et al. (1997): We also let W denote the corresponding (weak) limit In order to describe the the dynamics of the limiting matrix-valued measure W (t, x, ξ), first note that the complex 2 × 2 matrix-valued symbol of (11.1) is given by The two eigenvalues of −iP (x, ξ) are These eigenvalues λ n , n = 1, 2 govern the time evolution of the limiting measure W (t), as proved in Gérard et al. (1997). They act as the correct classical Hamiltonian function in phase space, corresponding to the two energy levels, respectively. In the following, we shall say that two energy levels cross at a point x * ∈ R 2 if λ + (x * ) = λ − (x * ). Such a crossing is called conical if the vectors ∇ x v 1 (x * ) and ∇ x v 2 (x * ) are linearly independent. In Figure 11.1 we give two examples of conical crossings for potentials of the form If all the crossings are conical, the crossing set is a sub-manifold of co-dimension two in R 2 : see Hagedorn (1994). Assume that the Hamiltonian flows with Hamiltonians λ n leave invariant the set For (x, ξ) ∈ Ω, denote by χ n (x, ξ) the column eigenvector corresponding to the eigenvalue λ n (x, ξ), and the matrix is the orthogonal projection onto the eigenspace associated to λ n (x, ξ). By Theorem 6.1 of Gérard et al. (1997), the matrix-valued Wigner measure W (t) commutes with the projectors Π n , outside the crossing set S, and can thus be decomposed as Since the eigenspaces are both one-dimensional, the decomposition is simplified to be The scalar functions W n (t, x, ξ) given by are then found to satisfy the Liouville equations (11.4) subject to initial data W n (0) = tr(Π n W ), (x, ξ) ∈ Ω. (11.5) The scalar functions W n , n = 1, 2, are the phase space probability densities corresponding to the upper and lower energy levels, respectively. One can recover from them the particle densities ρ n via (11.6) In other words, the Liouville equations (11.4) yield the propagation of the Wigner measures W 1 (t, ·) and W 2 (t, ·) on any given time interval, provided that their supports do not intersect the eigenvalue crossing set S. However, analytical and computational challenges arise when their supports intersect the set S. In S the dynamics of W 1 , W 2 are coupled due to the non-adiabatic transitions between the two energy levels, and an additional hopping condition is needed (analogous to the interface condition considered in Section 10 above).

The Landau-Zener formula
Lasser, Swart and Teufel (2007) give a heuristic derivation of the nonadiabatic transition probability. The derivation is based on the Hamiltonian system corresponding to the Liouville equations (11.4), i.e., x n (t) = ∇ ξ λ n (t) = ξ n (t), ξ n (t) = −∇ x λ n (t), n = 1, 2. (11.7) The basic idea is to insert the trajectories (x(t), ξ(t)) of the Hamiltonian systems (11.7) into the trace-free part of the potential matrix (11.2) to obtain a system of ordinary differential equations, given by The non-adiabatic transitions occur in the region where the spectral gap between the eigenvalues becomes minimal. The function measures the gap between the eigenvalues in phase space along the classical trajectory (x(t), ξ(t)), where ϑ(x) = (v 1 (x), v 2 (x)) and | · | denotes the Euclidean norm. The necessary condition for a trajectory to attain the minimal gap is given by where ∇ x ϑ(x(t) is the Jacobian matrix of the vector ϑ(x(t)), and ξ(t) =ẋ(t).
Hence, a crossing manifold in phase space containing these points is given by The transition probability when one particle hits S * is assumed to be given by (11.8) which is the famous Landau-Zener formula (Landau 1932, Zener 1932. Note that T decays exponentially in x and ξ and In other words, as ε → 0, the transition between the energy bands only occurs on the set S * , which is consistent with the result in the previous subsections.

Numerical approaches
Lagrangian surface hopping One widely used numerical approach to simulation of the non-adiabatic dynamics at energy crossings is the surface hopping method, first proposed by Tully and Preston (1971), and further developed in Tully (1990) and Sholla and Tully (1998). The basic idea is to combine the classical transports of the system on the individual potential energy surfaces λ ± (x) that follow (11.7) with an instantaneous transitions at S * from one energy surface to another. The rate of transition is determined by the Landau-Zener formula (11.8) whenever available, or computed by some quantum mechanical simulation locally around S * . The hoppings were performed in a Monte Carlo procedure based on the transition rates. For a review of surface hopping methods see Drukker (1999). More recently, surface hopping methods have generated increasing interest among the mathematical community. For molecular dynamical simulations, Horenko, Salzmann, Schmidt and Schütte (2002) adopted the partial Wigner transform to reduce a full quantum dynamics into the quantum-classical Liouville equation, and then the surface hopping is realized by approximating the quantum Liouville equation using phase space Gaussian wave packets. From the analytical point of view, Lasser and Teufel (2005) and Fermanian Kammerer and Lasser (2003) analysed the propagation through conical surface crossings using matrix-valued Wigner measures and proposed a corresponding rigorous surface hopping method based on the semiclassical limit of the time-dependent Born-Oppenheimer approximation. In Lasser et al. (2007) and Kube, Lasser and Weber (2009) they used a particle method to solve the Liouville equation, in which each classical trajectory was subject to a deterministic branching (rather than Monte Carlo) process. Branching occurs whenever a trajectory attains one of its local minimal gaps between the eigenvalue surfaces. The new branches are consequently re-weighted according to the Landau-Zener formula for conical crossings These Lagrangian surface hopping methods are very simple to implement, and in particular, very efficient in high space dimensions. However, they require either many statistical samples in a Monte Carlo framework or the increase of particle numbers whenever hopping occurs. In addition, as is typical for Lagrangian methods, a complicated numerical re-interpolation procedure is needed whenever the particle trajectories diverge, in order to maintain uniform accuracy in time.

Eulerian surface hopping
The Eulerian framework introduced in Jin, Qi and Zhang (2011) consists of solving the two Liouville equations (11.4), with a hopping condition that numerically incorporates the . Note that the Schrödinger equation (11.1) implies conservation of the total mass, which in the semiclassical limit ε → 0, locally away from S * , yields (11.9) For this condition to hold for all x, ξ, the total flux in the direction normal to S * needs to be continuous across S * . To ensure this, the Landau-Zener transition at S * should be formulated as a continuity condition for the total flux in the normal direction e n . Define the flux function j n (x, ξ) ∈ R 2 for each eigenvalue surface via Assume, before hopping, that the particle remains on one of the eigenvalue surfaces, i.e., Then the interface condition is given by where (x ± , ξ ± ) are the (pre-and post-hopping) limits to (x 0 , ξ 0 ) ∈ S * along the direction e n .
Remark 11.2. There is a restriction of this approach based on the Liouville equation and the Landau-Zener transition probability. The interference effects generated when two particles from different energy levels arrive at S * at the same time are thereby not accounted for, and thus important quantum phenomena, such as the Berry phase, are missing. One would expect that Gaussian beam methods could handle this problem, but this remains to be explored.

Schrödinger equations with periodic potentials
So far we have only considered ε-independent potential V (x). The situation changes drastically if one allows for potentials varying on the fast scale y = x/ε. An important example concerns the case of highly oscillatory potentials V Γ (x/ε), which are periodic with respect to a periodic lattice Γ ⊂ R d . In the following we shall therefore consider Schrödinger equations of the form Here V ∈ C ∞ denotes some smooth, slowly varying potential and V Γ is a rapidly oscillating potential (not necessarily smooth). For definiteness we shall assume that for some orthonormal basis i.e., Γ = (2πZ) d .
Remark 12.1. Equations of this form arise in solid state physics, where they are used to describe the motion of electrons under the action of an external field and a periodic potential generated by the ionic cores. This problem has been extensively studied from both physical and mathematical points of view: see, e.g., Ashcroft and Mermin (1976), Bensoussan, Lions and Papanicolaou (1978), Teufel (2003) and the references given therein. One of the most striking dynamical effects due to the inclusion of a periodic potential V Γ is the occurrence of so-called Bloch oscillations. These are oscillations exhibited by electrons moving in a crystal lattice under the influence of a constant electric field V (x) = F · x, F ∈ R d (see Section 12.2 below).

Emergence of Bloch bands
In order to better understand the influence of V Γ , we recall here the basic spectral theory for periodic Schrödinger operators of the form (Reed and Simon 1976): With V Γ obeying (12.2), we have the following.
(1) The fundamental domain of the lattice Γ (2) The dual lattice Γ * can then be defined as the set of all wave numbers k ∈ R for which plane waves of the form e ik·x have the same periodicity as the potential V Γ . This yields Γ * = Z d in our case.
(3) The fundamental domain of the dual lattice, Y * , i.e., the (first) Brillouin zone, is the set of all k ∈ R closer to zero than to any other dual lattice point. In our case Y * = − 1 2 , 1 2 d , equipped with periodic boundary conditions, i.e., Y * T d .
By periodicity, it is sufficient to consider the operator H per on the fundamental domain Y only, where we impose the following quasi-periodic boundary conditions: It is well known (Wilcox 1978) that under mild conditions on V Γ , the operator H admits a complete set of eigenfunction {ψ m (y, k)} m∈N , parametrized by k ∈ Y * . These eigenfunctions provide, for each fixed k ∈ Y * , an orthonormal basis in L 2 (Y ). Correspondingly there exists a countable family of real eigenvalues {E(k)} m∈N , which can be ordered as where the respective multiplicities are accounted for in the ordering. The set {E m (k) : k ∈ Y } ⊂ R is called the mth energy band of the operator H per . Concerning the dependence on k ∈ Y * , it has been shown (Wilcox 1978) that for any m ∈ N there exists a closed subset X ⊂ Y * such that E m (k) is analytic and ψ m (·, k) can be chosen to be a real analytic function for all k ∈ Y * \X. Moreover, If this condition indeed holds for all k ∈ Y * , then E m (k) is called an isolated Bloch band. Moreover, it is known that meas X = meas {k ∈ Y * : E n (k) = E m (k), n = m} = 0.
(12.4) This set of Lebesgue measure zero consists of the so-called band crossings. See Figure 12.1 for an example of bands for Mathieu's model with potential V Γ = cos y.
Note that due to (12.3) one can rewrite ψ m (y, k) as for some 2π-periodic function χ m (·, k), usually called Bloch functions. In terms of χ m (y, k) the spectral problem for H per becomes (Bloch 1928) Most importantly, the spectral data obtained from (12.6) allow us to decompose the original Hilbert space L 2 (R d ) into a direct sum of so-called band spaces: L 2 (R d ) = ∞ m=1 H m . This is the well-known Bloch decomposition method , which implies that (12.8) The corresponding projection of f ∈ L 2 (R d ) onto the mth band space H m is thereby given via (Reed and Simon 1976) (12.9) In the following, we shall also denote (12.10) the coefficient of the Bloch decomposition.

Two-scale WKB approximation
Equipped with the basic theory of Bloch bands, we recall here an extension of the WKB method presented in Section 2.2 to the case of highly oscillatory periodic potentials. Indeed, it has been shown in Bensoussan et al. (1978) and Guillot, Ralston and Trubowitz (1988) that solutions to (12.1) can be approximated (at least locally in time) by where χ m is parametrized via k = ∇S(t, x). The phase function S thereby solves the semiclassical Hamilton-Jacobi equation The corresponding semiclassical flow X sc t : y → x(t, y) is given bẏ The wave vector k ∈ Y * is usually called crystal momentum. In the case of a constant electric field V = F · x, equation (12.13b) yields Note that since k ∈ Y * T d , this yields a periodic motion in time of x(t, y), the so-called Bloch oscillations.
In addition, the leading-order amplitude in (12.11) is found to be the solution of the semiclassical transport equation (Carles, Markowich and Sparber 2004), denotes the so-called Berry phase term (Carles et al. 2004, Panati, Spohn andTeufel 2006), which is found to be purely imaginary β m (t, x) ∈ (iR) d . It is, importantly, related to the quantum Hall effect (Sundaram and Niu 1999). The amplitude a is therefore necessarily complex-valued and exhibits a non-trivial phase modulation induced by the geometry of V Γ ; see also Shapere and Wilczek (1989) for more details. Note that (12.14) yields the following conservation law for ρ = |a| 2 : The outlined two-scale WKB approximation again faces the problem of caustics. Furthermore, there is an additional problem of possible bandcrossings at which ∇ k E m is no longer defined. The right-hand side of (12.11) can therefore only be regarded as a valid approximation for (possibly very) short times only. Nevertheless, it shows the influence of the periodic potential, which can be seen to introduce additional high-frequency oscillations (Γ-periodic) within u ε .
Remark 12.2. These techniques have also been successfully applied in weakly nonlinear situations (Carles et al. 2004).

Wigner measures in the periodic case
The theory of Wigner measures discussed in Section 3 can be extended to the case of highly oscillatory potentials. The theory of the so-called Wigner band series was developed in Markowich, Mauser and Poupaud (1994) and Gérard et al. (1997). The basic idea is to use Bloch's decomposition and replace the continuous momentum variable ξ ∈ R d by the crystal momentum k ∈ Y * .
A more general approach, based on space adiabatic perturbation theory (Teufel 2003), yields in the limit ε → 0 a semiclassical Liouville equation of the form

Bloch-decomposition-based time-splitting method
The introduction of a highly oscillatory potential V Γ x ε poses numerical challenges for the numerical computation of semiclassical Schrödinger equations. It has been observed by Gosse (2006) and Gosse and Markowich (2004) that conventional split-step algorithms do not perform well. More precisely, to guarantee convergence of the scheme, time steps of order O(ε) are required. In order to overcome this problem, a new time-splitting algorithm based on Bloch's decomposition method has been proposed in Huang, Jin, Markowich and Sparber (2007) and further developed in Huang, Jin, Markowich and Sparber (2008) and Huang, Jin, Markowich and Sparber (2009). The basic idea is as follows.
Step 1. For t ∈ [t n , t n+1 ] one first solves The main point is that by using the Bloch decomposition method, Step 1 can be solved exactly, i.e., only up to numerical errors. In fact, in each band space H m , equation (13.1) is equivalent to where u ε m ≡ P ε m u ε is the (appropriately ε-scaled) projection of u ε ∈ L 2 (R d ) onto H m defined in (12.9) and E m (−i∇) is the Fourier multiplier corresponding to the symbol E m (k). Using standard Fourier transformation, equation (13.2) is easily solved by where F −1 is the inverse Fourier transform. In other words, one can solve (13.1) by decomposing u ε into a sum of band-space functions u ε m , each of which is propagated in time via (13.3). After resummation this yields u ε (t n+1 , x). Once this is done, we proceed as usual to take into account V (x).
Step 2. On the same time interval as before, we solve the ODE where the solution obtained in Step 1 serves as initial condition for Step 2. In this algorithm, the dominant effects from the dispersion and the periodic lattice potential are computed in one step. It thereby maintains their strong interaction and treats the non-periodic potential as a perturbation. Because the split-step error between the periodic and non-periodic parts is relatively small, the time steps can be chosen considerably larger than a conventional time-splitting algorithm (see below). As in a conventional splitting method (see Section 5), the numerical scheme conserves the particle density ρ ε = |u ε | 2 on the fully discrete level. More importantly, if V (x) = 0, i.e., no external potential, the algorithm preserves the particle density (and hence the mass) in each individual band space H m .
Remark 13.1. Clearly, the algorithm given above is only first-order in time, but this can easily be improved by using the Strang splitting method (see Section 5). In this case, the method is unconditionally stable and comprises spectral convergence for the space discretization as well as secondorder convergence in time.

Numerical calculation of Bloch bands
In the numerical implementation of this algorithm a necessary prerequisite is the computation of Bloch bands E m (k) and Bloch eigenfunction χ m (y, k). This requires us to numerically solve the eigenvalue problem (12.6). In one space dimension d = 1 we proceed as in Gosse and Markowich (2004), by expanding V Γ ∈ C 1 (R) in its Fourier series Clearly, if V Γ ∈ C ∞ (R) the corresponding Fourier coefficients V (λ) decay faster than any power, as λ → ±∞, in which case we only need to take into account a few coefficients to achieve sufficient accuracy. Likewise, expand the Bloch eigenfunction χ m (·, k) in its respective Fourier series For λ ∈ {−Λ, . . . , Λ − 1} ⊂ Z, one consequently approximates the spectral problem (12.6) by the algebraic eigenvalue problem where the 2Λ × 2Λ matrix H(k) is given by The matrix H(k) has 2Λ eigenvalues. Clearly, this number has to be large enough to have sufficiently many eigenvalues E m (k) for the simulation, i.e., we require m ≤ 2Λ. Note, however, that the number Λ is independent of the spatial grid (in particular, independent of ε), thus the numerical costs of this eigenvalue problem are often negligible compared to those of the evolutionary algorithms (see below for more details). In higher dimensions d > 1, computing the eigenvalue problem (12.6) along these lines becomes numerically too expensive to be feasible. In many physical applications, however, the periodic potential splits into a sum of one-dimensional potentials, i.e., where y = (y 1 , y 2 , . . . , y d ) ∈ R d . In this case, Bloch's spectral problem can be treated separately (using a fractional step-splitting approach) for each coordinate y j ∈ R, as outlined before.
Remark 13.2. In practical applications, the accurate numerical computation of Bloch bands is a highly non-trivial task. Today, though, there exists a huge amount of numerical data detailing the energy band structure of the most important materials used in the design of semiconductor devices, for example. In the context of photonic crystals the situation is similar. Thus, relying on such data one can in principle completely avoid the above eigenvalue computations and their generalizations to more dimensions. To this end, one should also note that, given the energy bands E m (k), we do not need any knowledge about V Γ to solve (12.1) numerically. Also, we remark that it was shown in Huang et al. (2009) that the Bloch-decomposition-based time-splitting method is remarkably stable with respect to perturbations of the spectral data.

Implementation of the Bloch-decomposition-based time-splitting method
In the numerical implementation we shall assume that V Γ admits the decomposition (13.1). In this case we can solve (12.1) by using a fractional step method, treating each spatial direction separately, i.e., one only needs to study the one-dimensional equation on the time interval [t n , t n+1 ]. This equation will be considered on a one-dimensional computational domain (a, b) ⊂ R, equipped with periodic boundary conditions (necessary to invoke fast Fourier transforms). We suppose that there are L ∈ N lattice cells within (a, b) and numerically compute u ε at L × R grid points, for some R ∈ N. In other words we assume that there are R grid points in each lattice cell, which yields the discretization and thus u n ≡ u(t n ) are evaluated at the grid points x ,r = ε(2π( − 1) + y r ). (13.8) Note that in numerical computations one can use R L whenever ε 1, i.e., only a few grid points are needed within each cell.
Keeping in mind the basic idea of using Bloch's decomposition, one has the problem that the solution u ε of (13.6) does not in general have the same periodicity properties as ϕ m . A direct decomposition of u ε in terms of this new basis of eigenfunctions is therefore not possible. This problem can be overcome by invoking the following unitary transformation for f ∈ L 2 (R): with the properties f (y + 2π, k) = e 2 iπk f (y, k), f (y, k + 1) = f (y, k).
In other words f (y, k) admits the same periodicity properties with respect to k and y as the eigenfunction ψ m (y, k). In addition, the following inversion formula holds: (13.9) Moreover, one easily sees that the Bloch coefficient, defined in (12.10), can be equivalently written as (13.10) which, in view of (12.5), resembles a Fourier integral. In fact, all of these formulae can be easily implemented by using the fast Fourier transform.
The numerical algorithm needed to perform Step 1 outlined above is then as follows.
Step 1.1. First compute u ε at time t n by where x ,r is as in (13.8).
Step 1.2. Compute the coefficient C ε m (t n , k ) via (13.10): Step 1.3. Evolve C ε m (t n , k) up to time t n+1 according to (13.2), Step 1.4. u ε can be obtained at time t n+1 by summing over all band contributions, Step 1.5. Perform the inverse transformation (13.9): This concludes the numerical procedure performed in Step 1.
The Bloch-decomposition-based time-splitting method was found to be converging for ∆x = O(ε) and ∆t = O(1); see Huang et al. (2007) for more details. In other words, the time steps can be chosen independently of ε, a huge advantage in comparison to the standard time-splitting method used in Gosse (2006), for instance. Moreover, the numerical experiments done in Huang et al. (2007) show that, of only a few Bloch bands, E m (k) are sufficient to achieve very high accuracy, even in cases where V (x) is no longer smooth (typically m = 1, . . . , M, with M ≈ 8 is sufficient). Applications of this method are found in the simulation of lattice Bose-Einstein condensates (Huang et al. 2008) and of wave propagation in (disordered) crystals (Huang et al. 2009).
Remark 13.3. For completeness, we recall the numerical complexities for the algorithm outlined above: see Huang et al. (2007). The complexities of Steps 1.1 and 1.5 are O(RL ln L), the complexities of Steps 1.2 and 1.4 are O(MLR ln R), and for Step 1.3 it is O(ML). The complexity of the eigenvalue problem (13.5) is O(Λ 3 ). However, since Λ (or R) is independent of ε and since (13.5) needs to be solved only once (as a preparatory step), the computational costs for this step are negligible. In addition, since M and R are independent of ε, one can choose R L and M L, whenever ε 1. Finally, one should notice that the complexities in each time step are comparable to the usual time-splitting method.

Moment closure in Bloch bands
It is straightforward to adapt the moment closure method presented in Section 6 to the case of periodic potentials. To this end, one considers the semiclassical Liouville equation (12.16), i.e., and close it with the following ansatz for the Wigner measure: where we let δ denote the Γ * -periodic delta distribution, i.e., By following this idea, Gosse and Markowich (2004) showed the applicability of the moment closure method in the case of periodic potentials in d = 1 (see also Gosse (2006)). In addition, self-consistent Schrödinger-Poisson systems were treated in Gosse and Mauser (2006). As mentioned earlier, extending this method to higher space dimensions d > 1 is numerically challenging.

Gaussian beams in Bloch bands
The Gaussian beam approximation, discussed in Sections 8 and 9, can also be extended to the Schrödinger equation with periodic potentials. To this end, one adopts the Gaussian beam within each Bloch band of the Schrödinger equation (12.1). In the following, we shall restrict ourselves to the case d = 1, for simplicity.

Lagrangian formulation
As for the two-scale WKB ansatz (12.11), we define where the mth band Liouville operator L m is defined as (13.12) and Φ m is the complex-valued d-dimensional level set function for the velocity corresponding to the mth Bloch band. They also solve the following inhomogeneous Liouville equations for the phase S m and amplitude a m : where L m is defined by (13.12) and β m denotes the Berry phase term as given by (12.15). Here, the Hessian M m ∈ C is obtained from

Schrödinger equations with random potentials
Finally, we shall consider (small) random perturbations of the potential V (x). It is well known that in one space dimension, linear waves in a random medium get localized even when the random perturbations are small: see, e.g., Fouque, Garnier, Papanicolaou and Sølna (2007). Thus the analysis here is restricted to three dimensions. (The two-dimensional case is difficult because of criticality, i.e., the mean-field approximation outlined below is most likely incorrect.)

Scaling and asymptotic limit
Consider the Schrödinger equation with a random potential V R : Here V R (y) is a mean zero, stationary random function with correlation length of order one. Its correlation length is assumed to be of the same order as the wavelength. The √ ε-scaling given above is critical in the sense that the influence of the random potential is of the same order as the one given by V (x) (see also the remark below). We shall also assume that the fluctuations are statistically homogeneous and isotropic so that where · · denotes statistical average and R(|x|) is the covariance of random fluctuations. The power spectrum of the fluctuations is defined by When (14.2) holds, the fluctuations are isotropic and R is a function of |ξ| only.
Remark 14.1. Because of the statistical homogeneity, the Fourier transform of the random potential V R is a generalized random process with orthogonal increments If the amplitude of these fluctuations is large, then purely random scattering will dominate and waves will be localized: see Fröhlich and Spencer (1983).
On the other hand if the random fluctuations are too weak, they will not affect the transport of waves at all. Thus, in order to have scattering induced by the random potential and the influence of the slowly varying background V (x) affect the (energy transport of the) waves in comparable ways, the fluctuations in the random potential must be of order √ ε.
Using the ε-scaled Wigner transformation, we can derive the analogue of (3.2) in the following form: where the pseudo-differential operator Θ ε is given by (3.3). The behaviour of this operator as ε → 0 is very different from the case without V R , as can already be seen on the level of formal multiscale analysis: see Ryzhik, Papanicolaou and Keller (1996). Let y = x/ε be a fast spatial variable, and introduce an expansion of w ε (t) in the following form: w ε (t, x, ξ) = w(t, x, ξ) + ε 1/2 w (1) (t, x, y, ξ) + εw (2) (t, x, y, ξ) + · · · . (14.6) Note that we hereby assume that the leading term w does not depend on the fast scale. We shall also assume that the initial Wigner distribution w ε in (x, ξ) tends to a smooth, non-negative function w in (x, ξ), which decays sufficiently fast at infinity. Then, as ε → 0, one formally finds that w ε (t) , i.e., the averaged solution to (14.5), is close to a limiting measure w(t), which satisfies the following linear Boltzmann-type transport equation: Here, the linear scattering operator Q is given by ξ)) dp, (14.7) with differential scattering cross-section σ(k, p) = 4π R(ξ − p)δ(|ξ| 2 − |p| 2 ) (14.8) and total scattering cross-section Note that the transport equation (14.7) has two important properties. First, the total mass (or energy, depending on the physical interpretation) is conserved, i.e., Second, the positivity of w(t, x, ξ) is preserved. Rigorous mathematical results concerning the passage from (14.1) to the transport equation (14.7) can be found in Spohn (1977), Dell'Antonio (1983) and Ho, Landau and Wilkins (1993) (which contains extensive references), and, more recently, in Erdös and Yau (2000).

Coupling with other media
One can also study the problems when there are other media, including periodic media and interfaces (flat or random).
In the case of periodic media coupled with random media, one can use the above multiscale analysis combined with the Bloch decomposition method to derive a system of radiative transport equations: see Bal, Fannjiang, Papanicolaou and Ryzhik (1999a). In the limit system ε → 0, the transport part is governed by the semiclassical Liouville equation (12.16) in each Bloch band, while the right-hand side is a non-local scattering operator (similar to the one above) coupling all bands.
One can also consider random high-frequency waves propagating through a random interface. Away from the interface, one obtains the radiative transport equation (14.7). At the interface, due to the randomness of the interface, one needs to consider a diffusive transmission or reflection process, in which waves can be scattered in all directions (see Figure 10.1). To this end, a non-local interface condition has to be derived: see Bal et al. (1999b).
So far there have been few numerical works on random Schrödinger equations of the form (14.1). Bal and Pinaud (2006) studied the accuracy of the radiative transport equation (14.7) as an approximation to (14.1).  discretized a non-local interface condition for diffusive scattering, similar in spirit to the treatment described in Section 10.2. We also mention that the temporal resolution issue for time-splitting approximations of the Liouville equations with random potentials was rigorously studied in Bal and Ryzhik (2004). Finally, we refer to Fouque et al. (2007) and Bal, Komorowski and Ryzhik (2010) for a comprehensive reading on high-frequency waves in random media.

Nonlinear Schrödinger equations in the semiclassical regime
So far we have only considered linear Schrödinger equations. Nonlinear models, however, are almost as important, since they describe a large number of physical phenomena in nonlinear optics, quantum superfluids, plasma physics or water waves: see, e.g., Sulem and Sulem (1999) for a general overview. The inclusion of nonlinear effects poses new challenges for mathematical and numerical study.
The first case of nonlinearities has important applications in Bose-Einstein condensation and nonlinear optics, whereas the latter typically appears as a mean-field description for the dynamics of quantum particles, say electrons (in which case one usually has V 0 (x) = ± 1 |x| in d = 3). Concerning the existence of solutions to (15.1), we shall in the following only sketch the basic ideas. For more details we refer to Cazenave (2003) and Carles (2008). As a first step we represent the solution to (15.1) using Duhamel's formula: where U ε (t) is the unitary semi-group generated by the linear Hamiltonian H = − 1 2 ∆ + V (see the discussion in Section 2.1). The basic idea is to prove that the right-hand side of (15.2) defines a contraction mapping (in some suitable topology). To this end, one has to distinguish between the case of sub-and supercritical nonlinearities. The existence of a unique solution solutions (in the supercritical regime). Rather, we have to expect finitetime blow-up to occur in general. Clearly, this phenomenon will also have a significant impact on numerical simulations, in particular for ε 1.

WKB analysis of nonlinear Schrödinger equations
In order to understand the influence of nonlinear terms in the limit ε → 0, one can again invoke a WKB approximation of the form (2.10): ∼ a ε (t, x) e iS(t,x)/ε .
Let us assume for simplicity that α ∈ N. Then, plugging the ansatz given above into (15.1), and ordering equal powers in ε, we see that S solves We see that in the case α = 0 we can no longer solve the Hamilton-Jacobi equation for the phase S independently of the amplitude a. In other words, the amplitude influences the geometry of the rays or characteristics. This is usually referred to as supercritical geometric optics (Carles 2008), not to be confused with the supercritical regime concerning the existence of solutions, as outlined in Section 15.1 above.

Weakly nonlinear geometric optics
In contrast, the situation for α > 0 (subcritical geometric optics) is similar to the linear situation, in the sense that the rays of geometric optics are still given by (2.13) and thus independent of the nonlinearity. In this case, the method of characteristics guarantees the existence of a smooth S ∈ C ∞ ([−T, T ] × R d ) up to caustics, and one can proceed with the asymptotic expansion to obtain the following transport equation for the leading-order amplitude: ∂ t a + ∇S · ∇a + a 2 ∆S = 0 i fα > 1, −if (x, |a| 2 ) if α = 1. (15.5) One sees that if α > 1, nonlinear effects are asymptotically negligible for ε 1. The solution is therefore expected to be roughly the same as in the linear situation (at least before caustics). For α = 1, however, nonlinear effects show up in the leading-order amplitude. Note, however, that by multiplying (15.5) withā and taking the real part, one again finds ∂ t ρ + div(ρ∇S) = 0, which is the same conservation law for ρ = |a| 2 as in the linear case. The nonlinear effects for α = 1 are therefore given solely by nonlinear phase modulations of the leading-order amplitude. In the case of a simple cubic nonlinearity, one explicitly finds (Carles 2000) a(t, x) = a 0 (y(t, x)) J t (y(t, x)) e iG(t,x) , |t| ≤ T, (15.6) where the slowly varying phase G is given by This regime is therefore often called weakly nonlinear geometric optics, and indeed it is possible to prove a rigorous approximation result analogous to Theorem 2.3 in this case too.

Supercritical geometric optics
On the other hand, the situation for α = 0 is much more involved. Indeed, it can be easily seen that a naive WKB approximation breaks down since the system of amplitude equation is never closed, i.e., the equation for a n , obtained at O(ε n ), involves a n+1 , and so on (a problem which is reminiscent of the moment closure problem discussed in Section 6 above). This difficulty was overcome by Grenier (1998) and Carles (2007aCarles ( , 2007b, who noticed that there exists an exact representation of u ε in the form u ε (t, x) = a ε (t, x) e iS ε (t,x)/ε , (15.7) with real-valued phase S ε (possibly ε-dependent) and complex-valued amplitude a ε . (Essentially, the right-hand side of (15.7) introduces an additional degree of freedom by invoking complex amplitudes.) Plugging (15.7) into (15.1) with α = 0, one has the freedom to split the resulting equations for a ε and S ε as follows. On plugging the ansatz given above into (15.1), and ordering equal powers in ε, one sees that S solves Formally, we expect the limit of this system as ε → 0 to give a semiclassical approximation of u ε , at least locally in time. Indeed this can be done by first rewriting these equations into a new system for ρ ε = |a ε | 2 and v ε = ∇S ε . Under some assumptions on the nonlinearity f (satisfied, for instance, in the cubic case), the obtained equations form a strictly hyperbolic system in which one can rigorously pass to the limit as ε → 0 to obtain ∂ t v + v · ∇v + V (x) + f (ρ) = 0, ∂ t ρ + div(ρv) = 0. (15.8) This is a system of (irrotational) Euler equations for a classical fluid with pressure law p(ρ) = f (ρ)/ρ. Following this approach, one can prove that as long as (15.8) admits (local-in-time) smooth solutions ρ, v, the quantum mechanical densities (2.6)-(2.7) indeed converge strongly in C ([0, T ]; L 1 (R d )), with ρ ε ε→0 −→ ρ, J ε ε→0 −→ ρ.
Reconstructing from these limits an approximate solution of the nonlinear Schrödinger equation in WKB form, i.e., u ε app (t, x) = ρ(t, x) e iS(t,x)/ε , v(t, x) := ∇S(t, x), generally requires some care (see Carles (2007b) for more details). Essentially one needs to take into account a higher-order corrector to ensure that the formal approximation is indeed correct up to errors of order O(ε), independent of t ∈ [0, T ].
The semiclassical limit for α = 0 for |t| > T , i.e., beyond the formation of shocks in (15.8), is a very challenging mathematical problem. In the one-dimensional case the only available result is given by Jin, Levermore and McLaughlin (1999), using the inverse scattering technique. The semiclassical limit of focusing NLS is more delicate, but can also be obtained by inverse scattering: see Kamvissis, McLaughlin and Miller (2003).
Remark 15.3. Let us close this subsection by noticing that the analysis for Hartree-type nonlinearities is found to require slightly less sophistication than for local ones (Alazard and Carles 2007) and that the intermediate case 0 < α < 1 can be understood as a perturbation of the situation for α = 0 (Carles 2007b). In addition, a threshold condition for global smooth solutions to (15.8) has been determined in Liu and Tadmor (2002).

Wigner measure techniques for nonlinear potentials
One might hope to extend the results for Wigner measures given in Section 3 to nonlinear problems. This would have the considerable advantage of avoiding problems due to caustics. Unfortunately, this idea has not been very successful so far, the reason being that Wigner measures are obtained as weak limits only, which in general is not sufficient to pass to the limit in nonlinear terms. Indeed, one can even prove an ill-posedness result for Wigner measures in the nonlinear case (Carles 2001).
A notable exception is the case of Hartree nonlinearities f = V 0 * |u ε | with smooth interaction kernels V 0 ∈ C 1 b (R d ) (Lions and Paul 1993). In this case the Wigner measure associated to u ε is found to be a solution of the self-consistent Vlasov equation: The physically more interesting case of non-smooth interaction kernels, such as V 0 ∼ 1 |x| , which describes the coupling to a Poisson equation, is not covered by this result and can only be established in the particular case of the fully mixed quantum state; see Lions and Paul (1993) and Markowich and Mauser (1993) for more details, and also Zhang (2002) for a similar result valid for short times only.

Numerical challenges
Due to the introduction of nonlinear effects, the numerical difficulties discussed in Sections 4 and 5 are enhanced. The main numerical obstacles are the formation of singularities in focusing nonlinear Schrödinger equations and the creation of new scales at caustics for both focusing and defocusing nonlinearities. Basic numerical studies of semiclassical NLS were conducted in Bao et al. (2002) and Bao, Jin and Markowich (2003b), with the result that finite difference methods typically require prohibitively fine meshes to even approximate observables well in semiclassical defocusing and focusing NLS computations. In the case when these very restrictive meshing constraints are bypassed, the usual finite difference schemes for NLS can deliver incorrect approximations in the classical limit ε → 0, without any particular sign of instability (Carles and Gosse 2007). Time-splitting spectral schemes are therefore the preferred method. To this end, we refer to Jin, Markowich and Zheng (2004) for the application of the time-splitting spectral method to the Zakharov system, to Huang, Jin, Markowich, Sparber and Zheng (2005) for the numerical solution of the Dirac-Maxwell system, and to Bao et al. (2003b) for numerical studies of nonlinear Schrödinger equations. In addition, let us mention Markowich (2003a, 2004), where numerical simulations of the cubically nonlinear Gross-Pitaevskii equation (appearing in the description of Bose-Einstein condensates) are given using time-splitting trigonometric spectral methods. A numerical study of ground state solutions of the Gross-Pitaevskii equation can be found in Bao, Wang and Markowich (2005).
Due to the nonlinear creation of new highly oscillatory scales in the limit ε → 0, time-splitting methods suffer from more severe meshing restrictions for NLS than for linear Schrödinger equations, particularly after the appearance of the first caustics in the corresponding WKB approximation; see Bao et al. (2003b) and Carles and Gosse (2007) for more details. In the weakly nonlinear regime the following meshing strategy is sufficient: ∆x = O(ε), ∆t = O(ε) (to be compared with (5.5)). In the regime of supercritical geometric optics, however, one typically requires (even for quadratic observable densities) i.e., a severe restriction on the time step. In addition, one may need to invoke the Krasny filter (Krasny 1986), i.e., high Fourier-mode cut-offs, to avoid artifacts (such as symmetry breaking) in focusing NLS computations (Bao et al. 2003b). The latter, however, violates the conservation of mass: a clear drawback from the physics point of view. In order to overcome this problem, higher-order methods (in time) have to be deployed, such as exponential time-differencing or the use of integrating factors, and we refer to Klein (2008) for a comparison of different fourth-order methods for cubic NLS in the semiclassical regime.
Remark 15.5. In the closely related problem of the complex Ginzburg-Landau equation in the large space and time limit, the situation is known to be slightly better, due to the dissipative nature of the equation; see Degond, Jin and Tang (2008) for a numerical investigation. Finally, we note that the cubic NLS in d = 1 is known to be fully integrable by means of inverse scattering. This feature can be used in the design of numerical algorithms, as has been done in Zheng (2006), for instance. A generalization to higher dimensions or more general nonlinearities seems to be out of reach at present.
The lack of a clear mathematical understanding of the asymptotic behaviour of solutions to the semiclassical NLS beyond the formation of caustics has so far hindered the design of reliable asymptotic schemes. One of the few exceptions is the case of the Schrödinger-Poisson equation in d = 1, which can be analysed using Wigner measures, and which has recently been studied numerically in Jin, Wu and Yang (2010a) using a Gaussian beam method. In addition, moment closure methods have been employed for this type of nonlinearity, since it is known that the underlying classical problem, i.e., the Euler-Poisson system, allows for a construction of multivalued solutions. Numerical simulations for the classical system have been conducted in Wöhlbier et al. (2005). In addition, the case of the Schrödinger-Poisson equation with periodic potential is treated in Gosse and Mauser (2006).