Semi-classical Models for the Schrödinger Equation with Periodic Potentials and Band Crossings

The Bloch decomposition plays a fundamental role in the study of quantum mechanics and wave propagation in periodic media. Most of the homogenization theory developed for the study of high frequency or semi-classical limit for these problems assumes no crossing of the Bloch bands, resulting in a classical Liouville equation in the limit along each Bloch band. 
   
 In this article, we derive semi-classical models for the Schrodinger equation in periodic media that take into account band crossings, which is important to describe quantum transitions between Bloch bands. Our idea is still based on the Wigner transform (on the Bloch eigenfunctions), but in taking the semi-classical approximation, we retain the off-diagonal entries of the Wigner matrix, which cannot be ignored near the points of band crossings. This results in coupled inhomogeneous Liouville systems that can suitably describe quantum tunneling between bands that are not well-separated. We also develop a domain decomposition method that couples these semi-classical models with the classical Liouville equations (valid away from zones of band crossings) for a multiscale computation. Solutions of these models are numerically compared with those of the Schrodinger equation to justify the validity of these new models for band-crossings.

Appendix C The integration method of a simple model system 32

Abstract
The Bloch decomposition plays a fundamental role in the study of quantum mechanics and wave propagation in periodic media. Most of the homogenization theory developed for the study of high frequency or semiclassical limit for these problems assumes no crossing of the Bloch bands, resulting in classical Liouville equations in the limit along each Bloch band.
In this article, we derive semiclassical models for the Schrödinger equation in periodic media that take into account band crossing, which is important to describe quantum transitions between Bloch bands. Our idea is still based on the Wigner transform (on the Bloch eigenfunctions), but in taking the semiclassical approximation, we retain the off-diagonal entries of the Wigner matrix, which cannot be ignored near the point of band crossing. This results in coupled inhomogenious Liouville systems that can suitably describe quantum tunnelling between bands that are not well-separated. We also develop a domain decomposition method that couples these semiclassical models with the classical Liouville equations (valid away from zones of band crossing) for a multiscale computation. Solutions of these models are numerically compared with those of the Schröding equation to justify the validity of these new models for band-crossings.

