Bloch dynamics with second order Berry phase correction

Abstract

We derive the semiclassical Bloch dynamics with 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.

Keywords

Bloch dynamics Berry phase second-order correction WKB analysis

1. 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,25,33] and the references therein). In recent years, many works such as [2,3,7,14,15,30,32–34] 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 validity of the physics models and provide insight for possible generalizations (see for example [4,12,21,23] 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 non-dimensionalization, the dynamics is given by $\begin{matrix} (1.1) & i ε \frac{\partial}{\partial t} ψ (t, x) = H ψ (t, x) = (- \frac{ε^{2}}{2} Δ_{x} + V (\frac{x}{ε}) + U (x)) ψ (t, x), \end{matrix}$ where $ψ : [0, \infty) \times R^{d} \to 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 $\begin{matrix} (1.2) & i ℏ \frac{\partial}{\partial t} ψ = - \frac{ℏ^{2}}{2 m} Δ ψ + V (x) ψ + U (x) ψ, \end{matrix}$ following the procedures as in, e.g., [12]. 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 / (m L^{2})$ . In this paper we only consider the distinguished limit when $ε = h ≪ 1$ . In such setting, the typical wavelength of the wave function is comparable with the size of the lattice, which are much smaller than the observation scales in space and time. The external potential $U (x)$ is slowly-varying because $U (\cdot)$ only depends on the macroscopic length scale x, not the quantum length scale $z = x / ε$ . Such scaling has been vastly adopted in mathematical analysis of Bloch electrons, and the role of the slowly varying potential has been explored, see [4,17,28].

It is well known that such semi-classical Schrödinger equations propagate oscillations of order $O (ε)$ both in space and time. With this model, the relevant physical scale translates to the case when the typical wavelength is comparable to the period of the medium, and both of which are assumed to be small on the length-scale of the 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,12] 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 [25], Markowich, Mauser and Poupaud in [21] 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 [4,11,18] nonlinear interactions were present.

Berry phase is an important object that appears during the adiabatic limit of quantum dynamics, as some slow-changing variables enter the quantum evolution as parameters, see [2,26]. As observed by Simon in [27], 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 [15,33] and related references. Panati, Spohn and Teufel later gave a rigorous derivation of such Bloch dynamics in [23,24] 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 [12].

Recently, Gao, Yang and Niu in [15,16] constructed a second order semi-classical theory for Bloch electrons under uniform electromagnetic fields, which was based on the semiclassical Lagrangian approach originated in [31]. 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 [7,14,30,32,34]. In particular, some of the applications are connected with second order Bloch dynamics with only the external electric field as perturbation, see e.g. [14,32]. 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 two-scale approximation theory [4] serves as the basis of our work. The final results that we obtain agree with the recent paper [15] 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 [13] using Weyl quantization of operator valued symbols (see [29] for a related work).

More specifically, with a two-scaled WKB ansatz for equation (1.1), we derive the phase equation with the second order corrections, and correspondingly the perturbation in Hamiltonian and in Bloch energy. The truncated WKB solution is proved to be a valid approximation to the exact solution. The phase equation with the second order corrections can still be viewed as a Hamilton-Jacobi equation when the leading order phase equation and the amplitude equation are solved in the pre-processing stage. Note that, since we study here WKB type solutions to the Schrödinger equation, the asymptotic solution we derive is valid only before caustics. If long time validity of the asymptotic solution is desired, one needs to consider instead for example the Gaussian beam methods [10,19], semiclassically scaled wave packets [5,6], the Wigner functions [20,28], or the frozen Gaussian approximation for periodic media [8,9]. 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 [15,16] 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 - t h$ order term by $f^{n}$ , and the sum of the first $n + 1$ terms by $T^{n} f$ , namely, $\begin{matrix} f^{ε} = f^{0} + ε f^{1} + \dots + ε^{n} f^{n} + O (ε^{n + 1}) = T^{n} f + O (ε^{n + 1}) . \end{matrix}$ Note that here we use the superscripts n to indicate the $n - t h$ order term in an asymptotic expansion, since the subscript will be reserved to index different eigenvalues and eigenfunctions in the sequel.

2. 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 $\begin{matrix} (2.1) & i ε \frac{\partial}{\partial t} ψ (t, x) = H ψ (t, x) = (- \frac{ε^{2}}{2} Δ_{x} + V (\frac{x}{ε}) + U (x)) ψ (t, x) . \end{matrix}$ In the absence of the external potential U, the Hamiltonian, after a change of variable $z = x / ε$ , is given by $\begin{matrix} H_{per} = - \frac{1}{2} Δ_{z} + V (z) . \end{matrix}$ 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., [25]. In particular, we have the periodic Bloch wave functions $Ψ_{n}^{0} (z, p)$ , given as the eigenfunctions of 1

¹
The superscript 0 stands for “unperturbed”, as will be clear in the sequel.

\begin{matrix} (2.2) & H^{0} (p) Ψ_{n}^{0} (z, p) : = (\frac{1}{2} {(- i \nabla_{z} + p)}^{2} + V (z)) Ψ_{n}^{0} (z, p) = E_{n}^{0} (p) Ψ_{n}^{0} (z, p) \end{matrix}

on Γ with periodic boundary conditions. Here Γ is the unit cell of lattice

L

and

p \in Γ^{*}

is the crystal momentum, where

Γ^{*}

denotes the first Brillouin zone (unit cell of the reciprocal lattice). For each fixed

p \in Γ^{*}

, the Bloch Hamiltonian

H^{0} (p)

is a self-adjoint operator with compact resolvent, the spectrum of which is given by

\begin{matrix} σ (H^{0} (p)) = {E_{n}^{0} (p) | n \in Z_{+}} \subset R, \end{matrix}

where the eigenvalues

E_{n}^{0} (p)

(counting multiplicity) are increasingly ordered

E_{1}^{0} (p) ⩽ \dots ⩽ E_{n}^{0} (p) ⩽ E_{n + 1}^{0} (p) ⩽ \dots

. It is shown by Nenciu [22] that for any

n \in Z_{+}

there exists a closed set

C_{n} \subset Γ^{*}

of measure zero such that

\begin{matrix} E_{n - 1}^{0} (p) < E_{n}^{0} (p) < E_{n + 1}^{0} (p), p \in Γ^{*} ∖ C_{n}, \end{matrix}

and moreover

E_{n}^{0} (p)

and

Ψ_{n}^{0} (\cdot, p)

are analytic on

p \in Γ^{*} ∖ 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} = \emptyset

). As a consequence,

E_{n}^{0} (p)

and

