{'fraction_non_alphanumeric': 0.04120638260564615, 'fraction_numerical': 0.014228100498484508, 'mean_word_length': 4.282784173613868, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bloch dynamics with second order Berry phase correction

We derive the semiclassical Bloch dynamics with the second-order Berry phase correction in the presence of the slow-varying scalar potential as perturbation. Our mathematical derivation is based on a two-scale WKB asymptotic analysis. For a uniform external electric field, the bi-characteristics system after a positional shift introduced by Berry connections agrees with the recent result in previous works. Moreover, for the case with a linear external electric field, we show that the extra terms arising in the bi-characteristics system after the positional shift are also gauge independent.


INTRODUCTION
The understanding of dynamics of Bloch electrons and their response to external electromagnetic fields plays an important role in solid state physics (see for example [1,2,21,29] and the references therein). In recent years, many works such as [2,3,5,11,12,26,[28][29][30] have explored the significant role of the Berry phase in Bloch dynamics and vast related fields. There have been series of important mathematical works in the direction as well, which are devoted to rigorously justify the the validity of the physics models and provide insight for possible generalizations (see for example [4,9,17,19] and the references therein). Under the single-particle approximation, the dynamics of an electron is treated as an independent particle on the effective periodic potential generated by ions and other electrons (as a mean-field) in the crystal. After nondimensionalization, the dynamics is given by where ψ : [0, ∞) × R d → C is a single particle wave function, ε is the semi-classical parameter, V (z) is the lattice potential which is periodic with respect to L, and U (x) is the slowly-varying scalar potential.
The dimensionless equation (1.1) is obtained from the original Schrödinger equation in physical units following the procedures as in, e.g., [9]. Here, m is the mass and ħ is the reduced Planck constant. We introduce l as the lattice constant, and L and T as the macroscopic length and time scales, respectively. Following the calculations in Appendix A), we obtain two dimensionless parameters ε = l/L, h = ħT /(mL 2 ). In this paper we only considered physical domain. This consequently leads us to a problem involving two-scales where from now on we shall denote by 0 < ε ≪ 1 the small dimensionless parameter describing the microscopic/macroscopic scale ratio.
We remark that, equation (1.1) can also be derived from the Schrödinger equation in physical units by introducing certain rescaling, which we shall omit in this paper. The readers may refer to [4,9] for such calculations.
The electronic dynamics in crystals have been studied for many years in the semi-classical regime, where the Liouville equations replace the role of the Schrödinger equation in the limit when the rescaled Planck constant tends to zero. With the help of the Bloch-Floquet theory [21], Markowich, Mauser and Poupaud in [17] derived the semi-classical Liouville equation for describing the propagation of the phase-space density for an energy band, which controls the macroscopic dynamic behavior of the electrons. Later these results were generalized to the case when a weak random potential in [1] and in [14] nonlinear interactions were present.
Berry phase is an important object that appears during the adiabatic limit of quantum dynamics, as some slowchanging variables enter the quantum evolution as parameters, see [2,22]. As observed by Simon in [23], the adiabatic Berry phase has an elegant mathematical interpretation as the holonomy of a certain connection, the Berry connection, in the appropriate fiber bundle. This setup gives rise to the Berry curvature, which is gauge invariant and can be considered as a physical observable. It has been used in the Bloch dynamics to explain various important phenomena in crystals, see for example [12,29] and related references. Panati, Spohn and Teufel later gave a rigorous derivation of such Bloch dynamics in [19,20] by writing down the effective Hamiltonian with the help of the Weyl quantization. A simple derivation of the Bloch dynamics with Berry phase correction based on WKB asymptotics was also given by E, Lu and Yang in [9].
Recently, Gao, Yang and Niu in [12,13] constructed a second order semi-classical theory for Bloch electrons under uniform electromagnetic fields, which was based on the semiclassical Lagrangian approach originated in [27].
The second order correction terms in Bloch dynamics are obtained, where the first order correction to the Berry curvature is derived. The second order semi-classical theory can be used to explain some important physical phenomenon, such as electric polarizability, magnetic susceptibility, and magnetoelectric polarizability, etc., which are investigated in the recent literature [5,11,26,28,30]. In particular, some of the applications are connected with the second order Bloch dynamics with only the external electric field as perturbation, see e.g. [11,28]. This provides the motivation to our current study: we aim to extend the mathematical derivation of Bloch dynamics to include second order corrections in the presence of an external electric field.
The main purpose of this paper is to give a rigorous derivation of the effective Bloch dynamics in crystals based on asymptotic analysis up to the second order. The justification of the corrected dynamics is not based on the semiclassical Lagrangian approach, but the WKB analysis. Since with the proper incorporation of the Bloch function, the WKB solution with high order corrections is proved in [4] to approximate the true solution of (1.1) up to a quantifiable error in an arbitrary order of ε, this approximation theory acts as the foundation of our derivation.
The final results that we obtain agree with the recent paper [12] in the situation of uniform electric field, but our approach is mathematically rigorous and is able to handle more general potentials. We also note an independent derivation of the Bloch dynamics with second order correction [10] using Weyl quantization of operator valued symbols (see [25] for a related work).
More specifically, with a two-scaled WKB ansatz for equation ( [8,15], the Wigner functions [16,24], or the frozen Gaussian approximation for periodic media [6,7]. The derivations of Bloch dynamics with Berry phase corrections using these approaches are interesting future directions. We further derive the characteristic equations of the phase equation with the second order corrections. In two specific cases: when the electric field is uniform or linear in space, we perform a physics-inspired change of variable and show that the all the terms in the resulting equations are gauge independent except those from the extra wave packet energy. In the former case, our result essentially agrees with the recent results [12,13] and thus provides a rigorous justification of these results. In the latter case of linear electric field, we find some additional gauge independent terms which appear to be physically relevant, but previously not discovered in the literature.
In this paper, we do not include the external magnetic field in the Schrödinger equation. Although the second order correction in the Bloch dynamics from the magnetic field results in many interesting applications, rigorous justification of such terms remains a challenge in mathematics, which we leave for future studies.
The rest of the paper is organized in the following way. We present a brief review of the theory of Bloch decomposition and introduce the framework of the perturbation method to the Bloch wave function in Section 2. In Section 3, we carry out a systematic two-scaled WKB analysis to the Schrödinger equation with a lattice potential and a slow-varying external potential, where the phase equation up to second order corrections has been derived and the validity of the WKB ansatz has been justified. At last, we show in Section 4 the derivation of the characteristic equations of the phase equation with second order corrections, and under certain physical assumptions, the characteristic equation reduces to the bi-characteristic equations with corrected Berry curvature.
Throughout this paper, we assume the following the convention in notation. If an ǫ dependent function f ε admits an asymptotic expansion, we denote the n−th order term by f n , and the sum of the first n +1 terms by f (n) , namely, Also we use notations as A (n) to stress that it is a n−th order tensor.

PRELIMINARIES AND THE STATIC PERTURBATION
2.1. Bloch decomposition. Recall the Schrödinger equation with a periodic lattice potential and a slow-varying scalar potential In the absence of the external potential U , the Hamiltonian, after a change of variable z = x/ε, is given by It is translational invariant with respect to the lattice L. As a result, the spectrum of the Hamiltonian can be understood by the Bloch-Floquet theory, see e.g., [21]. In particular, we have the periodic Bloch wave functions Ψ 0 n (z, p), given as the eigenfunctions of 1 on Γ with periodic boundary conditions. Here Γ is the unit cell of lattice L and p ∈ Γ * is the crystal momentum, where Γ * denotes the first Brillouin zone (unit cell of the reciprocal lattice). For each fixed p ∈ Γ * , the Bloch Hamiltonian H 0 (p) is a self-adjoint operator with compact resolvent, the spectrum of which is given by where the eigenvalues E 0 n (p) (counting multiplicity) are increasingly ordered E 0 1 (p) · · · E 0 n (p) E 0 n+1 (p) · · · . It is shown by Nenciu [18] that for any n ∈ Z + there exists a closed set C n ⊂ Γ * of measure zero such that and moreover E 0 n (p) and Ψ 0 n (·, p) are analytic on p ∈ Γ * \C n . In this paper, we stick to the adiabatic regime and assume the energy band n of interest is an isolated Bloch band (i.e., C n = ). As a consequence, E 0 n (p) and Ψ 0 n (·, p) are analytic with respect to p in Γ * . Moreover, as we will focus on the particular single band, we will suppress the subscript n unless otherwise indicated.
Given f ∈ L 2 (R d ), we recall the Bloch transform, which is an isometry from L 2 (R d ) to L 2 (Γ × Γ * ) where |Γ| denotes the volume of the unit cell of the lattice L. The inverse transform is given by e i p·z f (z, p) dp.

Perturbed Bloch wave functions.
With the external potential U , our analysis needs perturbation of the Bloch waves. For this, let us recall the perturbation theory in the context of Bloch wave functions. We assume a family of Hamiltonian H ε (x, p) on L 2 z (Γ) parametrized by x and p admits the static asymptotic expansion This expansion is called static because it does not capture the dynamical information in the time propagation. We assume the leading order term to be just given by the Bloch we assume the asymptotic expansions Here, we call E ε s the static Bloch energy, in contrast to the dynamic Bloch energy which we shall define later. But, we suppress the subcript in E 0 since E 0 is well-defined by the unperturbed Hamiltonian. The leading (zeroth) order terms in ε yields which is just the unperturbed eigenvalue problem. In particular, E 0 s and Ψ 0 are independent of the parameter x. On the other hand, the higher order terms of the Hamiltonian will depend on x explicitly, since these terms (to be specified later) will be used to capture the inhomogeneous influence by the slow-varying scalar potential U (x).
Thus, the higher order expansions of E ε and Ψ ε depend on both p and x.
Collecting the terms in the next order, one gets Taking the inner product with Ψ 0 and using the leading order equation, one gets Rewrite the first order equation as Given E 1 s and H 1 , Ψ 1 is then determined up to an arbitrary constant multiple of Ψ 0 . To fix the arbitrariness, we will take 〈Ψ 1 , Ψ 0 〉 = 0, which turns out to simplify some of the calculations in our analysis.
The equation of the second order terms reads To get E 2 s , it suffices to take the inner product with Ψ 0 (2.9) where the second equality follows from 〈Ψ 1 , Ψ 0 〉 = 0. Moreover, Ψ 2 can be solved from (2.8) given E 2 s . This procedure can be continued to higher orders, which we will omit as only the corrections up to the second order will be considered in this paper.
To derive explicit formulas for the static expansion, we consider the Hamiltonian in (2.1). To treat the slowvarying potential U (x), we assume we consider x in the neighborhood of a point x c , with x = x c + εz satisfied.
Then, we do Taylor expansion of U (x) around x c and get Clearly, εz ·∇U (x c ) corresponds the the first order correction to the unperturbed Hamiltonian, and ε 2 1 2 z ·∇ 2 U (x c )z accounts for the second order correction. We remark that, the perturbation in Hamiltonian in the form of εz · ∇U (x c ) accounts for for many physics phenomenon, such as the Wannier-Stark ladders, see e.g., [18].
To see how a perturbation like z · ∇U (x c ) can be related to crystal momentum p, we note that using Bloch transformation, we have z f (z) =ˆΓ * ze i p·z f (z, p) dp =ˆΓ * −i ∇ p e i p·z f (z, p) dp =ˆΓ * i e i p·z ∇ p f (z, p) dp.
Hence, we arrived at the first order correction to the Hamiltonian and according to (2.7), the first order correction to the static Bloch energy is Similarly, by (2.9), we obtain the second order correction to the static Bloch energy We remark that the term i 〈Ψ 0 , ∇ p Ψ 1 〉 in (2.11) is related to the first order correction to the Berry curvature, which we define in the next section. For convenience, we also define the real part of E 2 which naturally shows up in the WKB analysis in Section 3.

Deriving the first order correction of the Berry connection.
As we have mentioned in Section (2.2), some important terms in the static Bloch energy are known as the Berry connections. Here, we give the definition of the Berry connection and we shall reformulate the first order Berry connection to an expression in terms of the leading order Bloch wave functions. In Section (3), we shall see that those terms are also involved in the second order Bloch dynamics.
Define the leading order Berry connection as and the first order correction of the Berry connection as We also introduce the velocity operator, which is also used by the physicists (see [12]), It is easy to show that the operatorV in equation (2.15) is equivalent to ∇ p H 0 . In fact, and after some simplification, we getV On the other hand, from (2.2) we have and thus we concludeV = ∇ p H 0 in the operator sense.
Next, by differentiating the equation (2.2) with respect to p, we have We replace ∇ p H 0 withV , take the inner product with Ψ 0 n (n = m), and obtain For the rest of this section, we assume that the external potential U is linear, and thus we have −∇U = E , where E is the constant electric field. Please note that, the assumption is taken only to reach the agreement with the results in the physics literature. Then, the first order perturbation of the Hamiltonian is simplified to Without loss of generality, we focus on the first order correction of the Berry connection associated with the first Bloch band. Collecting the terms in the first order of the Bloch equation, Due to the completeness of the leading order Bloch functions, and the fact that 〈Ψ 0 1 , Ψ 1 1 〉 = 0, we write Take the inner product with Ψ 0 k (k = 1), we obtain With equation (2.18), we can further rewrite Substitute (2.22) and (2.25) into equation (2.14), we obtain the first order correction of the Berry connection as an expression of the leading order Bloch functions and Bloch energies, The equivalent index form is given by In fact, this result is standard from the perspective of the perturbation theory, which exactly agrees with equation (1) of [12] in the absence of the external magnetic field. Whereas, investigating the role of the Berry connection with its first order correction in dynamics is more challenging.
The work [12] proposed a derivation of the effective Lagrangian for the wave packet dynamics when the external fields are constant, and its Euler-Lagrangian equations with a positional shift lead to effective semi-classical dynamics with the second order corrections. In the next section, we aim to apply a two-scale WKB asymptotic analysis to obtain the Bloch dynamics with the second order corrections. With the same positional shift, we shall see our results essentially agree with the results in [12] and thus provide a mathematical justification.

WKB ASYMPTOTIC ANALYSIS
3.1. The ansatz and the zeroth order equation. In this section, we carry out a two-scaled mono-kinetic WKB analysis to the Schrödinger equation. The starting point is the following ansatz, which is a natural extension of the ones applied in [9], to the Schrödinger equation (2.1).
where χ ε (z, t , p, x) is the modified Bloch waves with the asymptotic expansion We will take which does not depend on t or x, and we expect χ k (z, t , p, x) contains Ψ k (z, p, x) and necessary modification terms to be specified. We emphasize that, even though Ψ k (z, p, x) is time-independent, χ k (z, t , p, x) might be timedependent due to the modification terms. We also assume the asymptotic expansions for the phase and amplitude in the ansatz (3.1) We emphasize that here the phase function series {S k } and the amplitude function series {A k } are real-valued, while the functions χ ε are complex-valued, all yet to be determined. Note that, since we will focus on a particular band throughout the analysis, we have suppressed the energy band subscript k to simplify the notation. The validity and the accuracy of this ansatz will be discussed later.
Note that, we can rewrite the ansatz in the following way is the total amplitude, which is clearly complex-valued. We remark that, the ansatz (3.1) is special in the sense that, the total amplitude a ε depend on the χ n in a restrictive way. To be more specific, since A k are all real-valued, it contains no phase information at all while all the phase information is contained in the term exp i (S 1 + εS 2 + · · · ) .
This assumption is essential to guarantee that the canonical variables have a unique trajectory. This is also why this ansatz is called mono-kinetic.
We also emphasize that, as our purpose is to derive Bloch dynamics and its corrections, we do not aim to find a general approximate solution to the Schrödinger equation up to all time, but rather a specific solution which describes the behavior of a wave packet propagating under (2.1). In particular, (3.1) poses restrictions on the initial condition, which we will discuss further below.
A straightforward calculation yields After sorting the terms in order of ε, we obtain where we have introduced the short-hand notations Combining (3.2) and (3.3) and matching by order of ε, to the leading order, we get the following Hamilton-Jacobi equation for the leading term of the phase function where we have used the identity (2.2). Recall that, for an isolated Bloch band, E 0 (p) is analytic for all p ∈ Γ * . Thus, the equation of S 0 can be solved by the method of characteristics, where the characteristic flow is determined by the following Hamiltonian equations:Q with initial conditions We remark that, the characteristic lines obtained by the above Hamiltonian flow are interpreted as the rays of geometric optics. Given the initial phase S 0 (0, x), the Hamiltonian system locally defines a flow map. Caustics may appear at some finite time when the characteristics initiated at difent locations intersect. In the event of caustics, It is natural to interpret Q and P as canonical variables of the Hamilton-Jacobi equation. In Bloch dynamics, one considers a localized wave packet, in which state the expectations of the position operator and the momentum operator are called the classical variables. In the semiclassical limit, to the leading order, the classical variables agree with the canonical variables (see e.g., [19,29]). As we will see in later sections, within the Bloch dynamical equations with corrections there is position shift introduced by Berry connections.
3.2. The first order corrections.

Useful identities for Bloch waves.
To study the first order correction, we shall first derive some useful identities from the leading order Bloch eigenvalue problem: H 0 (p)χ 0 (z, p) = E 0 (p)χ 0 (z, p).
difentiating this equation with respect to p, we get The inner product with χ 0 gives (since (H 0 p − E 0 (p))χ 0 = 0) Let us difentiate (3.5) once again with respect to p, take a product with ∇ 2 x S ε and sum over indices (i.e., a Frobenius inner product for matrices), we arrive at which implies after taking inner product with χ 0 , We shall use these identities in our asymptotic derivation.
We take inner product of both hand sides with χ 0 and simplify the result using identities (3.6) and (3.8), we obtain By separating the real and imaginary parts of the above equation, we get where we denote by S (1) = S 0 +εS 1 . Let us consider the structure of these two equations. The leading order equation then the correction term in the phase equation is clearly E 1 (∇ x S (1) , x). This correction term agrees with exactly with the first order correction to static Bloch energy, and we will identity this term with the first order correction to Bloch energy in a dynamic picture in Section 3.2.3. If we further difentiate (3.10) with respect to x, we get This implies the bi-characteristics in canonical variables with the notation P = ∇ x S (1) Here, i 〈χ 0 , ∇ p χ 0 〉 = A (p) is known as the Berry connection or Berry vector potential. This quantity is gaugedependent, which means if one chooses difent phase factors for χ 0 , the resulting A (p) are actually difent. Namely, for an arbitrary smooth function ζ(p), if the so-called gauge transformation is performed, the corresponding change happens for the Berry connection, By a change of variables that takes into account the position shift between classical variables and canonical variables by the Berry connection, we obtain, after dropping the tilde,Q = ∇ P E 0 (P ) + ε∇ Q U × ∇ P × A (P ), Here, ∇ p ×A (p) is the Berry curvature, which is gauge independent. This implies, in the corrected Bloch dynamics, an anomalous velocity is introduced by the response of Bloch electrons to the external electric field. In other words, the characteristic speed has been modified due to the correction term i ε〈χ 0 , ∇ p χ 0 〉 · ∇ x U in the phase equation.
We remark that, so far our results agree with previous works on first order corrections to the Bloch dynamics, see [9,19,29]. The focus of this paper is to extend this to second order corrections to the Bloch dynamics.

3.2.3.
Derivation of H 1 and E 1 . In this part, we aim to derive the specific expressions of H 1 (p, x) and E 1 (p, x), keeping in mind that they should satisfy the equation (2.1) and the solution ansatz (3.1). Also, we will establish the relation between χ 1 with Ψ 0 and Ψ 1 .
We observe that, since equation (3.10) and equation (3.11) have been derived, there is a difent perspective to view equation (3.9). Substitute (3.10) and (3.11) into (3.9), we obtain The main idea here is to view this equation as the first order perturbation equation for the dynamic Bloch eigenvalue problem, as it has the same structure as (2.6), in particular, all the terms with time derivatives in (3.9) are canceled.
We decompose the perturbation χ 1 into Unlike v, w contains terms that cannot be written explicitly in χ 0 , which is given by Note that this has the same structure of (2.6) (recalled here for convenience) Note here H 1 is the first order correction to the unperturbed Hamiltonian. Clearly, the first order corrections to the Bloch energy in the static expansion is of the same form as the first order correction to the phase equation, and without confusion, hereafter, we use E 1 to denote both meanings. However, in the dynamic picture, the eigenfunction χ 1 picks up an extra time-dependent part v besides Ψ 1 . This will have impact on the next order correction, as will be discussed below.
As we discussed in Section 2.2, it is natural to enforce the orthogonality condition 〈χ 0 , w〉 = 〈Ψ 0 , Ψ 1 〉 = 0. Moreover, this condition will simplify the results in Section 3.3.1.
We remark that χ 1 depend on x as a parameter, as clearly seen from the definition of v in (3.15) and the equation (3.16) for w. In addition, although w does not depend on t , v depends on t through S 0 and log A 0 , and we have and hence χ 1 depends on t .
3.3. The second order corrections.

More useful identities.
For the second order correction, we need some more identities for the Bloch waves, following similar strategies as in Section 3.2.1 applied on (3.16). We conclude with the the following identities, we omit the details for the straightforward derivations.
Besides, from the definition of v, one gets Then, by direct substitution and simplification, we obtain the following two identities,

Derivation of the corrected phase equation. Now we collect the terms from (3.2) and (3.3) up to O(ε 2 ), this
gives Take the inner product with χ 0 of the above equation, the real part of the resulting identity reveals 2 where E 0 and E 1 are of the form as before but are evaluated at ∇ x S (2) instead, and we denote by S (2) = S 0 +εS 1 +ε 2 S 2 .
Here E 2 is the second order correction term to the phase equation. Using the identities we derived before, we find that E 2 takes the following form Note that, here And by using equation (3.4) and (3.11), we can directly compute that Note that, the last two lines are actually real-valued, they do not enter the corrected phase equation, so we have Note, in the leading order phase equation, the unperturbed Bloch energy E 0 naturally shows up, and the first order correction term to the phase equation happens to be the first order correction of the static Bloch energy E 1 .
However, the second order correction to the phase equation E 2 does not necessarily agree with the second order correction of the static Bloch energy E 2 s . Due to the WKB ansatz (3.1) we have used, the phase equation and the amplitude equation are no longer decoupled with second order corrections, and thus some amplitude dependent terms, such as ∆ x A 0 2A 0 , enter the phase equation with second order corrections. Besides, in the first order correction analysis, we have learned that it is necessary to incorporate the term v in χ 1 , v related terms also contribute to the second order correction to the phase function.
Recall that we have derived the second order correction of the Bloch energy (2.12) in the static expansion, repeated here for convenience As we have mentioned, this is not necessarily the same as E 2 defined in (3.23), but E 2 contains all the terms inẼ 2 s (2.12).

3.3.3.
Well-posedness of initial data and accuracy of WKB approximation. Before we turn to the characteristic equations from the phase equation (3.22), let us summarize this section by discussing the necessary assumptions such that the ansatz (3.1) is a valid approximation and the corresponding accuracy before the formation of caustics.
To this end, we reformulate (3.1) in the following form in order to apply the convergence results in [4] for twoscaled WKB analysis.
and g j is determined by asymptotically matching the terms in (3.1). For example, Taking t = 0, the ansatz above imposes well-preparedness requirement on the initial condition to guarantee the accuracy of this approximation. For simplicity, we denote the truncated J −th order WKB approximation by To make the leading order WKB approximation valid, we make the following assumptions, similar to those used in [4].
Assumption. The initial wave function ψ ε I is in the Schwartz space S (R d ), and is of the WKB-type, i.e.
x ε e iφ I (x)/ε + εϕ ε I (x), . The function ϕ ε I is to be specified later. The amplitude u I (x, y) is assumed to be concentrated on a single isolated Bloch band E k (p) corresponding to a simple eigenvalue of H 0 , i.e.
where a I ∈ S (R d ; C) is a given amplitude.
One can check that, with the assumptions above, relation (3.26) is always satisfied at t = 0, and when t = 0, A 0 and S 0 and S 1 are uniquely determined consequently. Moreover, from equation (3.14), χ 1 is also determined at initial time.
Next, recall that we focus on wave functions which oscillate at the scale of O(ε), so we define the following spaces: for s ∈ N, and define X s ε as: Now, we are ready to state the additional assumption for the next order WKB approximation.
Assumption. (Well-prepared initial data.) The initial conditions ψ ε I (x) satisfy the assumptions above, and the leading part of the perturbation ϕ ε I is of the particular form, so that there exist solutions to A 1 ∈ R and S 2 ∈ R when t = 0, whereas the residual part of ϕ ε This assumption also implies that, initially each term in the asymptotic expansion (3.1) is uniformly bounded in X s ε for all s ∈ N . Note that, there are always initial conditions which satifiy all the assumptions above, which are WKB-type initial condition concentrated on one band with no tails. Then, by [4,Theorem 4.5], we obtain the second order approximation to the exact solution.
Theorem 3.1. Define Ψ ε (t ) to be the solution to equation (2.1) with initial conditions satisfying all the assumptions above. Assume there is no caustic formed before time t 0 , the second order WKB approximation ψ 2 W is valid up to any t < t 0 , and for all s ∈ N there exists a constant C such that

CHARACTERISTIC EQUATIONS WITH SECOND ORDER CORRECTIONS
In this section, we derive the bi-characteristic equations for the phase equation with the second order corrections. Recall that, the phase equation with the first order correction (3.10) is still a Hamilton-Jacobi equation, whose characteristic equation in canonical variables Q and P are given bẏ However, for the phase equation with the second order corrections (3.22), high order derivatives of S 0 as well as derivatives of log A 0 are involved. This means, the corrected phase equation does not admit a simple bicharacteristic structure. In other words, such a characteristic system is not yet closed.
Whereas, one observes that the higher order derivatives of S 0 and log A 0 are only contained in the terms of order O(ε 2 ). Therefore, we can proceed in two ways: (1) One can solve the leading order phase equation and amplitude equation in the pre-processing stage, then in the phase equation with second order corrections, all the higher order derivatives in the O(ε 2 ) terms are thus treated as known quantities. We take this treatment for the rest of the paper, and the physical meaning of the resulting characteristic equations are discussed in §4.2.
(2) One can close the system of trajectories by introducing the characteristic equations of those higher order derivatives appeared in the O(ε 2 ) correction. Since the results by this approach are lengthy but straightforward, we shall omit the details in this work.

The bi-characteristic equations in canonical variables.
Following [12,13], we call E 2 the second order wave packet energy and define E w (Q, P ) = E 2 (Q, P ) −Ẽ 2 s (Q, P ), where E w is interpreted as the extra wave packet energy besides the static Bloch energy due to the specific profile of thel wave function. In other words, E 2 consists of two parts, the second order correction to the static Bloch energy, and extra wave packet energy which may be time-dependent due to its dependence on the phase and the amplitude. We focus on the correction of the static Bloch energy partẼ 2 s (Q, P ) for the rest of this section. In particular, we prove that the bi-characteristic equation withẼ 2 s (Q, P ) is gauge invariant and analyze the physical interpretation of the bi-characteristic equations.
Let us writeẼ x U with the following notations , and We recall that A (1) is referred to as the Berry connection with the first order correction which responds to the first derivative of the scalar potential. While B is a new quantity which is yet to be explored. It responds to the second order derivative of the scalar potential and if only a uniform electric field is considered as a perturbation to periodic Hamiltonians as in [12,13],Ẽ 2 s loses the contribution from B because ∇ 2 x U vanishes. Then, we can write the bi-characteristic equations with the second order corrections aṡ We observe that, given that the auxiliary quantities involved in E w are known; the canonical coordinates still satisfy a Hamiltonian flow under a modified Hamiltonian We conclude this session by remarking that the bi-characteristic equations in canonical variables are valid for generic potential functions, whereas the correction terms might be gauge dependent. In the next, we proceed by considering special potential functions to investigate the bi-characteristic equations in physical variables.

Physical interpretation of the bi-characteristic equations.
To convert the characteristic equations in canonical variables to those in physical variable, one can apply a physics-motivated change of variables to guarantee that the correction terms in new variables are gauge invariant and thus physically meaningful. Such change of variables are far from being fully understood in the second order theory. When the external electrical field is uniform, the first order correction of the Berry connection introduces an additional positional shift in the change of variables (see [12,13]). We choose to take such change of variables, and discuss its applications for two scenarios: when the external electric field is uniform in space and when the field is linearly varying in space.

Uniform electric field.
To compare our results with those in recent work [12,13], we follow their setup to consider the special case taken in these papers where the electric field is uniform, or equivalently the scaler potential is linear. In this scenario, ∇ x U reduces to a constant vector, and all the higher order derivatives of U vanish.
Moreover, H 1 and E 1 no longer depend on x, and as a result, A ′ = A 0 + εA 1 becomes x−independent.
Next, we carry out the change of variables that incorporates into the position shift due to the Berry connection, Dropping the subscripts c, the bi-characteristic equations reduce tȯ At this stage, the characteristic equations have successfully captured the correction term the Berry curvature.
However, the expression for the extra wave packet energy E w is still complicated and may not be gauge invariant (in fact the expression of E w is not derived in previous works [12,13]).
We now show that it is possible to avoid the extra wave packet energy if we consider well prepared initial condition. The expression of the extra wave packet energy is explicitly given by (4.1) Let us consider the case where initial condition is a plane wave multiplied by a periodic Bloch wave, If we denote the potential function U (x) = c 0 + c 1 x, then, the exact solutions to equation (3.4) and equation (3.11) are In other words, A 0 (t , x) stays as an constant, and S 0 (t , x) remains be be a linear function in x. In this case, by direct calculations, the term v reduces to 0, and the extra wave packet energy E w simplifies to 0. Hence the bicharacteristic equations become This form essentially agrees with the recent results [12,13], although in the above special case, the extra wave packet energy E w has simplified to 0. However, our derivation is valid also for more general assumptions on the amplitude function and the phase function, and accordingly we would expect possibly different expression of the extra wave packet energy. Our results show that the second order correction to the Bloch dynamics can be derived rigorously in more general situations than that considered in [12,13].

4.2.2.
Linearly varying electric field. Next, we consider the case when the potential function U is quadratic in x, which corresponds to a linearly varying electric field. We prove in the following that the characteristic equations are gauge independent if we only keep the static Bloch energyẼ 2 s in the second order correction. The bi-characteristic equation which only include theẼ 2 s correction reduces tȯ Notice that ∇ 2 Q U = constant when Q is quadratic, and thus ∇ Q (B : ∇ 2 Q U ) = 0, then the characteristic equations becomeQ We carry out the same "positional-shift" change of variable as in the last case then after length calculations (which we provide in Appendix B), dropping the subscript c and ignoring higher order terms, we obtain the characteristic equations in physical variableṡ Compared with the uniform electric field case, more terms appear in the second order Bloch dynamics due to the extra correction term of the static Bloch energy partẼ 2 s when the potential function is quadratic in x. In fact, as shown in Appendix C, all the quantities in (4.6) and (4.7) are gauge invariant. This suggests that when the external electric field is linearly varying, the positional-shift change of variable can also transform the characteristic equations in canonical variable to those in physical variables. We remark that this scenario has not been previously considered in the literature, and investigation of physical interpretations of the extra terms is interesting research directions to be pursued.
Finally let us make some comments on the term E w in this case, in general it seems challenging to prepare initial condition to make it vanish (as in the uniform electric field case); it is also difficult to establish its gauge invariance.
In fact, since no a priori knowledge of the leading order Bloch energy E 0 is assumed, the Hamilton-Jacobi equation (3.4) does not admit a close-form solution unless the potential function U (x) is linear, which makes the analysis concerning E w intractable. While such extra wave packet energy is completely ignored in previous works, we see through the rigorous derivation that it naturally appears in the WKB analysis, and thus further understanding of the term is required to clarify its physical implications. This is beyond the scope of our current work.

APPENDIX A. PRELIMINARIES ON THE SCHRÖDINGER EQUATION
We follow the nondimensionalization steps as in [9] here, as we consider the same dimensionless form (1.1).
Given a crystal lattice L, we consider a rescaled Schrödinger equation (1.1). we rewrite it below where V (·) is a lattice potential, which is periodic with respect to L and U (x) is the slow-varying scalar potential.
(A.1) is a standard model for describing the motion of electrons in a perfect crystal when an external macroscopic potential is applied. In physical units, the equation is given by where m is the mass and ħ is the reduced Planck constant. As in [9], we introduce l as the lattice constant and τ = ml 2 /ħ as the small (quantum) time scale, and L and T as the macroscopic length and time scales, then one obtains after dropping the tildes, The change of variable (4.5) for the P equation is straightforward: For the Q equation, notice thatQ Thus we obtainQ Observe that thus we geṫ Finally, withQ c = ∇ P E 0 + O(ε) and some simplification, we get (B.2)Q c = ∇ P E 0 + εṖ × ∇ P × A (1) − ε 2 ∇ P A 0 · (∇ 2 Q U · A 0 ) + ε 2 ∇ P B : ∇ 2 Q U + ε 2 (∇ Q A 1 ) T · ∇ P E 0 + O(ε 3 ).
APPENDIX C. GAUGE INVARIANCE OF SOME IMPORTANT TERMS C.1. A 1 is gauge invariant. It suffices to show if we change χ into e i f (p) χ and then A 1 is still the same.
Recall that Thus when we change χ into e i f (p) χ, we also need to change w into e i f (p) w, and A 1 becomes the following function 1 2 (i < e i f (p) χ 0 , ∇ p e i f (p) w > +c.c.) = 1 2 (i < e i f (p) χ 0 , e i f (p) ∇ p w + i f ′ (p)e i f (p) w > +c.c.).
Hence, A 1 is gauge independent.
We then compare the expression above with ∇ 2 P f (P )(∇ 2 Q U A 0 ) and ∇ P A 0 · (∇ 2 Q U ∇ P f (P )). The k-th component of which is the opposite of part (B1).