Multiband Effects in Equations of Motion of Observables Beyond the Semiclassical Approach

The equations of motion for the position and gauge invariant crystal momentum are considered for multiband wave packets of Bloch electrons. For a localized packet in a subset of bands well-separated from the rest of the band structure of the crystal, one can construct an effective electromagnetic Hamiltonian with respect to the center of the packet. We show that the equations of motion can be obtained via a projection operator procedure, which is derived from the adiabatic approximation within perturbation theory. These relations explicitly contain information from each band captured in the expansion coefficients and energy band structure of the Bloch states as well as non-Abelian features originating from interband Berry phase properties. This general and transparent Hamiltonian-based approach is applied to a wave packet spread over a single band, a set of degenerate bands, and two linear crossing bands. The generalized equations of motion hold promise for novel effects in transport currents and Hall effect phenomena.

Introduction -Semiclassical theory has been a useful tool in describing various transport phenomena in materials and composites. This powerful theory assumes that under external fields, electrons essentially behave as free particles provided the energy band structure is used as the dispersion. Recent developments have shown that this picture is not adequate and a more rigorous treatment is needed. One of the most striking examples is the correction to the usual quasiparticle group velocity determined by the energy band structure [1]. Such an anomalous velocity correction comes from the Bloch states Berry curvature, which can lead to important modifications in many phenomena. A number of new features, such as the internal anomalous and spin Hall effects, anomalous thermoelectricity, and intrinsic magnetic moments of electronic wave packets among others have been demonstrated as a result of these new developments [2][3][4][5][6][7][8][9]. This transport formalism has also been applied to light propagation in photonic materials, where the Berry phase effects have been shown experimentally [10].
The anomalous velocity from the Berry phase is important in semiclassical transport and has been considered in several studies for the time dynamics in terms of characteristic equations of motion (EOM) of Bloch electrons in an electromagnetic field [1,5,[11][12][13][14][15]. Several phenomena linked to this anomalous velocity are especially pronounced in materials with nontrivial topology, such as Weyl and Dirac materials, where an array of intrinsic Hall effects have been demonstrated [1,16,17]. The EOM description typically assumes that the wave packet spreads out over several unit cells of a lattice and the electromagnetic waves are taken in the longwavelength approximation with a much larger spread than the packet. For such a situation, one can construct an effective Hamiltonian as a perturbative series with respect to the location of the center of the wave packet in real space [11][12][13][14]. A key issue here is the number of bands comprising the wave packet in reciprocal space. Several works have considered wave packets extending over a single band [6,8,[11][12][13], which is in line with the assumption from semiclassical the-ory that a single band far away from band degeneracies and crossings in the band structure dominates the contributions to transport. A limited number of studies have examined the dynamics of wave packets that extend over multiple degenerate bands. For example, such multiband situations have been found to influence various types of topological current transport, like parity polarization currents [14], or the time dynamics of the spin degree of freedom for degenerate bands [5].
The general scenario of describing the EOM for a wave packet extending over more than one band needs a thorough investigation, however. This is especially relevant for topological materials with spin textures, where interband effects can lead to non-Abelian gauge fields arising from an SU(N) invariance of the dynamics in a subspace of N degenerate bands. In this context, several other important issues need to be studied. For example, in deriving the wave packet dynamics, one typically invokes a projection procedure of operators onto a subset of the wave packet bands. This could be onto a single band as is the case in [13,15] or a subset of degenerate bands as is the case in [5,14]. A thorough justification of such a projection is for the most part lacking, however. Another issue is that the wave packet dynamics derived from the effective Hamiltonian is a perturbative series. The notion of order coming from this perturbative series and its effect on the EOM of the position and gauge invariant crystal momentum needs further understanding. EOM derived for bands of different dispersions and presented in a straightforward manner is also unavailable. Thus more work is needed in order to understand multiband effects in the EOM.
Here we present a direct approach for deriving EOM for the position and gauge invariant crystal momentum of multiband Bloch electron wave packets in the long-wavelength limit. Within this method, we use a perturbative series of the general electromagnetic Hamiltonian with respect to the location of the wave packet [1,13]. It is shown that for slowly time varying fields and a multiband subset of bands that are degenerate or nearly degenerate at some point in reciprocal space near the Fermi level and sufficiently separated in energy from the rest of band structure, one can apply perturbation theory with adiabaticity. From this the projected operators onto this band subset follow naturally. Full expressions for the EOM in terms of multiband quantities, like Berry curvatures and connections, are given explicitly and applied to the cases of N degenerate bands and a two-band model with a linear dispersion. In addition to its transparency, this method gives the basic framework to study multiband electron wave packet dynamics, which is to be used in future quantum kinetic transport models to capture multiband effects in transport properties.
Theoretical Model -In the presence of electromagnetic fields, quantum mechanical processes can be described using where p is the momentum operator and A(x) and ϕ(x) are the vector and scalar potentials, respectively. In a solid, the fields propagate in a periodic environment set by the crystal lattice with Bloch states giving the natural basis in reciprocal space. Transport in solids, as captured in the Boltzmann equation approach [18], involves semiclassical particles represented as wave packets [14] where ψ n (q) = e iq·x u n (q) are Bloch states with Bloch functions u n (q) for wave vector q in the n th band with expansion coefficients c n with z n = a n e −iγn(q) . Thus the wave packet is characterized by a phase γ n (q) and an amplitude ρ(q−q c ), which specifies the localization of the packet about the expectation value of its center q = q c in reciprocal space. The probability that the particle is found in the n th band is a 2 n . We assume the packet is well-localized in reciprocal space and the time scale of changing its shape is larger than the time scale of the dynamics of its center. In Fig. 1, we schematically show how such a packet spreads out over several unit cells in real space due to its narrow localization in reciprocal space. For transport, one typically works with single band wave packets in the long-wavelength approximation, where the wavelengths λ of the fields Fourier components are much longer than the spread in real space ∆x, thus = ∆x λ 1.
Neglecting the wave packet spread in reciprocal space, one then obtains the well-known semiclassical EOM [18][19][20][21] where x = x c and k = k c are the expectation values of position and gauge invariant crystal momentum k = q − e A respectively. Here ∇E(k c ) is the reciprocal space gradient of the single band dispersion of the bare crystal Hamiltonian H 0 (p, x) = p 2 2m + V (x) with H 0 (p + q, x)u n (q) = E n (q)u n (q) (∇ is understood as differentiation with respect to k c ) and E and B are the external electric and magnetic fields respectively. The EOM can then be utilized within the Boltzmann equation to determine electric and heat currents in a given material [18].
It has been recognized, however, that the interplay between the time dynamics of the fields and the material band structure has to be re-examined within the semiclassical approximation [1,7,[11][12][13]. Specifically, due to the long-wavelength approximation, one can effectively use a power series of the field potentials in the Hamiltonian about the real space wave packet center x c [13]. Retaining only up to the first order in this series, one finds Here is the vector potential and ϕ c = ϕ(x c ) is the scalar potential, which are all functions of the center of the wave packet x c . The last two terms in (3) determine the first order perturbation of the linearized Hamiltonian in the long-wavelength limit via the parameter = ∆x λ . This is realized by noting that ∂Hc ∂xc = j ∂Hc ∂Ac,j ∂Ac,j ∂xc results in a 1/λ factor from the Fourier transformed vector potential (see (S1) in the Supplementary Information), while (x − x c ) determines the wave packet spread in real space ∆x. Using a Lagrangian formalism within the semiclassical single band approximation [11][12][13], the EOM are found aṡ whereΩ n = ∇ ×Ā n is the Berry curvature andĀ n = i u n (k c )|∇|u n (k c ) is the Berry connection with Bloch functions u n (k c ) for the n th band. The energy of the single band packet that couples to the magnetic field B.
One notes that H c (p−eA c , x; x c ) in the Hamiltonian from (3) leads to the standard well-known EOM in (2). The perturbation term proportional to (x − x c ) in (3) is responsible for the appearance of the Berry curvature in˙ x . Relations (4) are considered more complete as compared to (2) since they account for geometrical features associated with the Berry curvature of the single band responsible for the transport in a given material. The relations in (4) have also been obtained by a Hamiltonian approach [14] which uses a projection procedure of operators onto a wave packet subspace consisting of FIG. 2. Schematics of the D energy band region with a band width δE bw (yellow) separated from the rest of the band structure (gray) with a large gap ∆E, such that δE bw ∆E. As an example, two linear bands in the D region with a crossing point in the vicinity of the Fermi level EF are shown. a single energy band. The resulting operators from this projection lead to noncanonical commutation relations giving rise to the Berry phase contributions in˙ x . Another approach [15], based on a semiclassical expansion and diagonalization of the electromagnetic Hamiltonian, has also resulted in the same EOM (in (4)). This method also relies on projected operators leading to noncanonical commutation relations [14,15,22].
These recent developments have shown that to take into account features from the Berry curvature in the EOM, one must carefully consider the slow time dynamics of the electromagnetic fields in the Hamiltonian. Another key issue is the application of the operator projection procedure as means to capture the role of a finite subset of bands that determines the transport in the solid. Interestingly, a rigorous derivation of such a projection has not been given yet, although some justification within a semiclasscal expansion has been offered in [23]. Also, the majority of the works have focused on a single band projection with the exception of [5,24] where a finite set of identical bands of a multiband packet was considered. We note that the EOM in (4) constitute a significant step towards a more complete understanding of transport in materials, however the outstanding questions of justifying the projection procedure and explicitly taking into account multiband wave packet with different dispersions must be resolved, especially for situations when the Fermi level is in close proximity to crossing points between bands.
Projected EOM -We first show that the projection procedure of operators follows naturally from perturbation theory and an adiabatic approximation. For this purpose, we consider the following conditions: (i) there is a Hamiltonian of the form H = H 0 (α(t)) + H 1 (α(t)) where the unperturbed part H 0 may be time-dependent due to a set of time-dependent parameters, collectively denoted as α(t), that may also be present in the perturbation H 1 (α(t)) with perturbation parameter , (ii) α(t) is turned on and thereafter varies slowly in time, and (iii) the eigenfunctions of H 0 (α(t)) are ψ 0 n (α(t)) and are well separated from the other eigenstates by a large energy gap for all time. The subspace of these states is denoted as D and indexed by I D , with ψ 0 n ∈ D and n ∈ I D . By applying perturbation theory in tandem with adiabaticity under these circumstances, the time dynamics of an observable O is found to have the following form, where P D is the projection operator onto D. Here c n are the coefficients for the expansion of the wave packet ψ = n∈I D c n ψ 0 n and their coupled dynamics is essentially given from the Schrödinger equation restricted to the subspace D. We emphasize that the projection operator is a natural consequence of the assumptions for the adiabatic approximation and the large energy gap separation condition as described above. Extensive details of the derivation for the dynamics of O and the expansion coefficients c n are shown in the Supplemental Information.
We wish to apply the projected dynamics in (5) to the Hamiltonian in (3) for a multiband wave packet centered in reciprocal space about wave vectors with energies in the vicinity of the Fermi level. In this situation, the subspace D consists of energy band states in the multiband wave packet near the Fermi level, which may also include crossing points. In the context of this perturbation theory with adiabaticity, the band width of the energy band states in D is δE bw , which is separated by a large energy gap ∆E from the rest of the band structure (δE bw ∆E), as schematically illustrated in Fig. 2. As explained in [13] and the Supplementary Information, the series expansion of the electromagnetic Hamiltonian in (3) can be considered as having an unperturbed term coming from the bare crystal potential H 0 = H c (p−eA c , x; x c ) and a perturbation H 1 = 1 2 ∂Hc ∂xc · (x − x c ) + h.c., associated with the derivatives of the electromagnetic potentials A c and ϕ c with respect to x c . The potentials A c and ϕ c serve as the collective parameter α(t) in conditions (i)-(iii). In addition to the perturbation as described, the slow time scale of these potentials meets conditions (i) and (ii). Furthermore, in accordance with the wave packet description, the initial wave function is a multiband packet consisting of Bloch states as shown in (1). The band width δE bw of states near a Fermi level tuned to a crossing point is well separated from the rest of the band structure (δE bw ∆E), leading to condition (iii). In Fig. 2, we schematically show this situation of two crossing bands with linear dispersion in the vicinity of the Fermi level. Given that conditions (i)-(iii) are met, the EOM for x and k subject to the perturbative Hamiltonian in (3) are obtained using (5) in what follows for multiband packets.
Results for the EOM -Relations (5) provide a rather straight forward (although tedious) path to obtain the projected EOM for an observable provided the Hamiltonian for the system is given explicitly. The problem essentially involves working out several commutators. Here (5) is applied to obtain the EOM for x and k for the electromagnetic Hamiltonian in its linearized form (in (3)), using the eigenstates ψ 0 n (q) = e iq·x u 0 n (q − eA c ). For clarity and concreteness we take the circular gauge with A = 1 2 (B × x) and ϕ = −E · x, where the fields may be slowly varying in time. After some calculations (details are given in the Supplementary Information), we finḋ where summation over all bands in the multiband wave packet for each band index n 1 , n 2 , and n 3 in each term in˙ x and˙ k is implied. We have defined the interband Berry connectionĀ n1n2 = i u n1 (k c )|∇|u n2 (k c ) and interband Berry curvatureΩ n1n2 = ∇ ×Ā n1n2 , as well as an interband magnetizationM n1n2 The non-Abelian coupling between interband Berry connections is denoted as Q n1n2n3 =Ā n1n2 ×Ā n2n3 . We note that x c and k c are expectation values calculated with the wave function for the wave packet. The coefficients c n are found from first solving for the z n in (6) and then using (1). Also, v = n∈I D a 2 n v n depends only on the wave packet distribution across bands and their dispersions, v n = 1 ∇E n , ν γn = ∇γ n , and ∇ is understood as differentiation with respect to k c everywhere.
A comparison of the above results with the singleband EOM in (4) shows that the results for˙ x are much more complicated for the multiband packet. The first term in (6a) is simply the standard group velocity weighted over all bands. The remaining terms on the first line and all terms on the second line capture various aspects of the different band dispersions. In the singleband case, these dispersion terms reduce to the ∇E(k c ) contribution in˙ x in (4). All terms in the third line and the terms containing v n in the fourth line correspond to the anomalous velocity, which now reflects the different dispersions and interband Berry curvatures and connections. Clearly, this is more complicated when comparing with the singleband˙ k ×Ω n anomalous velocity in (4). The rest of the terms on the last line consist of energy contributions of the magnetic field coupled to the interband magnetizations, which corresponds to theM n · B contribution to˙ x in (4). Surprisingly, the only difference between the˙ k in (6) and that in (4) is the v that appears in (6) instead of the full expression for˙ x . The reason that˙ k in (6) is similar to˙ k in (4) is because k is diagonal in the zeroth order Bloch basis. Physically, the real space wave packet doesn't "see" the distribution across bands in reciprocal space in regards to the momentum and we thus retain a Lorentz-like force expression for˙ k .
We examine (6) in some specific cases. Firstly, the single band wave packet EOM reduce to˙ , which is the same as (4) except for the v that appears in˙ k here as opposed to˙ x . We also consider the case of degenerate bands, i.e. E n = E for all N bands, for which (6) reduces to All notation and conventions for summation over band indices and ∇ differentiation are the same as in (6). It is noted that the above results are of the same form as in [14], where the wave packet localized in N degenerate bands was also considered.
Comparing the results in (7a) with the general expression in (6a) shows that the only terms that do not vanish in˙ x are those coming from the coupling between the magnetic field and interband magnetizations. Also, the anomalous velocity now contains only terms coupling˙ k toΩ n1n2 andQ n1n3n2 , which are manifestations of the non-canonical commutation relations of projected position operators onto the multiband subspace. Finally, (6) is applied to the case of a pair of linear bands with dispersion E ± (k) = ±vk, where v is the isotropic Fermi velocity (schematics in Fig. 2). Such a band structure corresponds to a Hamiltonian H = vσ · k (σ are the Pauli matrices) in the vicinity of a linear two band crossing in a Weyl semimetal or graphene, for example [16,17]. where k , which suggests that this term would vanish if the distribution across the two bands is equal. Because of the linear dispersion many terms in the general case in (6a) vanish or become modified. Specifically, there is only one term coming from the coupling between the interband magnetization, while dispersion terms with v n1 − v n2 account for differences in energy band slopes. The anomalous velocity (last term in (8a)) is now associated with the interband Berry curvature and the band dispersions. We note that there is an explicit dependence on a ± , γ ± in the exponents and in the ν γ± , thus to obtain more specific results, detailed knowledge of the wave packet properties is needed, which is beyond the scope of this paper.
Conclusions -In this study, the time dynamics of the position and gauge invariant crystal momentum subject to elec-tromagnetic fields in a periodic environment is considered using a general Hamiltonian constructed for a localized multiband wave packet. One of the main contributions here is showing that when the wave packet is spread over a subset of bands well-separated from the rest of the band structure of the crystal, a projection procedure of the Hamiltonian and EOM can be applied. We demonstrate that such a procedure follows naturally from an adiabatic approximation within perturbation theory, whereas this justification is mostly lacking in previous works. This projection procedure reflects that transport is dominated by a finite number of energy bands in the vicinity of the Fermi level. Another merit of our work is that this Hamiltonian-based approach is rather transparent and accounts for different energy band dispersions that make up the wave packet. Perhaps the most noteworthy result comes from the found interband effects in the anomalous velocity. These are captured via non-Abelian Berry characteristics giving new perspectives for nontrivial topology. The application of this general theory to several cases (including two linear crossing bands) shows that the generalized EOM depend strongly on the explicit band structure of the material and quantum mechanical effects beyond the single band wave packet description can lead to novel properties.
To further understand how specific features in the EOM for a multiband wave packet affect transport however, one must utilize a multiband generalization of the Boltzmann equation [18] in order to obtain specific transport properties, such as electric, heat, and spin currents. The EOM described in this work constitute much needed ground for such a generalized transport approach that will enable finding transport currents, Hall effects and other properties influenced by interband properties like the non-Abelian Berry curvature, which could lead to measurable signatures in the laboratory. We

ADIABATIC APPROXIMATION AND PERTURBATION THEORY
Here we demonstrate in detail how the equations of motion (EOM) projected onto a subset of locally degenerate or nearly degenerate bands arise from perturbation theory with an adiabatic condition ( (5) and (6) from the main text). We begin with a general situation, which will be related to the situation in the main text.
General Model: We consider a general Hamiltonian given as H(α(t)) = H 0 (α(t)) + H 1 (α(t)), where the time dependence is characterized by the composite parameter α(t). H 0 has a degenerate spectrum with degenerate subspace D and the perturbation H 1 (α(t)), is turned on at time t = 0. The eigenstates for the full Hamiltonian are the set {ψ n } with Hψ n (α(t)) = E n (α(t))ψ n (α(t)) and the eigenstates for the unperturbed Hamiltonian are the set {ψ 0 n } with H 0 ψ 0 n (α(t)) = E 0 n (α(t))ψ 0 n (α(t)). For the following, we utilize states ψ 0 n ∈ D that diagonalize H 1 in D. We require a large energy gap condition on D, meaning that the energies of eigenstates in D are separated from all other energies by a large energy gap. We also assume that a particle starts out in the degenerate subspace D of H 0 at t = 0 and the large energy gap condition is maintained for the instantaneous energies of H 0 while D changes in time.
Perturbation: It is noted that the Hamiltonian in the derivation of the EOM of operators in the main text involves the perturbation H 1 = 1 2 ∂Hc ∂xc · (x − x c ) + h.c., where the perturbation parameter is = ∆x λ with ∆x as the wave packet spread in real space and λ as the wavelength of the Fourier components of the vector potential. It's clear that ∆x comes from the operator (x − x c ). To show that ∂Hc ∂xc is proportional to 1 λ , we first simply observe that ∂Hc ∂xc = j ∂Hc ∂Ac,j ∂Ac,j ∂xc . Now we write A c as a Fourier transform, Here q = 2π λ is the magnitude of the Fourier wave vector component with corresponding wavelength λ and x c is the location of the center of the wave packet in real space. In the long-wavelength limit, only those q with large λ contribute appreciably to the Fourier transform. Therefore,Ã(q) is non-negligible only in the vicinity of q ≈ 0. Additionally, one finds for the spatial derivative ∂Ac,j ∂x c,k = i q kÃj (q)e −iq·xc d 3 q ∝ 1 λ . Expansion Coefficients: Our goal is to study the dynamics of Bloch wave packets in the vicinity of crossing points in the band structure for the linearized Hamiltonian in the main text (Eq. (3)). If the Fermi level of our system is tuned to such a crossing point, then the Bloch wave packet consists of these degenerate states. The degenerate states gives rise to the degenerate subspace D that is separated from all other band states by a large energy gap. Therefore, the goal to describe the dynamics of such Bloch wave packets meets the requirements of the general approach outlined above.
Returning to the general situation, we start by expanding the full wave function ψ(t) for the particle in terms of the full Hamiltonian's instantaneous eigenstates ψ n (t), where the sum is over the index set I H for a complete set of instantaneous eigenstates in the Hilbert space H. One notes that each coefficient c j m,n for j > 0 contains factors of the form for some states ψ 0 k , ψ 0 l and energies E 0 k , E 0 l . Since states in D are separated from all other states by a large energy gap, such terms are negligible when either ψ 0 k / ∈ D and ψ 0 l ∈ D or vice versa, and thus (using the index set I D for the set of basis eigenstates that span D) The above result implies that ψ m is approximated as simply ψ 0 m for m ∈ I D and is entirely contained outside of D when m / ∈ I D . Next, the k th order of the expansion coefficients from (S3) are considered and can be cast aṡ It is realized that in the zeroth order term,ċ 0 n = − i c 0 n E 0 n − m∈I H c 0 m ψ 0 n |ψ 0 m , one can make use of the adiabatic approximation in the term containing ψ 0 n |ψ 0 m , givinġ (S7) From the above relation it is realized that since there is no perturbation at t = 0, then Consequently, there is no zeroth order contribution to the coefficients c n for n / ∈ I D as expected from adiabatic theory. To zeroth order, the particle stays in the degenerate space D.
We now claim that c n ≈ 0 to all orders for n / ∈ I D , which can be shown by induction. Suppose that c j n = 0 for all j < k when n / ∈ I D . Thus using (S6), for n / ∈ I D , By invoking the adiabatic approximation, the only term in the first sum that is non-negligible occurs when m = n. By invoking the large energy gap condition and (S5), the terms in the second sum are negligible. Therefore, the non-negligible terms giveċ However, there is no is perturbation at t = 0 and thus c k n = 0, n / ∈ I D . (S11) Consequently, there is no k th order contribution to the coefficients c n for n / ∈ I D . Coupled to the fact that c 0 n ≈ 0 for n / ∈ I D , induction shows that c n = 0 to all orders for n / ∈ I D . Therefore, (S3) transforms tȯ where the summation now runs over I D instead of I H . We are interested in (S12) when n ∈ I D . By using the large energy gap condition and (S5), (S12), it is found thatċ One further notes that ψ n = ψ 0 n for n ∈ I D from (S5), thus the instantaneous energies can be written as Hψ 0 n = E n ψ 0 n , n ∈ I D .
(S14) Figure S1. Schematic representation of the degenerate subspace D (red) with states ψ 0 m and the larger spaceD (orange) that includes states ψ 0 n that are nearly degenerate to the states ψ 0 m ∈ D. The full Hilbert space is denoted H.
Finally, by multiplying both sides of (S21) by √ ρ and integrating over k, one finds These are the coupled dynamics for the degrees of freedom z n in the Bloch wave packet. The argument for each quantity on the right hand side is simply k c .

NEAR DEGENERACY
The above results require the subspace D to be degenerate. Next we consider the near-degenerate case, where states in given Block band eigenstates with different wave vectors have close in energy eigenvalues. Schematically this situation can be given in Fig. S1, where the degenerate subspace D and the near degenerate subspaceD are shown. The eigenstates of the unperturbed Hamiltonian H 0 are ψ 0 n ∈ D and ψ 0 m ∈D. Since H 0 is nearly degenerate inD, then |E 0 n − E 0 m | < ψ 0 m |H 1 |ψ 0 n meaning that So, even in cases of near degeneracy, the projected EOM hold.

ELEMENTS IN THE EOM CALCULATIONS
The projected EOM in (5) and (6) of the main text within first order of the perturbation H 1 can now be used to obtain the time dynamics of the position x and the gauge invariant wave vector k = q − e A c found in (7) of the main text.
The wave packet subspace D is given by the Bloch states ψ n (k) with n ∈ I D , which is the natural basis for the calculations that follow. In this representation, the elements of position x and gauge invariant wave vector k, needed for P D (x) and P D (k), are found as where A n1n2 (k) = i u n1 (k)|∇|u n2 (k) . One also needs to consider P D (H 1 ) where H 1 = 1 2 ∂Hc ∂xc · (x − x c ) + (x − x c ) · ∂Hc ∂xc . We find that