Ψ_{n}^{0} (\cdot, 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 \in L^{2} (R^{d})$ , we recall the Bloch transform, which is an isometry from $L^{2} (R^{d})$ to $L^{2} (Γ \times Γ^{*})$ $\begin{matrix} (2.3) & \tilde{f} (z, p) = \frac{| Γ |}{{(2 π)}^{d}} \sum_{X \in L} f (z + X) e^{- i p \cdot (z + X)}, \end{matrix}$ where $| Γ |$ denotes the volume of the unit cell of the lattice $L$ . The inverse transform is given by $\begin{matrix} (2.4) & f (z) = \int_{Γ^{*}} e^{i p \cdot z} \tilde{f} (z, p) d p . \end{matrix}$

2.2. 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_{z}^{2} (Γ)$ parametrized by x and p admits the static asymptotic expansion $\begin{matrix} H^{ε} (x, p) = H^{0} (p) + ε H^{1} (x, p) + ε^{2} H^{2} (x, p) + O (ε^{3}) . \end{matrix}$ 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 Hamiltonian $H^{0} (p) = \frac{1}{2} {(- i \nabla_{z} + p)}^{2} + V (z)$ , which is independent of x. For the eigenvalue problem $\begin{matrix} H^{ε} Ψ^{ε} = E_{s}^{ε} Ψ^{ε}, \end{matrix}$ we assume the asymptotic expansions $\begin{array}{l} E_{s}^{ε} (p, x) = E^{0} (p) + ε E_{s}^{1} (p, x) + ε^{2} E_{s}^{2} (p, x) + O (ε^{3}), \\ Ψ^{ε} (z, p, x) = Ψ^{0} (z, p) + ε Ψ^{1} (z, p, x) + ε^{2} Ψ^{2} (z, p, x) + O (ε^{3}) . \end{array}$ 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 $\begin{matrix} (2.5) & H^{0} (p) Ψ^{0} (\cdot, p) = E^{0} (p) Ψ^{0} (\cdot, p), \end{matrix}$ which is just the unperturbed eigenvalue problem. In particular, $E_{s}^{0}$ 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 $\begin{matrix} (2.6) & H^{0} Ψ^{1} + H^{1} Ψ^{0} = E^{0} Ψ^{1} + E_{s}^{1} Ψ^{0} . \end{matrix}$ Taking the inner product with $Ψ^{0}$ and using the leading order equation, one gets $\begin{matrix} (2.7) & E_{s}^{1} (p, x) = ⟨ Ψ^{0} (\cdot, p), H^{1} (p, x) Ψ^{0} (\cdot, p) ⟩ . \end{matrix}$ Rewrite the first order equation as $\begin{matrix} (H^{1} - E_{s}^{1}) Ψ^{0} = (E^{0} - H^{0}) Ψ^{1} . \end{matrix}$ Given $E_{s}^{1}$ 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 $\begin{matrix} (2.8) & H^{0} Ψ^{2} + H^{1} Ψ^{1} + H^{2} Ψ^{0} = E^{0} Ψ^{2} + E_{s}^{1} Ψ^{1} + E_{s}^{2} Ψ^{0} . \end{matrix}$ To get $E_{s}^{2}$ , it suffices to take the inner product with $Ψ^{0}$ $\begin{matrix} (2.9) & \begin{matrix} E_{s}^{2} (p, x) = & ⟨ Ψ^{0} (\cdot, p), H^{1} (p, x) Ψ^{1} (\cdot, p, x) ⟩ - E_{s}^{1} (p, x) ⟨ Ψ^{0} (\cdot, p), Ψ^{1} (\cdot, p, x) ⟩ \\ + ⟨ Ψ^{0} (\cdot, p), H^{2} (p, x) Ψ^{0} (\cdot, p) ⟩ \\ = & ⟨ Ψ^{0} (\cdot, p), H^{1} (p, x) Ψ^{1} (\cdot, p, x) ⟩ + ⟨ Ψ^{0} (\cdot, p), H^{2} (p, x) Ψ^{0} (\cdot, p) ⟩, \end{matrix} \end{matrix}$ where the second equality follows from $⟨ Ψ^{1}, Ψ^{0} ⟩ = 0$ . Moreover, $Ψ^{2}$ can be solved from (2.8) given $E_{s}^{2}$ . 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 slow-varying 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 $\begin{matrix} U (x) = U (x_{c} + ε z) = U (x_{c}) + ε z \cdot \nabla U (x_{c}) + ε^{2} \frac{1}{2} z \cdot \nabla^{2} U (x_{c}) z + O (ε^{3}) . \end{matrix}$ Clearly, $ε z \cdot \nabla U (x_{c})$ corresponds the first order correction to the unperturbed Hamiltonian, and $ε^{2} \frac{1}{2} z \cdot \nabla^{2} U (x_{c}) z$ accounts for the second order correction. We remark that, the perturbation in Hamiltonian in the form of $ε z \cdot \nabla U (x_{c})$ accounts for many physics phenomenon, such as the Wannier-Stark ladders, see e.g., [22].

To see how a perturbation like $z \cdot \nabla U (x_{c})$ can be related to crystal momentum p, we note that using Bloch transformation, we have $\begin{array}{l} z f (z) = \int_{Γ^{*}} z e^{i p \cdot z} \tilde{f} (z, p) d p = \int_{Γ^{*}} - i \nabla_{p} e^{i p \cdot z} \tilde{f} (z, p) d p = \int_{Γ^{*}} i e^{i p \cdot z} \nabla_{p} \tilde{f} (z, p) d p . \end{array}$ Hence, we arrived at the first order correction to the Hamiltonian $\begin{matrix} H^{1} (p, x_{c}) = i \nabla U (x_{c}) \cdot \nabla_{p}, \end{matrix}$ and according to (2.7), the first order correction to the static Bloch energy is $\begin{matrix} (2.10) & E_{s}^{1} (p, x_{c}) = \nabla U (x_{c}) \cdot i ⟨ Ψ^{0}, \nabla_{p} Ψ^{0} ⟩ . \end{matrix}$ Similarly, by (2.9), we obtain the second order correction to the static Bloch energy $\begin{matrix} (2.11) & E_{s}^{2} (p, x_{c}) = \nabla U (x_{c}) \cdot i ⟨ Ψ^{0}, \nabla_{p} Ψ^{1} ⟩ + \frac{1}{2} \nabla^{2} U (x_{c}) : ⟨ Ψ^{0}, \nabla_{p}^{2} Ψ^{0} ⟩ . \end{matrix}$ We remark that the term $i ⟨ Ψ^{0}, \nabla_{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_{s}^{2}$ $\begin{matrix} (2.12) & {\tilde{E}}_{s}^{2} = \frac{1}{2} \nabla U (x_{c}) \cdot (i ⟨ Ψ^{0}, \nabla_{p} Ψ^{1} ⟩ + c . c .) + \frac{1}{4} \nabla^{2} U (x_{c}) : (⟨ Ψ^{0}, \nabla_{p}^{2} Ψ^{0} ⟩ + c . c .), \end{matrix}$ which naturally shows up in the WKB analysis in Section 3.

2.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 $\begin{matrix} (2.13) & A^{0} (p) = i ⟨ Ψ^{0}, \nabla_{p} Ψ^{0} ⟩, \end{matrix}$ and the first order correction of the Berry connection as $\begin{matrix} (2.14) & A^{1} (t, p, x) = \frac{1}{2} (i ⟨ Ψ^{0}, \nabla_{p} Ψ^{1} ⟩ + c . c .) . \end{matrix}$ We also introduce the velocity operator, which is also used by the physicists (see [15]), $\begin{matrix} (2.15) & \hat{V} = - i [z, H^{0}] = - i (z H^{0} - H^{0} z) . \end{matrix}$ It is easy to show that the operator $\hat{V}$ in equation (2.15) is equivalent to $\nabla_{p} H^{0}$ . In fact, $\begin{array}{l} \hat{V} f (z) & = - i (z H^{0} - H^{0} z) f (z), \\ = - i (z (\frac{1}{2} {(- i \nabla_{z} + p)}^{2} + V (z)) - (\frac{1}{2} {(- i \nabla_{z} + p)}^{2} + V (z)) z) f (z), \end{array}$ and after some simplification, we get $\begin{matrix} \hat{V} f (z) = (- i \nabla_{z} + p) f (z) . \end{matrix}$ On the other hand, from (2.2) we have $\begin{matrix} \nabla_{p} H^{0} (p) = - i \nabla_{z} + p, \end{matrix}$ and thus we conclude $\hat{V} = \nabla_{p} H^{0}$ in the operator sense.

Next, by differentiating the equation (2.2) with respect to p, we have $\begin{matrix} (2.16) & \nabla_{p} H^{0} Ψ_{m}^{0} + H^{0} \nabla_{p} Ψ_{m}^{0} = \nabla_{p} E_{m}^{0} Ψ_{m}^{0} + E_{m}^{0} \nabla_{p} Ψ_{m}^{0} . \end{matrix}$ We replace $\nabla_{p} H^{0}$ with $\hat{V}$ , take the inner product with $Ψ_{n}^{0}$ ( $n \neq m$ ), and obtain $\begin{matrix} (2.17) & ⟨ Ψ_{n}^{0}, \hat{V} Ψ_{m}^{0} ⟩ + E_{n}^{0} ⟨ Ψ_{n}^{0}, \nabla_{p} Ψ_{m}^{0} ⟩ = E_{m}^{0} ⟨ Ψ_{n}^{0}, \nabla_{p} Ψ_{m}^{0} ⟩ . \end{matrix}$ Now we have, $\begin{matrix} (2.18) & ⟨ Ψ_{n}^{0}, \nabla_{p} Ψ_{m}^{0} ⟩ = \frac{⟨ Ψ_{n}^{0}, \hat{V} Ψ_{m}^{0} ⟩}{E_{m}^{0} - E_{n}^{0}}, \end{matrix}$ or we can write it in the index form $\begin{matrix} (2.19) & ⟨ Ψ_{n}^{0}, \partial_{p_{j}} Ψ_{m}^{0} ⟩ = \frac{⟨ Ψ_{n}^{0}, {\hat{V}}_{j} Ψ_{m}^{0} ⟩}{E_{m}^{0} - E_{n}^{0}}, \end{matrix}$ where ${\hat{V}}_{j}$ is the j-th component of $\hat{V}$ .

For the rest of this section, we assume that the external potential U is linear, and thus we have $- \nabla U = \vec{E}$ , where $\vec{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 $\begin{matrix} (2.20) & H^{1} = - i \vec{E} \cdot \nabla_{p} . \end{matrix}$

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, $\begin{matrix} (2.21) & H^{0} Ψ_{1}^{1} + H^{1} Ψ_{0}^{1} = E_{0} Ψ_{1}^{1} + E^{1} Ψ_{0}^{1} . \end{matrix}$ Due to the completeness of the leading order Bloch functions, and the fact that $⟨ Ψ_{1}^{0}, Ψ_{1}^{1} ⟩ = 0$ , we write $\begin{matrix} (2.22) & Ψ_{1}^{1} = \sum_{k = 2}^{\infty} C_{k} Ψ_{k}^{0}, \end{matrix}$ and the equation (2.21) becomes $\begin{matrix} (2.23) & \sum_{k = 2}^{\infty} C_{k} E_{k}^{0} Ψ_{k}^{0} - i \vec{E} \cdot \nabla_{p} Ψ_{1}^{0} = \sum_{k = 2}^{\infty} C_{k} E_{1}^{0} Ψ_{k}^{0} + E_{1}^{1} Ψ_{0}^{0} . \end{matrix}$ Take the inner product with $Ψ_{k}^{0}$ ( $k \neq 1$ ), we obtain $\begin{matrix} (2.24) & C_{k} = \frac{i \vec{E} \cdot ⟨ Ψ_{k}^{0}, \nabla_{p} Ψ_{1}^{0} ⟩}{E_{k}^{0} - E_{1}^{0}} . \end{matrix}$

With equation (2.18), we can further rewrite $\begin{matrix} (2.25) & C_{k} = \frac{i \vec{E} \cdot ⟨ Ψ_{k}^{0}, \hat{V} Ψ_{1}^{0} ⟩}{{(E_{k}^{0} - E_{1}^{0})}^{2}} = \sum_{j} \frac{i {(\vec{E})}_{j} ⟨ Ψ_{k}^{0}, {\hat{V}}_{j} Ψ_{0}^{0} ⟩}{{(E_{k}^{0} - E_{1}^{0})}^{2}} . \end{matrix}$ 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, $\begin{matrix} (2.26) & A^{1} (t, p, x) = Re (- \sum_{k = 2}^{\infty} \frac{\vec{E} \cdot ⟨ Ψ_{k}^{0}, \hat{V} Ψ_{1}^{0} ⟩ ⟨ Ψ_{1}^{0}, \hat{V} Ψ_{k}^{0} ⟩}{{(E_{k}^{0} - E_{1}^{0})}^{3}}) \end{matrix}$ The equivalent index form is given by $\begin{matrix} A_{i}^{1} (t, p, x) = \sum_{j} Re (\sum_{k = 2}^{\infty} \frac{{(\vec{E})}_{j} ⟨ Ψ_{k}^{0}, {\hat{V}}_{j} Ψ_{1}^{0} ⟩ ⟨ Ψ_{1}^{0}, {\hat{V}}_{i} Ψ_{k}^{0} ⟩}{{(E_{1}^{0} - E_{k}^{0})}^{3}}) \end{matrix}$ In fact, this result is standard from the perspective of the perturbation theory, which exactly agrees with equation (1) of [15] 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 [15] 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 [15] and thus provide a mathematical justification.

3. 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 [12], to the Schrödinger equation (2.1). $\begin{matrix} (3.1) & ψ_{w} = A^{ε} (t, x) χ^{ε} (\frac{x}{ε}, t, \nabla_{x} S^{ε}, x) exp (\frac{i}{ε} S^{ε} (t, x)), \end{matrix}$ where $χ^{ε} (z, t, p, x)$ is the modified Bloch waves with the asymptotic expansion $\begin{matrix} χ^{ε} (z, t, p, x) = χ^{0} (z, p) + ε χ^{1} (z, t, p, x) + ε^{2} χ^{2} (z, t, p, x) + O (ε^{3}) . \end{matrix}$ We will take $\begin{matrix} χ^{0} (z, p) = Ψ^{0} (z, p), \end{matrix}$ 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 time-dependent due to the modification terms. We also assume the asymptotic expansions for the phase and amplitude in the ansatz (3.1) $\begin{array}{l} S^{ε} (t, x) = S^{0} (t, x) + ε S^{1} (t, x) + ε^{2} S^{2} (t, x) + O (ε^{3}), \end{array}$ and $\begin{matrix} A^{ε} (t, x) = A^{0} (t, x) + ε A^{1} (t, x) + ε^{2} A^{2} (t, x) + O (ε^{3}) . \end{matrix}$ 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 $\begin{matrix} ψ_{w} = a^{ε} (t, \frac{x}{ε}, x) exp (i S^{0} (t, x) / ε), \end{matrix}$ where $\begin{matrix} a^{ε} (t, \frac{x}{ε}, x) = A^{ε} (t, x) χ^{ε} (\frac{x}{ε}, t, \nabla_{x} S^{ε}, x) exp (i (S^{1} + ε S^{2} + \dots)), \end{matrix}$ 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} + \dots))$ . 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 $\begin{array}{l} e^{- i S^{ε} / ε} \partial_{t} ψ_{w} = \partial_{t} A^{ε} χ^{ε} + A^{ε} \partial_{t} χ^{ε} + \frac{i}{ε} A^{ε} χ^{ε} \partial_{t} S^{ε} + A^{ε} \nabla_{p} χ^{ε} \cdot \nabla_{x} \partial_{t} S^{ε}, \\ e^{- i S^{ε} / ε} \nabla_{x} ψ_{w} = \nabla_{x} A^{ε} χ^{ε} + \frac{i}{ε} A^{ε} χ^{ε} \nabla_{x} S^{ε} + A^{ε} \nabla_{x} χ^{ε} + A^{ε} \nabla_{x}^{2} S^{ε} \nabla_{p} χ^{ε} + \frac{1}{ε} A^{ε} \nabla_{z} χ^{ε}, \end{array}$ and $\begin{array}{l} e^{- i S^{ε} / ε} Δ_{x} ψ_{w} = & Δ_{x} A^{ε} χ^{ε} + 2 \frac{i}{ε} \nabla_{x} A^{ε} \cdot \nabla_{x} S^{ε} χ^{ε} + 2 \nabla_{x} A^{ε} \cdot \nabla_{x} χ^{ε} + 2 \nabla_{x} A^{ε} \cdot \nabla_{x}^{2} S^{ε} \nabla_{p} χ^{ε} \\ + 2 \frac{1}{ε} \nabla_{x} A^{ε} \cdot \nabla_{z} χ^{ε} + \frac{i}{ε} A^{ε} χ^{ε} Δ_{x} S^{ε} - \frac{1}{ε^{2}} {| \nabla_{x} S^{ε} |}^{2} A^{ε} χ^{ε} + 2 \frac{i}{ε} A^{ε} \nabla_{x} S^{ε} \cdot \nabla_{x} χ^{ε} \\ + 2 \frac{i}{ε} A^{ε} \nabla_{x} S^{ε} \cdot \nabla_{x}^{2} S^{ε} \nabla_{p} χ^{ε} + 2 \frac{i}{ε^{2}} A^{ε} \nabla_{x} S^{ε} \cdot \nabla_{z} χ^{ε} + A^{ε} Δ_{x} χ^{ε} \\ + 2 A^{ε} \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{x} χ^{ε} + 2 \frac{1}{ε} \nabla_{x} \cdot \nabla_{z} χ^{ε} A^{ε} + \nabla_{p} χ^{ε} \cdot (\nabla_{x} \cdot \nabla_{x}^{2} S^{ε}) A^{ε} \\ + A^{ε} {(\nabla_{x}^{2} S^{ε} \nabla_{p})}^{2} χ^{ε} + 2 \frac{1}{ε} A^{ε} \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{z} χ^{ε} + \frac{1}{ε^{2}} A^{ε} Δ_{z} χ^{ε} . \end{array}$ After sorting the terms in order of ε, we obtain $\begin{array}{l} (3.2) & i ε e^{- i S^{ε} / ε} \partial_{t} ψ_{w} = T_{0} + ε T_{1}, \\ (3.3) & e^{- i S^{ε} / ε} (- \frac{ε^{2}}{2} Δ_{x} + V (\frac{x}{ε}) + U (x)) ψ_{w} = F_{0} + ε F_{1} + ε^{2} F_{2}, \end{array}$ where we have introduced the short-hand notations $\begin{array}{l} T_{0} = - A^{ε} χ^{ε} \partial_{t} S^{ε}, T_{1} = i \partial_{t} A^{ε} χ^{ε} + i A^{ε} \partial_{t} χ^{ε} + i A^{ε} \nabla_{p} χ^{ε} \cdot \nabla_{x} \partial_{t} S^{ε}, \\ \begin{matrix} F_{0} & = \frac{1}{2} {| \nabla_{x} S^{ε} |}^{2} A^{ε} χ^{ε} - i A^{ε} \nabla_{x} S^{ε} \cdot \nabla_{z} χ^{ε} - \frac{1}{2} A^{ε} Δ_{z} χ^{ε} + V (\frac{x}{ε}) A^{ε} χ^{ε} + U (x) A^{ε} χ^{ε} \\ = (H^{0} (p) (\frac{x}{ε}, \nabla_{x} S^{ε}) + U (x)) χ^{ε} A^{ε}, \end{matrix} \\ \begin{matrix} F_{1} = & - \frac{i}{2} (Δ_{x} S^{ε} + 2 (\nabla_{x} S^{ε} - i \nabla_{z}) \cdot \nabla_{x}^{2} S^{ε} \nabla_{p}) χ^{ε} A^{ε} - i (\nabla_{x} S^{ε} - i \nabla_{z}) \cdot \nabla_{x} χ^{ε} A^{ε} \\ - i \nabla_{x} A^{0} \cdot (\nabla_{x} S^{ε} - i \nabla_{z}) χ^{ε}, \end{matrix} \\ \begin{matrix} F_{2} = & - \frac{1}{2} Δ_{x} A^{ε} χ^{ε} - \nabla_{x} A^{ε} \cdot \nabla_{x} χ^{ε} - \nabla_{x} A^{ε} \cdot \nabla_{x}^{2} S^{ε} \nabla_{p} χ^{ε} - \frac{1}{2} A^{ε} Δ_{x} χ^{ε} - A^{ε} \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{x} χ^{ε} \\ - \frac{1}{2} \nabla_{p} χ^{ε} \cdot (\nabla_{x} \cdot \nabla_{x}^{2} S^{ε}) A^{ε} - \frac{1}{2} A^{ε} {(\nabla_{x}^{2} S^{ε} \nabla_{p})}^{2} χ^{ε} . \end{matrix} \end{array}$ 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 $\begin{matrix} (3.4) & - \partial_{t} S^{0} = E^{0} (\nabla_{x} S^{0}) + U (x), \end{matrix}$ where we have used the identity (2.2). Recall that, for an isolated Bloch band, $E^{0} (p)$ is analytic for all $p \in Γ^{*}$ . 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: $\begin{matrix} \dot{Q} = P, \dot{P} = - \nabla_{P} E^{0} (Q), \end{matrix}$ with initial conditions $\begin{matrix} Q (0) = x, P (0) = \nabla_{x} S^{0} (0, x) . \end{matrix}$

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, the Hamiltonian system no long has classical solutions and the WKB solutions to the Schrödinger equation (2.1) breaks down as well. But since we aim to derive corrections to the phase equation and to Bloch dynamics, it suffices to consider the WKB solutions before caustics formation.

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., [23,33]). 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