Introduction
The linear Schrödinger equation with a periodic potential is an important model in solid state physics. It describes the motion of electrons in a crystal with a lattice structure. We consider the following one dimensional Schrödinger equation Here, φ ε is the complex-valued wave function, ε is the dimensionless rescaled Planck constant, V (x) is a periodic potential with the lattice L = 2πZ, so V (x + ν) = V (x) for any x ∈ R and ν ∈ L. U (x) is a smooth external potential function. We denote the dual lattice of L by L * and L * = Z. The fundamental domain of L is (0, 2π) and the first Brillouin zone is B = − 1 2 , 1 2 . In the semi-classical regime ε ≪ 1, the Schrödinger equation has a semi-classical limit governed by the Liouville system [31], which was justified rigorously in [10,21,22]. It has been shown that the electrons remain in a certain quantum subsystem, "move along the m-th band" and the dynamics is given byẋ = ∂ k E m (k), where E m is the energy corresponding to the m-th Bloch band [2]. Higher order correction relevant to Berry phase can be included, see [24,25,8]. All of these results use the adiabatic assunption, namely, different Bloch bands are well-separated and there is no band-crossing.
The non-adiabatic or diabatic effect should be considered whenever the transitions between energy levels of the quantum system play an important role. This may happen when the gap between the energy levels becomes small enough in comparison to the scaled Planck constant ε. The well-known Landau-Zener formula [18,32] describes the asymptotic effect of avoided crossings in various specific situations. The study of such "qunatum tunnellings" is important in many applications, from quantum dynamics in chemical reaction [30], semiconductors to Bose-Einstein condensation [3]. While in the case of band separation there have been significant mathematical progress in understanding the semiclassical limit [1,10,22,24,25], as well as numerical methods that utilizes the Bloch decomposition [14], there has been little mathematical and computational works for the band-mixing case. In the context of "surface hopping method", associated with the Born-Oppenheimer approximation, the Landau-Zener phenomenon has been studied computationally by Tully etc. [30,29,26,7] and mathematically [12,9,20,19]. For a Liouville equation based computational approach, see [16], and a quantum-classical model for surface hopping, see [23,13].
In this paper, we use the Wigner-Bloch theory to derive semiclassical models for the Schrödinger equation with periodic potentials (1.1) that account for band-crossing. Without band-crossing, the classical Liouville equation can be obtained along each Bloch band [1,10,22], as ε → 0. This only includes the diagonal entries of the Wigner-Bloch matrix (the Wigner transform of the Bloch functions) which is valid away from the crossing zone.
In our semiclassical models we include the leading order of the off-diagonal entries, resulting in system of complex, inhomogeneous Liouville equations. These systems contain terms that describe for transitions between bands, as well as Berry phase information which is related to the quantum Hall effect [28]. These systems are still hyperbolic, with oscillatory forcing terms that converge (in the weak sense) to zero in the semiclassical limit away from band-crossing zones so the classical Liouville equation can be recovered. Several numerical experiments using these models produce numerical solutions that adequately describe the quantum transitions between bands, when compared with the solutions of the original Schrödinger equation (1.1).
The computational cost of these semiclassical models, while considerably lower than that of the original Schröginer equation, is higher than the adiabatic Liouville equations without band-crossing. In order to further reduce the computational cost, a hybrid method that couples the classical Liouville equation away from the crossing zone to the new semiclassical models to be used in the crossing zones is introduced. A couple of other more efficient computational approaches for linear periodic potentials, are also introduced in section 5. This paper is organized as follows. In Section 2, we introduce the Wigner transform and the Bloch eigenvalue problem that are the two tools to study the semiclassical limit of the Schrödinger equation with periodic potentials. We also review the classical limit without band-crossing which gives rise to the classical Liouville equation for each Bloch band.
In Section 3 we derive the new semi-classical Liouville systems that account for quantum transitions between different bands. In Section 4 we introduce a domain-decomposition based hybrid model that couples the classical Liouville systme away from the band crossing zone to the new semiclassical models used in the crossing zones. Some numerical examples are presented in Section 5 to validate these semiclassical models for quantum transions between Bloch bands. More efficient numerical methods were introduced for simpler linear potentials. Section 6 concludes this paper.
2 The classical limit without band-crossing

The Wigner transform
Define the asymmetric Wigner transformation as in [1], whereφ ε is the complex conjugate of φ ε , and φ ε is the solution to the Schrödinger equation (1.1), then one obtains the Wigner equation: whereÛ (ω) is the Fourier transform of U (y): Denote z = x ε as the fast variable. To separate the dependence on both the slow and the fast variables, one can write W ε (t, x, k) as W ε (t, x, z, k), and replace the spatial derivative ∂ ∂x by ∂ ∂x + 1 ε ∂ ∂z . Then (2.2) becomes: Assume U (x) is smooth enough, one throws away high order terms, and (2.5) becomes: where the skew symmetric operator L is given by Asymptotically one expands W ε by: and plugs this ansatz into (2.5). In the order of O( 1 ε ) and O(1) respectively, one gets: This implies that W 0 is in the kernel of L. We seek a good basis of kerL in section 2.2, and leave the exact formulation of W 0 to section 2.3.

The Bloch eigenvalue problem
The eigenfunctions of L are constructed by studying the following eigenvalue problem: With each specific p, one constructs a boundary condition (2.9b,2.9c) and solves the eigenvalue problem. Denote E m (p) as the m-th eigenvalue with multiplicity r m , and Ψ α m , with α = 1, ..., r m , as the associated α-th eigenfunction. Ψ m is usually called the m-th Bloch eigenfunction [2]. For this problem, one can easily check the following properties: (a) The eigenvalues E m (p) are L * −periodic in p and have constant finite multiplicity outside a closed zero-measure subset F of p ∈ R. Outside F , one orders the eigenvalues (c) For all φ ∈ L 2 (R), one has the following Bloch decomposition: c α m (p)Ψ α m (x, p) dp (2.11) where c α m is the Bloch coefficient: then Φ α m is a z−periodic function with period 2π. For any p ∈ − 1 2 , 1 2 , {Φ α m (·, p)} forms a complete orthonormal basis in L 2 (0, 2π).
For simplicity, throughout this paper we make the following assumption: Assumption 1. The multiplicity for each eigenvalue E m is 1 outside F . We will drop the superscript α when there is no ambiguity.