3.2.1. 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: $\begin{matrix} H^{0} (p) χ^{0} (z, p) = E^{0} (p) χ^{0} (z, p) . \end{matrix}$ difentiating this equation with respect to p, we get $\begin{matrix} (3.5) & (p - i \nabla_{z}) χ^{0} + H^{0} (p) \nabla_{p} χ^{0} = \nabla_{p} E^{0} χ^{0} + E^{0} \nabla_{p} χ^{0} . \end{matrix}$ The inner product with $χ^{0}$ gives (since $(H_{p}^{0} - E^{0} (p)) χ^{0} = 0$ ) $\begin{matrix} (3.6) & ⟨ χ^{0}, (p - i \nabla_{z}) χ^{0} ⟩ = \nabla_{p} E^{0} . \end{matrix}$ Let us difentiate (3.5) once again with respect to p, take a product with $\nabla_{x}^{2} S^{ε}$ and sum over indices (i.e., a Frobenius inner product for matrices), we arrive at $\begin{array}{l} Δ_{x} S^{ε} χ^{0} + 2 (- i \nabla_{z} + p) \cdot \nabla_{x}^{2} S^{ε} \nabla_{p} χ^{0} + H^{0} (p) \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{p} χ^{0} \\ (3.7) & = \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{p} E^{0} χ^{0} + 2 \nabla_{x}^{2} S^{ε} \nabla_{p} E^{0} \cdot \nabla_{p} χ^{0} + E^{0} \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{p} χ^{0}, \end{array}$ which implies after taking inner product with $χ^{0}$ , $\begin{matrix} (3.8) & ⟨ χ^{0}, (Δ_{x} S^{ε} + 2 (- i \nabla_{z} + p) \cdot \nabla_{x}^{2} S^{ε} \nabla_{p}) χ^{0} ⟩ = \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{p} E^{0} + 2 \nabla_{x}^{2} S^{ε} \nabla_{p} E^{0} \cdot ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩ . \end{matrix}$ We shall use these identities in our asymptotic derivation.