The classical limit without band-crossing
In this section, we review the classical limit without band-crossing [1]. For that, we need the following assumption: The set F is a null set, i.e. eigenvalues are strictly apart from each other By taking the Wigner transformation on the Bloch eigenfunctions, one obtains a basis on the phase space. Define the z−periodic functions Q mn (z, k) by: where k is an arbitrary real number and is decomposed as: (2.14) Lemma 2.1. Define the inner product ·, · : then for any p ∈ (−1/2, 1/2), {Q mn (·, ·, p)} forms a complete orthonormal basis in L 2 (0, 2π)× Proof. The orthonormal condition can be proved by simply using (2.10). To prove the completeness, it is sufficient to show that: if there exists an f ∈ L 2 (0, 2π) × ℓ 2 (L * ) , such that f, Q mn = 0 for all m, n ∈ N, Since {Φ n (·, p)} forms a complete orthonormal basis in L 2 (0, 2π), the above equality implies which implies that f (z, µ) ≡ 0.
A straightforward computation gives: Apparently, under Assumptions 1 and 2, (2.16) gives: with σ mm representing the expansion coefficients. To derive the equation for them, one plugs it back into (2.8b), and takes the inner product with Q mm on both sides, the right hand side vanishes due to the skew symmetry of L and (2.16): The left hand side, on the other hand, gives: Thus one can combine (2.19) and (2.20), and obtain: In the derivation of (2.20), the following equalities were used : The details can also be found in [1].
Notice that 2.21 is the classical Liouville equation for each Bloch band.

Asymptotic models for the band-to-band transition
This section is for the case when Assumption 2 is not satisfied. Physically, ε is small but nonzero, and electrons can tunnel across bands for all p. But this tunneling is negligible when bands are far away from each other, the so-called adiabatic assumption. This is the basis for the asymptotically expansion (2.18), which throws away the dependence on the bandtransition terms Q mn (m = n) of W 0 . However, if there exists m 0 = n 0 and a point p c , such , so the asymptotic expansion (2.8) does not hold any more. When this happens, physically one observes significant tunneling effect. We seek asymptotic models to handle the band-to-band transition in this section.

The derivation of a two-band semi-classical Liouville system using the asymmetric Wigner transform
In order to handle the band-to-band transition phenomena, we come back to the asymptotic model (2.6) and use the following expression (compare with (2.18)!): Without loss of generality, we tackle a two-band problem. Define p c = arg min p {|E 1 − E 2 |} and assume |E 1 (p c ) − E 2 (p c )| = O(ε) and p c = 0.
Plugging (3.1) into (2.6) and taking the inner product with Q mn as mentioned in Sec.
2.3, one gets a system for σ mn 's: This system can be written in vector form as: Where the superscript T denotes the matrix transpose, I is the identity matrix. Generally, ψ 12 =ψ 21 , φ 12 = −φ 21 are complex-valued quantities, and φ 11 , φ 22 are purely imaginary, so A = A † is Hermitian, and C = −C † is anti-Hermitian (where the superscript † denotes the matrix conjugate transpose). Note that ∂ x U φ mm is the so-called Berry-phase [4,25,28].
It can be shown the system (3.2) is a hyperbolic system (see the details in Appendix A). Since the source matrix C = −C † , all of its eigenvalues are purely imaginary. Thus this system has unique bounded solution if the initial data are bounded [6] uniformly in varepsilon.
For later convenience, we call the semi-classical Liouville system (3.2) Liouville-A system, and this is the system that will be discussed and numerically solved in the paper.
One also needs to equip it with appropriate initial condition. Choose the initial data of the Schrödinger equation as two wave packets along the two Bloch bands in the following form [4,14]: Then, the initial data of the Wigner function, for ε ≪ 1, has the approximation: then by ignoring the O(ε) term, and using the periodicity of Φ m (z, p) on z, one can change the integral into a summation of integrals from 0 to 2π: Applying the equality In the above derivation, from the second line to the third line, we use the fact that Without loss of generality, we assume that ∂ x S 0 ∈ (−1/2, 1/2), then (3.5) becomes Compare (3.6) with (3.4), one has the initial data for σ: 3.2 The semi-classical Liouville system using the symmetric Wigner transform All the analysis above is done in the framework of the asymmetric Wigner transform (2.1). One could also use the symmetric Wigner transform: The derivation is similar, thus we skip the details, and give a list of the results:

The Wigner equation corresponding (2.2) is
2. Corresponding to the asymptotic Wigner equation for asymmetrical transformation (2.6), one has: where the skew symmetric operator L s is given by 3. Corresponding to (2.8), one has the O( 1 ε ) and O(1) expansions: 4. Same as in (2.13), one has the following symmetrical definition for Q mm (3.14) By taking the inner product with Q s mm on both sides of (3.11b), one obtains the same classical Liouville equations for σ s mm as in (2.21) 6. If some bands touch at point p c , the solution to (3.10) is given by: In the two-band case, m, n = 1, 2, then σ s mn is governed by, where σ s = ( σ s 11 σ s 12 σ s 21 σ s 22 ) T , then B s = B, C s = C and D s = D are the same as the asymmetric case, while A s is given by and ψ mn , φ mn are given by (3.2f). Noted that A s = (A s ) † .
We call Away from p c , E 1 and E 2 are well-separated, and the i ε terms for the transition coefficients σ 12 and σ 21 in Liouville-A/S lead to high oscillations, thus as ε → 0, the system formally goes to its weak limit: the classical one (2.21).
For ε > 0, around p c , however, both σ 12 and σ 21 become significant, and the band-toband transition is no longer negligible.
In fact, based on the distance to the crossing point p c , one could obtain some asymptotic properties of the transition coefficients σ mn (m = n).
Assume that the initial data for the transition coefficients are all zero, and that U x does not change sign in time for all x, take −U x > 0 for example, then: Case 3. If p ≫ C 0 √ ε, σ 12 and σ 21 are highly oscillatory with mean 0.
We leave the justification for a simpler model problem to Appendix B.
The assumption σ mn (t = 0) = 0 (m = n) is a reasonable assumption. In fact, given arbitrary initial condition, one could check that away from p c , the weak limit of σ mn is always zero, as ε → 0, as can be seen in Appendix B. So numerically we treat the initial data for both σ 12 and σ 21 zero, given the initial velocity ∂ x S 0 (x) away from the crossing point, i.e.
This assumption is intuitive and empirical, but it does give us some convenience in solving the Liouville-A/S numerically. In fact, the numerical examples provided in Section 5 indeed show that the band-to-band transition is captured very well with initial data (4.1).

A domain decomposition method
Clearly, one has several observations in hand: • the classical Liouville is an approximation (in weak limit) to Liouville-A away from the crossing point;  In the interested regime ε ∼ δ 2 , o( √ ε) mesh is enough.

Numerical examples
In this section we solve the Liouville-A and Liouville-S numerically. Firstly in Section 5.1, we present a numerical method for the equation with a linear external potential U (x).
In this special case, we provide an efficient solver without using the domain decomposition method. In Section 5.2, the domain decomposition method is applied for general U (x).
For both examples, we use the Mathieu model, i.e. the periodic potential is V (z) = cos z. very close to each other around p = 0, ±0.5.
We will focus on the 4th and 5th bands 1 . Denote Ψ 1 and Ψ 2 as the Bloch functions corresponding to the 4th and 5th bands respectively.
For comparison, we will compare the numerical results to the ones given by the original Schrödinger equation, computed through the methods given in [14] with mesh size and time step much smaller than ε.

Linear U (x)
We deal with the linear external potential in this section: Then the Schrödinger equation is: with the initial data given as a wave packet along the 4th band: Correspondingly, the Liouville-A becomes: with R given by: One encounters two computational challenges here. Firstly, one needs to numerically resolve the rapid oscillation. We will present an efficient way to overcome this difficulty by following the characteristics in Section 5.1.1. Secondly, the initial data contains a delta function.
Usually one uses a Gaussian function with small variance to approximate it, and the error is related to this variance. As will be discussed in Section 5.1.2, in some special cases, this can be avoided by using the singularity decomposition idea of [15].
Remark 5.1. Here we only discuss Liouville-A. Liouville-S can be computed similarly.

A Fourier transform based integration method
Let the Fourier transform of f (t, x, p) with respect to x bê Taking this transform on Liouville-A (5.3a), one gets: In this special case, β is a constant so the characteristic line can be obtained analytically.
and consequently one can avoid using the mesh thus removes the difficulty due to the high oscillation introduced by the term iD ε , as will be clear in the following. We take the first time step t ∈ [0, ∆t] for example. Along the characteristic line p(t) = p 0 + βt, one evaluatê σ andR at (t, p(t), η) and has: dσ dt =Rσ. (5.6) Solution to this ODE system satisfieŝ where: For t small. an approximation to (5.7) iŝ Plug (5.8b) into (5.8a), to evaluateσ 11 (∆t) andσ 12 (∆t), one needs to compute: However, to compute F 0 , F 1 and F 2 is not easy since their integrands are highly oscillatory.
But if one chooses |β|∆t = ∆p, then at each time step one follows exactly the characteristics, so p(t) always lie on the grid points, thus F 0 , F 1 and F 2 are time independent, and one only needs to compute them once (with a highly resolved calculation). Similar analysis can be carried out forσ 21 andσ 22 .
We prove the stability of this method by using it on a simpler model problem, and it will be justified in Appendix C.

A singularity decomposition idea
To handle the delta function in the initial condition (5.3b), one usually approximates it with a Gaussian function with small variance, and numerical error was determined by the width of the Gaussian. As stated before, in some special cases, this error could be avoided, and the example we are discussing here is for when ∂ x S 0 (x) ≡ p * = const, for which, we apply the singularity decomposition method introduced in [15] to reduce the error. Write the ansatz ofσ mn (t, x, p) as: in which: • θ(t, p) = p − (p * + βt), which solves the Liouville equation • ω satisfies the same equation as σ: These can be proved by simple derivations. Formally, one has The equalities above should be understood in the distributional sense. The decomposition (5.10) enables one to solve for ω and θ separately with good (bounded) initial data |a 0 (x)| 2 and ∂ x S 0 (x) respectively. The equation ω satisfies is the same as the one for σ, thus the numerical method introduced in Section 5.1.1 can be used. In the final output, one needs to get back to σ using (5.10), so a discrete delta approximation is only needed at the output time, not during time evolution.

Numerical experiments
We show the numerical results of the Liouville-A/S with the following data We compute the density, the cumulative density function (c.d.f.) and mass in the 1st band, defined respectively by: where P n is the projection onto the nth band: The two integrals in (5.13) are calculated by the midpoint quadrature rule numerically.  reflecting the 4th-to-5th band transition. The experiment also shows that smaller ε gives smaller transition rate. Note that some small oscillations occur around the crossing region.
They are related to the interference phenomena, and are usually called the Stueckelberg oscillation [5,27,29].
Define L 1 error in the cumulative distribution function (c.d.f.) [11,17]: where ρ ε S and ρ ε L denote the density calculated by the Schrödinger equation and the Liouville system respectively. Numerically we compute (5.14) using the midpoint quadrature rule.

A domain decomposition computation
This section shows examples with varying U x . For this general case, the p-characteristic is no longer a straight line, and the fast solver in the previous section is no longer valid.
To numerically solve Liouville-A, we use the domain decomposition method. The classical finite volume method is used for the convection terms.
We compute the Liouville-A system with both a pure and a mixed state initial data with:

A pure state initial data
In this example, we use the same pure state initial data as in the previous example (5.2), (5.12). Correspondingly, the initial data for the Liouville-A system is given by (5.3b).
Numerically a Gaussian function centered at p * with variance of ε/16 is used to approximate the δ−function.

Remark 5.2.
As U x varies with x, the wave packet becomes decoherent. This weakens the interference phenomenon [5,27,29]. As one can see in Figure 5.6, the Stueckelberg oscillations around the crossing region is much weaker than those in the previous example.

A mixed state initial data.
This example is for the case when the initial data is a mixed state: Correspondingly, the initial data for the semi-classical Liouville system should be: Since p * is away from the crossing point and σ 12 and σ 21 weakly converge to zero as ε → 0, numerically, we regard them as zero and use T as the initial condition, where δ G (p) is a Gaussian function centered at zero.
The density and the c.d.f. are computed for the Liouville-A and the Schrödinger, compared in Figure 5.8. Err ε as a function of ε is shown in Figure 5.9.