3.2.2. Derivation of the phase equation with first order correction

Now we collect $O (1)$ and $O (ε)$ terms from (3.2) and (3.3) and get $\begin{matrix} (3.9) & \begin{matrix} - A^{0} \partial_{t} S^{ε} χ^{0} + ε i \partial_{t} A χ^{0} - ε A^{1} \partial_{t} S^{ε} χ^{0} - ε A^{0} \partial_{t} S^{ε} χ^{1} + i ε A^{0} \nabla_{p} χ^{0} \cdot \nabla_{x} \partial_{t} S^{ε} \\ = (H^{0} (p) (\frac{x}{ε}, \nabla_{x} S^{ε}) + U (x)) χ^{0} A^{0} \\ + ε (H^{0} (p) (\frac{x}{ε}, \nabla_{x} S^{ε}) + U (x)) χ^{1} A^{0} + ε (H^{0} (p) (\frac{x}{ε}, \nabla_{x} S^{ε}) + U (x)) χ^{0} A^{1} \\ - \frac{i ε}{2} (Δ_{x} S^{ε} + 2 (\nabla_{x} S^{ε} - i \nabla_{z}) \cdot \nabla_{x}^{2} S^{ε} \nabla_{p}) χ^{0} A^{0} \\ - i ε \nabla_{x} A^{0} \cdot (\nabla_{x} S^{ε} - i \nabla_{z}) χ^{0} . \end{matrix} \end{matrix}$ We take inner product of both hand sides with $χ^{0}$ and simplify the result using identities (3.6) and (3.8), we obtain $\begin{matrix} i ε \partial_{t} A^{0} - A^{0} \partial_{t} S^{ε} - i ε A^{0} ⟨ χ^{0}, \nabla_{q} χ^{0} ⟩ \cdot \nabla_{x} U = (E^{0} + U (x)) A^{0} - i ε \nabla_{x} A^{0} \cdot \nabla_{p} E^{0} - \frac{i ε}{2} \nabla_{x} \cdot \nabla_{p} E^{0} A^{0} . \end{matrix}$ By separating the real and imaginary parts of the above equation, we get $\begin{array}{l} (3.10) & \partial_{t} (T^{1} S) + U (x) + [E^{0} + i ε ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩ \cdot \nabla_{x} U] |_{p = \nabla_{x} (T^{1} S)} = 0, \\ (3.11) & \partial_{t} A^{0} = - \nabla_{x} A^{0} \cdot \nabla_{p} E^{0} - \frac{1}{2} \nabla_{x} \cdot \nabla_{p} E^{0} A^{0} \end{array}$ where we denote by $T^{1} S = S^{0} + ε S^{1}$ . Let us consider the structure of these two equations. The leading order equation (3.11) for the amplitude function is a transport equation. The phase equation (3.10) with the first order correction is still a Hamilton-Jacobi type equation. If we define $\begin{matrix} E^{1} (p, x) = i ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩ (p) \cdot \nabla_{x} U (x), \end{matrix}$ then the correction term in the phase equation is clearly $E^{1} (\nabla_{x} (T^{1} S), 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 $\begin{matrix} \partial_{t} \nabla_{x} (T^{1} S) + \nabla_{x}^{2} (T^{1} S) \cdot \nabla_{p} (E^{0} + ε E^{1}) + \nabla_{x} U + i ε \nabla_{x}^{2} U \cdot ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩ = 0 . \end{matrix}$ This implies the bi-characteristics in canonical variables with the notation $P = \nabla_{x} (T^{1} S)$ $\begin{array}{l} (3.12) & \dot{Q} & = \nabla_{p} (E^{0} + ε E^{1}) |_{p = P, x = Q} \\ (3.13) & \dot{P} & = - \nabla_{x} U (Q) - i ε \nabla_{x}^{2} U \cdot ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩ |_{p = P, x = Q} . \end{array}$ Here, $i ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩ = A (p)$ is known as the Berry connection or Berry vector potential. This quantity is gauge-dependent, 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 $\begin{matrix} χ^{0} \to exp (i ζ (p)) χ^{0}, \end{matrix}$ is performed, the corresponding change happens for the Berry connection, $\begin{matrix} A (p) \to A (p) - \nabla_{p} ζ (p) . \end{matrix}$ By a change of variables that takes into account the position shift between classical variables and canonical variables by the Berry connection, $\begin{matrix} Q = \tilde{Q} - ε A (\tilde{P}), P = \tilde{P}, \end{matrix}$ we obtain, after dropping the tilde, $\begin{array}{l} \dot{Q} & = \nabla_{P} E^{0} (P) + ε \nabla_{Q} U \times \nabla_{P} \times A (P), \\ \dot{P} & = - \nabla_{Q} U (Q) . \end{array}$ Here, $\nabla_{p} \times 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}, \nabla_{p} χ^{0} ⟩ \cdot \nabla_{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 [12,23,33]. 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 $\begin{array}{l} A \cdot \nabla_{x} U χ^{0} A^{0} - i \nabla_{p} χ^{0} \cdot \nabla_{x} U A^{0} \\ = + i \nabla_{x} A^{0} \cdot (H^{0} (p) - E^{0}) \nabla_{p} χ^{0} + \frac{i}{2} (H^{0} (p) - E^{0}) \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{p} χ^{0} A^{0} \\ (3.14) & + (H^{0} (p) - E^{0}) χ^{1} A^{0} . \end{array}$ 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 $χ^{1} = w + v$ , where $\begin{matrix} (3.15) & v = - \frac{i}{2} \nabla_{x}^{2} S^{0} \nabla_{p} \cdot \nabla_{p} χ^{0} - i \nabla_{x} log A^{0} \cdot \nabla_{p} χ^{0} . \end{matrix}$ Unlike v, w contains terms that cannot be written explicitly in $χ^{0}$ , which is given by $\begin{matrix} (3.16) & A \cdot \nabla_{x} U χ^{0} - i \nabla_{p} χ^{0} \cdot \nabla_{x} U = (H^{0} - E^{0}) w . \end{matrix}$ Note that this has the same structure of (2.6) (recalled here for convenience) $\begin{matrix} E^{1} Ψ^{0} - H^{1} Ψ^{0} = (H^{0} - E^{0}) Ψ^{1}, \end{matrix}$ if we identify $\begin{matrix} Ψ^{0} = χ^{0}, Ψ^{1} = w, E^{1} (p, x) = A (p) \cdot \nabla_{x} U (x), and H^{1} (p, x) f = i \nabla_{x} U (x) \cdot \nabla_{p} f . \end{matrix}$ 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 $\begin{matrix} \partial_{t} v (z, t, p, x) = - \frac{i}{2} \nabla_{x}^{2} \partial_{t} S^{0} \nabla_{p} \cdot \nabla_{p} χ^{0} - i \nabla_{x} \partial_{t} log A^{0} \cdot \nabla_{p} χ^{0}, \end{matrix}$ and hence $χ^{1}$ depends on t.

3.3. The second order corrections

3.3.1. 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 following identities, we omit the details for the straightforward derivations. $\begin{array}{l} (3.17) & ⟨ χ^{0}, (- i \nabla_{z} + p) w ⟩ + ⟨ χ^{0}, H^{1} (p) \nabla_{p} χ^{0} ⟩ = \nabla_{p} E^{1} + E^{1} ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩, \\ (3.18) & \begin{matrix} ⟨ χ^{0}, (Δ_{x} S + 2 (- i \nabla_{z} + p) \cdot \nabla_{x}^{2} S^{ε} \nabla_{p}) w ⟩ + ⟨ χ^{0}, H^{1} \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{p} χ^{0} ⟩ \\ = 2 \nabla_{x}^{2} S^{ε} \nabla_{p} E^{0} \cdot ⟨ χ^{0}, \nabla_{p} w ⟩ + \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{p} E^{1} \\ + 2 \nabla_{p} E^{1} \cdot ⟨ χ^{0}, \nabla_{x}^{2} S^{ε} \nabla_{p} χ^{0} ⟩ + E^{1} ⟨ χ^{0}, \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{p} χ^{0} ⟩, \end{matrix} \end{array}$ and $\begin{matrix} (3.19) & ⟨ χ^{0}, (- i \nabla_{z} + p) \cdot \nabla_{x} w ⟩ + ⟨ χ^{0}, \nabla_{x} \cdot (H^{1} \nabla_{p}) χ^{0} ⟩ = \nabla_{x} \cdot \nabla_{p} E^{1} + \nabla_{x} E^{1} \cdot ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩ . \end{matrix}$

Besides, from the definition of v, one gets $\begin{matrix} \frac{i}{2} \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{p} χ^{0} = - v - i \nabla_{x} log A^{0} \cdot \nabla_{p} χ^{0} . \end{matrix}$ Then, by direct substitution and simplification, we obtain the following two identities, $\begin{array}{l} (3.20) & \frac{i}{2} A^{0} E^{1} ⟨ χ^{0}, \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{p} χ^{0} ⟩ = - A^{0} E^{1} ⟨ χ^{0}, v ⟩ - i E^{1} \nabla_{x} A^{0} \cdot ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩ . \\ (3.21) & - \frac{i}{2} A^{0} ⟨ χ^{0}, H^{1} \nabla_{x}^{2} S^{ε} \nabla_{p} \cdot \nabla_{p} χ^{0} ⟩ = i A^{0} ⟨ χ^{0}, \nabla_{p} v ⟩ \cdot \nabla_{x} U - \nabla_{x} A^{0} \cdot ⟨ χ^{0}, \nabla_{p}^{2} χ^{0} ⟩ \cdot \nabla_{x} U . \end{array}$

3.3.2. Derivation of the corrected phase equation

Now we collect the terms from (3.2) and (3.3) up to $O (ε^{2})$ , this gives $\begin{array}{l} i ε \partial_{t} A^{0} χ^{0} + i ε^{2} \partial_{t} A^{1} χ^{0} + i ε^{2} \partial_{t} A^{0} χ^{1} + i ε^{2} \partial_{t} v A^{0} \\ - A^{0} \partial_{t} S^{ε} χ^{0} - ε A^{1} \partial_{t} S^{ε} χ^{0} - ε A^{0} \partial_{t} S^{ε} χ^{1} - ε^{2} A^{0} \partial_{t} S^{ε} χ^{2} - ε^{2} A^{1} \partial_{t} S^{ε} χ^{1} - ε^{2} A^{2} \partial_{t} S^{ε} χ^{0} \\ + i ε A^{0} \nabla_{p} χ^{0} \cdot \nabla_{x} \partial_{t} S^{ε} + i ε^{2} A^{1} \nabla_{p} χ^{0} \cdot \nabla_{x} \partial_{t} S^{ε} + i ε^{2} A^{0} \nabla_{p} χ^{1} \cdot \nabla_{x} \partial_{t} S^{ε} \\ = A^{0} (H^{0} (p) + U) χ^{0} + ε A^{1} (H^{0} (p) + U) χ^{0} + ε A^{0} (H^{0} (p) + U) χ^{1} \\ + ε^{2} A^{2} (H^{0} (p) + U) χ^{0} + ε^{2} A^{1} (H^{0} (p) + U) χ^{1} + ε^{2} A^{0} (H^{0} (p) + U) χ^{2} \\ - i ε \nabla_{x} A^{0} \cdot (\nabla_{x} S^{ε} - i \nabla_{z}) χ^{0} - i ε^{2} \nabla_{x} A^{1} \cdot (\nabla_{x} S^{ε} - i \nabla_{z}) χ^{0} - i ε^{2} \nabla_{x} A^{0} \cdot (\nabla_{x} S^{ε} - i \nabla_{z}) χ^{1} \\ - \frac{i ε}{2} A^{0} (Δ_{x} S^{ε} + 2 (- i \nabla_{z} + \nabla_{x} S^{ε}) \cdot \nabla_{x}^{2} S^{ε} \nabla_{p}) χ^{0} \\ - \frac{i ε^{2}}{2} A^{1} (Δ_{x} S^{ε} + 2 (- i \nabla_{z} + \nabla_{x} S^{ε}) \cdot \nabla_{x}^{2} S^{ε} \nabla_{p}) χ^{0} \\ - \frac{i ε^{2}}{2} A^{0} (Δ_{x} S^{ε} + 2 (- i \nabla_{z} + \nabla_{x} S^{ε}) \cdot \nabla_{x}^{2} S^{ε} \nabla_{p}) χ^{1} \\ - i ε^{2} (- i \nabla_{z} + \nabla_{x} S^{ε}) \cdot \nabla_{x} χ^{1} A^{0} - \frac{ε^{2}}{2} Δ_{x} A^{0} χ^{0} - ε^{2} \nabla_{x} A^{0} \cdot \nabla_{x}^{2} S^{ε} \nabla_{p} χ^{0} \\ - \frac{ε^{2}}{2} A^{0} (\nabla_{x} \cdot \nabla_{x}^{2} S^{ε}) \cdot \nabla_{p} χ^{0} - \frac{ε^{2}}{2} A^{0} {(\nabla_{x}^{2} S^{ε} \nabla_{p})}^{2} χ^{0} . \end{array}$ Take the inner product with $χ^{0}$ of the above equation, the real part of the resulting identity reveals2

²
The imaginary part leads to the equation that $A_{1}$ satisfies, which does not contribute to the second order corrections to Bloch dynamics, and will hence be neglected in this paper.

\begin{matrix} (3.22) & \partial_{t} (T^{2} S) + U (x) + [E^{0} + ε E^{1} + ε^{2} E^{2}] |_{p = \nabla_{x} (T^{2} S)} = 0, \end{matrix}

where

E^{0}

and

E^{1}

are of the form as before but are evaluated at

\nabla_{x} (T^{2} S)

instead, and we denote by

T^{2} S = 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

\begin{array}{rcl} E^{2} (p, x) & = & \frac{1}{2} (i ⟨ χ^{0}, \nabla_{p} w ⟩ + c . c .) \cdot \nabla_{x} U + ⟨ \nabla_{p} χ^{0}, \nabla_{x}^{2} U \nabla_{p} χ^{0} ⟩ - \frac{Δ_{x} A^{0}}{2 A^{0}} + \frac{1}{2} {| \nabla_{x}^{2} S^{0} \nabla_{p} χ^{0} |}^{2} \\ + \nabla_{x} log A^{0} \cdot Im (⟨ χ^{0}, (p - i \nabla_{z} - \nabla_{p} E^{0}) v ⟩) - \frac{1}{2} \nabla_{x}^{2} S^{0} \nabla_{p} \cdot \nabla_{p} E^{0} Im (⟨ χ^{0}, v ⟩) \\ - \nabla_{x}^{2} S^{0} \nabla_{p} E^{0} \cdot Im (⟨ χ^{0}, \nabla_{p} v ⟩) + \frac{1}{2} Im (⟨ χ^{0}, (Δ_{x} S^{0} + 2 (p - i \nabla_{z}) \cdot \nabla_{x}^{2} S^{0} \nabla_{p}) v ⟩) \\ (3.23) & + Im (⟨ χ^{0}, (p - i \nabla_{z}) \cdot \nabla_{x} v ⟩) + Im (⟨ χ^{0}, \partial_{t} v ⟩) . \end{array}

Note that, here

\begin{matrix} \nabla_{x} v = \nabla_{x} v (z, t, p, x) |_{z = \frac{x}{ε}, p = \nabla_{x} S_{2}} . \end{matrix}

Since

χ^{0}

does not depend on x, we have

\begin{matrix} \nabla_{x} v = - \frac{i}{2} \nabla_{x}^{3} S^{0} : \nabla_{p}^{2} χ^{0} - i \nabla_{x}^{2} log A^{0} \cdot \nabla_{p} χ^{0} . \end{matrix}

And by using equation (3.4) and (3.11), we can directly compute that

\begin{array}{l} ⟨ χ^{0}, \partial_{t} v ⟩ = & - \frac{i}{2} \nabla_{x}^{2} \partial_{t} S^{0} : ⟨ χ^{0}, \nabla_{p}^{2} χ^{0} ⟩ - i \nabla_{x} \partial_{t} log A^{0} \cdot ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩, \\ = & \frac{i}{2} (\nabla_{x}^{3} S^{0} \cdot \nabla_{p} E^{0} + (\nabla_{x}^{2} S^{0} \cdot \nabla_{p}) \otimes (\nabla_{x}^{2} S^{0} \cdot \nabla) E^{0} + \nabla_{x}^{2} U) : ⟨ χ^{0}, \nabla_{p}^{2} χ^{0} ⟩ \\ + i (\nabla_{x}^{2} log A^{0} \cdot \nabla_{p} E^{0} + \nabla_{x} log A^{0} (\nabla_{x}^{2} S^{0} : \nabla_{p}^{2} E^{0})) \cdot ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩ \\ + i (\frac{1}{2} \nabla_{x}^{3} S^{0} \cdot \nabla_{p} \cdot \nabla_{p} \cdot E^{0} + \frac{1}{2} (\nabla_{x}^{2} S^{0} \cdot \nabla_{p} \cdot \nabla_{p}) (\nabla_{x}^{2} S^{0} \cdot \nabla_{p}) E^{0}) \cdot ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩ . \end{array}

Note that, the last two lines are actually real-valued, they do not enter the corrected phase equation, so we have

\begin{matrix} Im (⟨ χ^{0}, \partial_{t} v ⟩) = \frac{1}{2} (\nabla_{x}^{3} S^{0} \cdot \nabla_{p} E^{0} + (\nabla_{x}^{2} S^{0} \cdot \nabla_{p}) \otimes (\nabla_{x}^{2} S^{0} \cdot \nabla) E^{0} + \nabla_{x}^{2} U) : Re ⟨ χ^{0}, \nabla_{p}^{2} χ^{0} ⟩ . \end{matrix}

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_{s}^{2}$ . 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 $\frac{Δ_{x} A^{0}}{2 A^{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 $\begin{matrix} (3.24) & \begin{matrix} {\tilde{E}}_{s}^{2} & = \frac{1}{2} (i ⟨ χ^{0}, \nabla_{p} w ⟩ + c . c .) \cdot \nabla_{x} U + \frac{1}{2} ⟨ \nabla_{p} χ^{0}, \nabla_{x}^{2} U \nabla_{p} χ^{0} ⟩ \\ = \frac{1}{2} (i ⟨ χ^{0}, \nabla_{p} w ⟩ + c . c .) \cdot \nabla_{x} U - \frac{1}{4} (⟨ χ^{0}, \nabla_{p}^{2} χ^{0} ⟩ + c.c.) : \nabla_{x}^{2} U . \end{matrix} \end{matrix}$ 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 ${\tilde{E}}_{s}^{2}$ (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 two-scaled WKB analysis. $\begin{array}{l} ψ_{w} & = \sum_{j = 0}^{\infty} ε^{j} g_{j} (t, \frac{x}{ε}, x) exp (i \sum_{m = 1}^{\infty} ε^{m - 1} S^{m} (t, x)) exp (\frac{i}{ε} S^{0} (t, x)) \\ (3.25) & : = \sum_{j = 0}^{\infty} ε^{j} u_{j} (t, \frac{x}{ε}, x) exp (\frac{i}{ε} S^{0} (t, x)), \end{array}$ where $\begin{matrix} u_{j} = g_{j} exp (i \sum_{m = 1}^{\infty} ε^{m - 1} S^{m} (t, x)), \end{matrix}$ and $g_{j}$ is determined by asymptotically matching the terms in (3.1). For example, $\begin{matrix} (3.26) & g_{0} = A^{0} (t, x) χ^{0} (x / ε, \nabla_{x} S^{0} (t, x)), \end{matrix}$ and $\begin{matrix} (3.27) & \begin{matrix} g_{1} = & A^{0} (t, x) χ^{1} (t, x, x / ε, \nabla_{x} S^{0} (t, x)) + A^{1} (t, x) χ^{0} (x / ε, \nabla_{x} S^{0} (t, x)) \\ + A^{0} (t, x) \nabla_{x} S^{1} (t, x) \cdot \nabla_{p} χ^{0} (x / ε, \nabla_{x} S^{0} (t, x)) . \end{matrix} \end{matrix}$

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 $\begin{matrix} ψ_{w}^{J} : = \sum_{j = 0}^{J - 1} ε^{j} u_{j} (t, \frac{x}{ε}, x) exp (\frac{i}{ε} S^{0} (t, x)) . \end{matrix}$

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. $\begin{matrix} ψ_{I}^{ε} (x) = u_{I} (x, \frac{x}{ε}) e^{i ϕ_{I} (x) / ε} + ε φ_{I}^{ε} (x), \end{matrix}$ with $ϕ_{I} \in C^{\infty} (R^{d}, R)$ , $u_{I} \in S (R^{d} \times Γ; C)$ . 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. $\begin{matrix} u_{I} (x, y) = a_{I} (x) χ_{k} (y, \nabla_{x} ϕ_{I} (x)), \end{matrix}$ where $a_{I} \in 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 \in N$ , $\begin{matrix} {‖ f^{ε} ‖}_{X_{ε}^{s}} = \sum_{| α | + | β | ⩽ s} {‖ x^{α} {(ε \partial)}^{β} f^{ε} ‖}_{L^{2}}, \end{matrix}$ and define $X_{ε}^{s}$ as: $\begin{matrix} X_{ε}^{s} = {f^{ε} \in L^{2} (R^{d}); sup_{0 < ε ⩽ 1} {‖ f^{ε} ‖}_{X_{ε}^{s}}} . \end{matrix}$ 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} \in R$ and $S^{2} \in R$ when $t = 0$ , whereas the residual part of $φ_{I}^{ε}$ is $O (ε)$ in $X_{ε}^{s}$ for all $s \in N$ .

This assumption also implies that, initially each term in the asymptotic expansion (3.1) is uniformly bounded in $X_{ε}^{s}$ for all $s \in N$ .

Note that, there are always initial conditions which satisfy 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 $ψ_{W}^{2}$ is valid up to any $t < t_{0}$ , and for all $s \in N$ there exists a constant C such that $\begin{matrix} sup_{0 < τ < t} {‖ Ψ^{ε} (τ) - ψ_{w}^{2} (τ) ‖}_{X_{ε}^{s}} ⩽ C ε^{2} . \end{matrix}$

4. 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 by $\begin{array}{l} \dot{Q} & = \nabla_{p} (E^{0} + ε E^{1}) |_{p = P, x = Q}, \\ \dot{P} & = - \nabla_{x} U (x) - i ε \nabla_{x}^{2} U \cdot ⟨ χ^{0}, \nabla_{p} χ^{0} ⟩ |_{p = P, x = Q} . \end{array}$ 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 bi-characteristic 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 following way. 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 can be treated as known quantities. We take this approach for the rest of the paper, and the physical meaning of the resulting characteristic equations are discussed in §4.2.

4.1. The bi-characteristic equations in canonical variables

Following [15,16], we call $E^{2}$ the second order wave packet energy and define $\begin{matrix} E_{w} (Q, P) = E^{2} (Q, P) - {\tilde{E}}_{s}^{2} (Q, P), \end{matrix}$ 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 ${\tilde{E}}_{s}^{2} (Q, P)$ for the rest of this section. We prove that the bi-characteristic equations with ${\tilde{E}}_{s}^{2} (Q, P)$ is gauge invariant and analyze the physical interpretation of the bi-characteristic equations. On the other hand, the extra wave packet energy may be time-dependent due to its dependence on the phase and the amplitude, and in the WKB analysis, it is given by $\begin{array}{l} E_{w} = & \frac{1}{2} ⟨ \nabla_{p} χ^{0}, \nabla_{x}^{2} U \nabla_{p} χ^{0} ⟩ - \frac{Δ_{x} A^{0}}{2 A^{0}} + \frac{1}{2} {| \nabla_{x}^{2} S^{0} \nabla_{p} χ^{0} |}^{2} \\ + \nabla_{x} log A^{0} \cdot Im (⟨ χ^{0}, (p - i \nabla_{z} - \nabla_{p} E^{0}) v ⟩) - \frac{1}{2} \nabla_{x}^{2} S^{0} \nabla_{p} \cdot \nabla_{p} E^{0} Im (⟨ χ^{0}, v ⟩) \\ - \nabla_{x}^{2} S^{0} \nabla_{p} E^{0} \cdot Im (⟨ χ^{0}, \nabla_{p} v ⟩) + \frac{1}{2} Im (⟨ χ^{0}, (Δ_{x} S^{0} + 2 (p - i \nabla_{z}) \cdot \nabla_{x}^{2} S^{0} \nabla_{p}) v ⟩) \\ + Im (⟨ χ^{0}, (p - i \nabla_{z}) \cdot \nabla_{x} v ⟩) + Im (⟨ χ^{0}, \partial_{t} v ⟩) . \end{array}$

Let us write $\begin{matrix} {\tilde{E}}_{s}^{2} = A^{1} \cdot \nabla_{x} U + B : \nabla_{x}^{2} U \end{matrix}$ with the following notations $\begin{matrix} B = - \frac{1}{4} (⟨ χ^{0}, \nabla_{p}^{2} χ^{0} ⟩ + c . c .), \end{matrix}$ and $\begin{matrix} T^{1} A = A^{0} + ε A^{1} . \end{matrix}$ We recall that $T^{1} A$ 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 [15,16], ${\tilde{E}}_{s}^{2}$ loses the contribution from $B$ because $\nabla_{x}^{2} U$ vanishes.

Then, we can write the bi-characteristic equations with the second order corrections as $\begin{array}{l} \dot{Q} & = \nabla_{P} E^{0} + ε \nabla_{P} (A^{0} \cdot \nabla_{Q} U) + ε^{2} \nabla_{P} {\tilde{E}}_{s}^{2} (Q, P) + ε^{2} \nabla_{P} E_{w}, \\ \dot{P} & = - \nabla_{Q} U - ε \nabla_{Q}^{2} U \cdot A^{0} - ε^{2} \nabla_{Q} {\tilde{E}}_{s}^{2} (Q, P) - ε^{2} \nabla_{Q} E_{w} . \end{array}$ 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 $\begin{matrix} T^{2} \tilde{H} = E^{0} (P) + U (Q) + ε E^{1} (Q, P) + ε^{2} {\tilde{E}}_{s}^{2} (Q, P) + ε^{2} E_{w} (Q, P) . \end{matrix}$ 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.

4.2. 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 [15,16]). 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.

4.2.1. Uniform electric field

To compare our results with those in recent work [15,16], we consider the special case taken in these papers where the electric field is uniform, or equivalently the scaler potential is linear. In this scenario, $\nabla_{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, $\begin{matrix} Q = Q_{c} - ε A^{0} - ε^{2} A^{1}, P = P_{c} . \end{matrix}$ Dropping the subscripts c, the bi-characteristic equations reduce to $\begin{array}{l} \dot{Q} & = \nabla_{P} E^{0} - ε \dot{P} \times \nabla_{P} \times A^{'} + ε^{2} \nabla_{P} E_{w}, \\ \dot{P} & = - \nabla_{Q} U - ε^{2} \nabla_{Q} E_{w} . \end{array}$

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. We remark that, the papers [15,16] did not give an explicit expression of $E_{w}$ , neither did they provide any analysis to this term. If we prepare the initial condition in some special form, the extra wave packet energy can be further simplified. For example, if the initial condition is a plane wave multiplied by a periodic Bloch wave, $\begin{matrix} A^{0} (0, x) = A_{0}, S^{0} (0, x) = S_{0} + K_{0} x . \end{matrix}$ 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 $\begin{matrix} A^{0} (t, x) = A_{0}, S^{0} (t, x) = b_{0} (t) + b_{1} (t) x, \end{matrix}$ where $\begin{matrix} b_{1} (t) = K_{0} - c_{1} t, b_{0} (t) = S_{0} - \int_{0}^{t} E^{0} (K_{0} - c_{1} s) d s - c_{0} t . \end{matrix}$ In other words, $A^{0} (t, x)$ stays as an constant, and $S^{0} (t, x)$ remains 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 bi-characteristic equations become $\begin{matrix} (4.1) & \{\begin{matrix} \dot{Q} & = \nabla_{P} E^{0} - ε \dot{P} \times \nabla_{P} \times T^{1} A, \\ \dot{P} & = - \nabla_{Q} U . \end{matrix} \end{matrix}$

This form essentially agrees with the recent results [15,16], 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. The works [15,16] make heuristic assumptions of the localization of the wave packet around the classical variables, which we do not need in our derivation. Our results provide a rigorous derivation of the second order correction to the Bloch dynamics in more general situations than that considered in the physics work [15,16].

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 ${\tilde{E}}_{s}^{2}$ in the second order correction. We remark that, although we have an explicit expression for $E_{w}$ depending only on leading order quantities, it is unclear yet whether $E_{w}$ is gauge invariant for general potential functions and initial conditions.

If we ignore the extra wave packet energy $E_{w}$ for this moment, which we admit needs further analysis, the bi-characteristic equations reduce to $\begin{array}{l} (4.2) & \dot{Q} & = \nabla_{P} E^{0} + ε \nabla_{P} (A^{0} \cdot \nabla_{Q} U) + ε^{2} \nabla_{P} {\tilde{E}}_{s}^{2}, \\ (4.3) & \dot{P} & = - \nabla_{Q} U - ε \nabla_{Q}^{2} U \cdot A^{0} - ε^{2} \nabla_{Q} {\tilde{E}}_{s}^{2} . \end{array}$

Notice that $\nabla_{Q}^{2} U = constant$ when Q is quadratic, and thus $\nabla_{Q} (B : \nabla_{Q}^{2} U) = 0$ , then the characteristic equations become $\begin{array}{l} \dot{Q} & = \nabla_{P} E^{0} + ε \nabla_{P} A^{0} \cdot \nabla_{Q} U + ε^{2} \nabla_{P} A^{1} \cdot \nabla_{Q} U + ε^{2} \nabla_{P} B : \nabla_{Q}^{2} U, \\ \dot{P} & = - \nabla_{Q} U - ε \nabla_{Q}^{2} U \cdot A^{0} - ε^{2} \nabla_{Q} (A^{1} \cdot \nabla_{Q} U) . \end{array}$

We carry out the same “positional-shift” change of variable as in the last case $\begin{matrix} (4.4) & Q_{c} = Q + ε A^{0} + ε^{2} A^{1}, P_{c} = P, \end{matrix}$ 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 variables $\begin{array}{l} (4.5) & \dot{P} = - \nabla_{Q} U (Q) - ε^{2} \nabla_{Q} A^{1} \cdot \nabla_{Q} U, \\ (4.6) & \begin{matrix} \dot{Q} = & \nabla_{P} E^{0} + ε \dot{P} \times \nabla_{P} \times T^{1} A - ε^{2} \nabla_{P} A^{0} \cdot (\nabla_{Q}^{2} U \cdot A^{0}) \\ + ε^{2} \nabla_{P} B : \nabla_{Q}^{2} U + ε^{2} {(\nabla_{Q} A^{1})}^{T} \cdot \nabla_{P} E^{0} . \end{matrix} \end{array}$ 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 ${\tilde{E}}_{s}^{2}$ when the potential function is quadratic in x.

In fact, as shown in Appendix C, all the quantities in (4.5) and (4.6) are gauge invariant. This seems to suggest that when the external electric field is linearly varying, and if we ignore the contribution from the extra wave packet energy, the positional-shift change of variable can also transform the characteristic equations in canonical variable to those in physical variables. We finally remark that this scenario has not been previously considered in the literature, and we leave the investigation of physical interpretations of the extra terms to future studies.

Footnotes

Acknowledgements

This work is partially supported by the National Science Foundation under grants DMS-1312659, DMS-1454939 and RNMS11-07444 (KI-Net). J.L. is also partially supported by the Alfred P. Sloan Foundation. Z. Zhou is supported by NSFC grant No. 11801016. We thank Wolfgang Gaim and Stefan Teufel for helpful discussions.

Preliminaries on Schrödinger equation

We follow the nondimensionalization steps as in [12] 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 $\begin{matrix} (A.1) & i ε \frac{\partial}{\partial t} ψ (t, x) = H ψ (t, x) = (- \frac{ε^{2}}{2} Δ_{x} + V (\frac{x}{ε}) + U (x)) ψ (t, x), \end{matrix}$ where $V (\cdot)$ 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 $\begin{matrix} (A.2) & i ℏ \frac{\partial}{\partial t} ψ = - \frac{ℏ^{2}}{2 m} Δ ψ + V (x) ψ - U (x) ψ, \end{matrix}$ where m is the mass and ℏ is the reduced Planck constant. As in [12], we introduce l as the lattice constant and $τ = m l^{2} / ℏ$ as the small (quantum) time scale, and L and T as the macroscopic length and time scales, then $\begin{matrix} V (x) = \frac{m l^{2}}{τ^{2}} \tilde{V} (\frac{x}{l}), U (x) = \frac{m L^{2}}{T^{2}} \tilde{U} (\frac{x}{L}) . \end{matrix}$ By defining $\begin{matrix} \tilde{x} = \frac{x}{L}, \tilde{t} = \frac{t}{T}, ε = \frac{l}{L}, h = \frac{ℏ T}{m L^{2}}, \end{matrix}$ one obtains after dropping the tildes, $\begin{matrix} (A.3) & i h \frac{\partial}{\partial t} ψ = - \frac{h^{2}}{2} Δ ψ + \frac{h^{2}}{ε^{2}} V (\frac{x}{ε}) ψ - U (x) ψ \end{matrix}$ This equation has two small dimensionless parameters, ε and h. We will only consider the distinguished limit when $h = ε$ , then we get (A.1)

The change of variable of the bi-characteristic equations in Section 4.2.2

The change of variable (4.4) for the P equation is straightforward: $\begin{matrix} (B.1) & \dot{P} = - \nabla_{Q} U (Q_{c}) - ε^{2} \nabla_{Q} A^{1} \cdot \nabla_{Q} U + O (ε^{3}) . \end{matrix}$

For the Q equation, notice that $\begin{matrix} \dot{Q_{c}} = \dot{Q} + ε \dot{A^{0}} + ε^{2} \dot{A^{1}}, \end{matrix}$ and $\dot{P} = - \nabla_{Q} U (Q_{c}) + O (ε^{2})$ , we have $\begin{array}{l} (B.2) & \dot{A^{0}} = - {(\nabla_{P} A^{0})}^{T} \cdot \nabla_{Q} U (Q_{c}) + O (ε^{2}), \\ (B.3) & \dot{A^{1}} = - {(\nabla_{P} A^{1})}^{T} \cdot \nabla_{Q} U (Q_{c}) + {(\nabla_{Q} A^{1})}^{T} \cdot \dot{Q} + O (ε^{2}) . \end{array}$ Thus we obtain $\begin{array}{l} \dot{Q_{c}} = & \dot{Q} + ε \dot{A^{0}} + ε^{2} \dot{A^{1}} \\ = & \nabla_{P} E^{0} + ε \nabla_{P} A^{0} \cdot \nabla_{Q} U (Q) + ε^{2} \nabla_{P} A^{1} \cdot \nabla_{Q} U (Q) + ε^{2} \nabla_{P} B : \nabla_{Q}^{2} U \\ - ε {(\nabla_{P} A^{0})}^{T} \cdot \nabla_{Q} U (Q_{c}) - ε^{2} {(\nabla_{P} A^{1})}^{T} \cdot \nabla_{Q} U (Q_{c}) + ε^{2} {(\nabla_{Q} A^{1})}^{T} \cdot \dot{Q} + O (ε^{3}) . \end{array}$

Observe that $\begin{matrix} \nabla_{Q} U (Q) = \nabla_{Q} U (Q_{c}) - ε \nabla_{Q}^{2} U (Q_{c}) \cdot A^{0} + O (ε^{2}), A^{1} (P, Q) = A^{1} (P, Q_{c}) + O (ε), \end{matrix}$ thus we get $\begin{array}{l} \dot{Q_{c}} = & \nabla_{P} E^{0} + ε \nabla_{P} A^{0} \cdot \nabla_{Q} U (Q_{c}) - ε^{2} \nabla_{P} A^{0} \cdot (\nabla_{Q}^{2} U (Q_{c}) \cdot A^{0}) + ε^{2} \nabla_{P} A^{1} \cdot \nabla_{Q} U (Q_{c}) \\ + ε^{2} \nabla_{P} B : \nabla_{Q}^{2} U - ε {(\nabla_{P} A^{0})}^{T} \cdot \nabla_{Q} U (Q_{c}) - ε^{2} {(\nabla_{P} A^{1})}^{T} \cdot \nabla_{Q} U (Q_{c}) + ε^{2} {(\nabla_{Q} A^{1})}^{T} \cdot \dot{Q_{c}}, \end{array}$ Finally, with $\dot{Q_{c}} = \nabla_{P} E^{0} + O (ε)$ and some simplification, we get $\begin{array}{rcl} \dot{Q_{c}} & = & \nabla_{P} E^{0} + ε \dot{P} \times \nabla_{P} \times T^{1} A - ε^{2} \nabla_{P} A^{0} \cdot (\nabla_{Q}^{2} U \cdot A^{0}) + ε^{2} \nabla_{P} B : \nabla_{Q}^{2} U \\ (B.4) & + ε^{2} {(\nabla_{Q} A^{1})}^{T} \cdot \nabla_{P} E^{0} + O (ε^{3}) . \end{array}$

Gauge invariance of some important terms

References

Bal,

Fannjiang,

Papanicolaou and

Ryzhik, Radiative transport in a periodic structure, J. Statist. Phys. 95 (1999), 479–494. doi:10.1023/A:1004598015978.

M.V.

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

Bloch, Ultracold quantum gases in optical lattices, Nat. Phys. 1(1) (2005), 23–30. doi:10.1038/nphys138.

Carles,

P.A.

Markowich and

Sparber, Semiclassical asymptotics for weakly nonlinear Bloch waves, J. Stat. Phys. 117(1–2) (2004), 343–375. doi:10.1023/B:JOSS.0000044070.34410.17.

Carles and

Sparber, Nonlinear dynamics of semiclassical coherent states in periodic potentials, J. Phys. A (2012), 45.

Carles and

Sparber, Semiclassical wave packet dynamics in Schrödinger equations with periodic potentials, Discrete Cont. Dyn. Syst. – Series B 17 (2012).

Culcer,

A.C.

Keser,

Y.Q.

Li and

Tkachov, Transport in two-dimensional topological materials: Recent developments in experiment and theory, 2019, preprint, arXiv:1907.10058.

Delgadillo,

Lu and

Yang, Frozen Gaussian approximation for high frequency wave propagation in periodic media, 2015, preprint, arXiv:1504.08051.

Delgadillo,

Lu and

Yang, Gauge-invariant frozen Gaussian approximation method for the Schrödinger equation with periodic potentials, 2015, preprint, arXiv:1509.05552.

10.

Dimassi,

J.C.

Guillot and

Ralston, Semiclassical asymptotics in magnetic Bloch bands, J. Phys. A: Math. Gen. 35 (2002), 7597–7605. doi:10.1088/0305-4470/35/35/304.

11.

Dohnal and

Wahlers, Coupled mode equations and gap solitons in higher dimensions, J. Diff. Equ. (2020), 269.

12.

Lu and

Yang, Asymptotic analysis of quantum dynamics in crystals: The Bloch–Wigner transform, Bloch dynamics and Berry phase, Acta Mathematicae Applicatae Sinica 29(3) (2013), 465–476. doi:10.1007/s10255-011-0095-5.

13.

Gaim and

Teufel, 2015, Private communication.

14.

Gao and

Xiao, Nonreciprocal directional dichroism induced by the quantum metric dipole, 2018, preprint, arXiv:1810.02728.

15.

Gao,

Yang and

Niu, Field induced positional shift of Bloch electrons and its dynamical implications, Phys. Rev. Lett. 112 (2014), 166601. doi:10.1103/PhysRevLett.112.166601.

16.

Gao,

Yang and

Niu, Geometrical effects in orbital magnetism, Phys. Rev. B 91 (2015), 214405. doi:10.1103/PhysRevB.91.214405.

17.

Gerard,

P.A.

Markowich,

N.J.

Mauser and

Poupaud, Homogenization limits and Wigner transforms, Comm. Pure Appl. Math. 50(4) (1997), 323–380. doi:10.1002/(SICI)1097-0312(199704)50:4<323::AID-CPA4>3.0.CO;2-C.

18.

Giannoulis,

Mielke and

Sparber, Interaction of modulated pulses in the nonlinear Schrödinger equation with periodic potential, J. Diff. Equ. (2008), 245.

19.

Jin,

Wu,

Yang and

Z.Y.

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

20.

P.L.

Lions and

Paul, Sur les mesures de Wigner, Revista Matemática Iberoamericana 9(3) (1993), 553–618. doi:10.4171/RMI/143.

21.

P.A.

Markowich,

N.J.

Mauser and

Poupaud, A Wigner function approach to semi-classical limits: Electrons in a periodic potential, Journal of Mathematical Physics 35 (1994), 1066–1094. doi:10.1063/1.530629.

22.

Nenciu, Dynamics of band electrons in electric and magnetic fields: Rigorous justification of the effective Hamiltonians, Rev. Mod. Phys. 63(1) (1991), 91–127. doi:10.1103/RevModPhys.63.91.

23.

Panati,

Spohn and

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

24.

Panati,

Spohn and

Teufel, Motions of electrons in adiabatically perturbed periodic structures, analysis, modeling and simulation of multiscale problems, in: Analysis, Modeling and Simulation of Multiscale Problems, Springer, Berlin, 2006, pp. 595–617. doi:10.1007/3-540-35657-6_22.

25.

Reed and

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

26.

Shapere and

Wilczek, Geometric Phases in Physics, World Scientific, Singapore, 1989.

27.

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

28.

Sparber,

P.A.

Markowich and

N.J.

Mauser, Wigner functions versus WKB-methods in multivalued geometrical optics, Asymptot. Anal. 33(2) (2003), 153–187.

29.

H.-M.

Stiepan and

Teufel, Semiclassical approximations for Hamiltonians with operator-valued symbols, Commun. Math. Phys. 320 (2013), 821–849. doi:10.1007/s00220-012-1650-5.

30.

H.P.

Sun and

H.Z.

Lu, Quantum transport in topological semimetals under magnetic fields (II), 2018, preprint, arXiv:1812.10120.

31.

Sundaram and

Niu, Wave-packet dynamics in slowly perturbed crystals: Gradient corrections and Berry-phase effects, Phys. Rev. B 59(23) (1999), 14915–14925. doi:10.1103/PhysRevB.59.14915.

32.

Xiao,

Z.Z.

Du and

Niu, Theory of nonlinear Hall effects: Renewed semiclassics from quantum kinetics, 2019, preprint, arXiv:1907.00577.

33.

Xiao,

M.C.

Chang and

Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82(3) (2010), 1959.

34.

Z.M.

Yu,

Liu and

S.A.

Yang, Anomalous spatial shifts in interface electronic scattering, 2019, preprint, arXiv:1901.00294.

Bloch dynamics with second order Berry phase correction

Abstract

Keywords

1. Introduction

2. Preliminaries and the static perturbation

2.1. Bloch decomposition

1 The superscript 0 stands for “unperturbed”, as will be clear in the sequel.

2.3. Deriving the first order correction of the Berry connection

3. WKB asymptotic analysis

3.1. The ansatz and the zeroth order equation

3.2. The first order corrections

3.2.1. Useful identities for Bloch waves

3.2.2. Derivation of the phase equation with first order correction

3.2.3. Derivation of H 1 and E 1

3.3. The second order corrections

3.3.1. More useful identities

3.3.2. Derivation of the corrected phase equation

2 The imaginary part leads to the equation that A 1 satisfies, which does not contribute to the second order corrections to Bloch dynamics, and will hence be neglected in this paper.

4.1. The bi-characteristic equations in canonical variables

4.2. Physical interpretation of the bi-characteristic equations

4.2.1. Uniform electric field

4.2.2. Linearly varying electric field

Footnotes

Acknowledgements

Preliminaries on Schrödinger equation

The change of variable of the bi-characteristic equations in Section 4.2.2

Gauge invariance of some important terms

References

¹
The superscript 0 stands for “unperturbed”, as will be clear in the sequel.

3.2.3. Derivation of $H^{1}$ and $E^{1}$

²
The imaginary part leads to the equation that $A_{1}$ satisfies, which does not contribute to the second order corrections to Bloch dynamics, and will hence be neglected in this paper.