Conclusion
In this paper we derive semiclassical models for the linear Schrödinger equation with periodic potentials. These models take into account the band corssing, which is important to describe quantum transitions between different Bloch bands. Away from the band-crossing zones these models reduce (in the sense of weak limit) to the classical Liouville system for each Bloch band. We also couple these semiclassical models (to be used near the corssing zones) and the classical Liouville equation (used away from the crossing zones) for an efficient multiscale computation. Our numerical experiments show that these semiclassical models provide correct quantum transitions near the crossing zones when compared with the direct simulation of the original Schrödinger equation. The semi-classical Liouville system (Liouville-A) is and σ, A, C, D are defined in (3.2). Noted that A = A † and S = −S † . To check the hyperbolicity, we separate the real and imaginary parts.
σ = Re σ + i Im σ, with Re F and Im F denoting the real part and imaginary part of F respectively. Then one can get is a symmetric matrix, which implies the hyperbolicity of the system.
In addition, since S = −S † , then Re S = −Re S T , Im S = Im S T . Therefore the matrix Similarly, one can obtain the hyperbolicity of system Liouvill-S (3.17).
B Some basic analysis of the semi-classical Liouville systems To understand the asymptotic behavior of the solution to the Liouville-A system (3.2), as mentioned in Section 4.1, we look at a simpler model system:

B.1 Weak convergence
We consider the weak limit of the solution of (B.1) in this subsection. To do this, we introduce the inner product ·, · as x, p) dx dp dt.
Choose an arbitrary test function h ∈ C ∞ 0 (R + × R 2 ), take the inner product on both side of (B.1) w.r.t h, one gets The derivatives in the equation of (B.2) are acted on the smooth function h, and the left side is bounded. One gets cf, h → 0 as ε → 0 for all h ∈ C ∞ 0 (R + × R 2 ).
Given that c is almost surely nonzero, and f is bounded, one gets f ⇀ 0 weakly.
Combined with the first equation in (B.2), one gets ∂ t g + b(x)∂ p g ⇀ 0 weakly.

B.2 Strong convergence: for constant b
In these two subsections, we formally prove that before getting close to the crossing region, c(p) is assumed to be bigger than a constant c 0 that is unrelated to ε. In this region, f is constantly small and controlled by O(ε). This subsection is for the case when the speed on p direction is a constant: b(x) = β. Along the p-characteristic line p(t) = p 0 + βt, one applies the Fourier transform to the x-variable, and gets: where f (t, η) = (ĝ(t, η, p(t)),f (t, η, p(t)) ) T and The two eigenvalues of R(t) are both real, and thus the system above has a bounded solution satisfying |ĝ(t, η)| 2 + |f (t, η)| 2 = |ĝ I (η)| 2 + |f I (η)| 2 .
with r ε 22 (p) = r 22 (p) + i ε c(p), while r 11 and r 22 are purely imaginary, and c(p) real and positive. r 11 , r 22 and r 12 are independent on ε.
Set up mesh as t j = j∆t, and p i = − 1 2 + (i − 1)∆p, with ∆p and ∆t being the mesh size. Denote g j i and f j i as the numerical result at (t j , p i ), then (5. g(t, p(t)) = g j i + g(t, p(t)) r 11 (p(t)) + r 12 (p i ) t 0f (t j + s, p i + s) ds (C.1b) Plug (C.1a) into (C.1b), and evaluate them at (t j+1 , p i+1 ), one obtains: Given ε ≪ ∆t, the integrands of F 1 and F 2 are highly oscillatory, and one can see that |F 1 | ∼ O(ε) and |F 2 | ∼ O(ε 2 ). Simple calculation shows that M i can be written as With purely imaginary r 11 and r ε 22 , it is easy to prove that Ω ∞ ≤ 1, and M i ∞ ≤ (1 + O(∆t)), and thus M i ∞ ≤ (1 + O(∆t)). This implies asymptotic stability of the scheme (C.3) independent of ε → 0.