Multiband effects in equations of motion of observables beyond the semiclassical approach

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–10]. This transport formalism has also been applied to light propagation in photonic materials, where the Berry phase effects have been shown experimentally [11].

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 electron wave packets in an electromagnetic field [1, 5, 12–16]. 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, 17–20]. 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 long-wavelength 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 [12–15]. 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, 12–14], which is in line with the assumption from semiclassical theory 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 [15] 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 further investigation, however. This is especially relevant for topological materials with spin textures, where interband effects can lead to non-Abelian gauge fields. 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 [14, 16] or a subset of degenerate bands as is the case in [5, 15]. 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, 14]. 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.

In the presence of electromagnetic fields, quantum mechanical processes can be described using the Hamiltonian $H{\left({\bf{p}}-e{\bf{A}}({\bf{x}}),{\bf{x}})=\tfrac{1}{2m}({\bf{p}}-e{\bf{A}}({\bf{x}})\right)}^{2}+V({\bf{x}})+e\varphi ({\bf{x}})$, 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 [21], involves semiclassical particles represented as wave packets [15]

$$\begin{eqnarray}&&\psi =\displaystyle \sum _{n}\int {{\rm{d}}}^{3}{\bf{q}}{c}_{n}({\bf{q}}){\psi }_{n}({\bf{q}}),{c}_{n}({\bf{q}})=\sqrt{\rho ({\bf{q}}-{{\bf{q}}}_{c})}{z}_{n}({\bf{q}}),\end{eqnarray} \tag{ 1 }$$

where ${\psi }_{n}({\bf{q}})={{\rm{e}}}^{{\rm{i}}{\bf{q}}\cdot {\bf{x}}}{u}_{n}({\bf{q}})$ are Bloch states with Bloch functions u_n(q) for wave vector q in the nth band with expansion coefficients c_n with ${z}_{n}={a}_{n}{{\rm{e}}}^{-{\rm{i}}{\gamma }_{n}({\bf{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 $\langle {\bf{q}}\rangle ={{\bf{q}}}_{c}$ in reciprocal space. The probability that the particle is found in the nth band is a_n².

This is a usual assumption for a single band wave packet, since its collapse in reciprocal space associated with the breakdown of its wave function localization is governed by the much larger Ehrenfest time [22] evolving at a scale $\mathrm{ln}(1/h)$ with h being the smaller of $\tfrac{{\hslash }T}{m{\lambda }^{2}}$ or $\tfrac{l}{\lambda }$ values (T—macroscopic time of the slowly varying fields, λ—wavelength of the fields, m—mass of the particle, l—lattice constant). For multiband wave packets, especially when there are crossing points in the band structure, this issue is currently under investigation. Recent studies have shown that the collapse time of a packet initially localized in a single band whose evolution allows interband transitions after encountering a crossing point is also relatively long and it is determined by the 1/h scale [23, 24]. Also, it has been shown that a multiband description of a wave packet on a generic honeycomb potential with Dirac points accurately approximates the full solution to the Schroedinger equation for large times even when the packet is centered at the crossing point [25]. Although the problem of the wave packet collapse in reciprocal space due to loss of localization needs further independent investigations, these recent works give evidence that the wave packet description is meaningful even in the case of multiband composition containing crossing band structure points.

In figure 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 $\epsilon =\tfrac{{\rm{\Delta }}x}{\lambda }\ll 1$. Neglecting the small wave packet spread in reciprocal space, one then obtains the well-known semiclassical EOM [21, 26–28]

$$\begin{eqnarray}&&\dot{\langle {\bf{x}}\rangle }=\displaystyle \frac{1}{{\hslash }}{\rm{\nabla }}{ \mathcal E }({{\bf{k}}}_{c}),\dot{\langle {\bf{k}}\rangle }=\displaystyle \frac{e}{{\hslash }}({\bf{E}}+\dot{\langle {\bf{x}}\rangle }\times {\bf{B}}),\end{eqnarray} \tag{ 2 }$$

where $\langle {\bf{x}}\rangle ={{\bf{x}}}_{c}$ and $\langle {\bf{k}}\rangle ={{\bf{k}}}_{c}$ are the expectation values of position and gauge invariant crystal momentum ${\bf{k}}={\bf{q}}-\tfrac{e}{{\hslash }}{\bf{A}}$, respectively. Also, $\tfrac{{\rm{d}}}{{\rm{d}}{t}}\langle {\bf{x}}\rangle =\dot{\langle {\bf{x}}\rangle }$ and $\tfrac{{\rm{d}}}{{\rm{d}}{t}}\langle {\bf{k}}\rangle =\dot{\langle {\bf{k}}\rangle }$ with the same notation holding for all operators. Here ${\rm{\nabla }}{ \mathcal E }({{\bf{k}}}_{c})$ is the reciprocal space gradient of the single band dispersion of the bare crystal Hamiltonian ${H}_{0}({\bf{p}},{\bf{x}})=\tfrac{{{\bf{p}}}^{2}}{2m}+V({\bf{x}})$ with ${H}_{0}({\bf{p}}+{\hslash }{\bf{q}},{\bf{x}}){u}_{n}({\bf{q}})={{ \mathcal E }}_{n}({\bf{q}}){u}_{n}({\bf{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 [21].

Zoom In Zoom Out Reset image size

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, 12–14]. 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 [14]. Retaining only up to the first order in this series, one finds

$$\begin{eqnarray}&&H\approx {H}_{c}({\bf{p}}-e{{\bf{A}}}_{c},{\bf{x}};{{\bf{x}}}_{c})+\displaystyle \frac{1}{2}\displaystyle \frac{\partial {H}_{c}}{\partial {{\bf{x}}}_{c}}\cdot ({\bf{x}}-{{\bf{x}}}_{c})+{\rm{h}}.{\rm{c}}.\end{eqnarray} \tag{ 3 }$$

Here ${H}_{c}({\bf{p}}-e{{\bf{A}}}_{c},{\bf{x}};{{\bf{x}}}_{c})={H}_{0}({\bf{p}}-e{{\bf{A}}}_{c},{\bf{x}})+e{\varphi }_{c}=\tfrac{1}{2m}{\left({\bf{p}}-e{{\bf{A}}}_{c}\right)}^{2}\,+\,V({\bf{x}})+e{\varphi }_{c},{{\bf{A}}}_{c}={\bf{A}}({{\bf{x}}}_{c})$ 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 $\epsilon =\tfrac{{\rm{\Delta }}x}{\lambda }$. This is realized by noting that $\tfrac{\partial {H}_{c}}{\partial {{\bf{x}}}_{c}}={\sum }_{j}\tfrac{\partial {H}_{c}}{\partial {A}_{c,j}}\tfrac{\partial {A}_{c,j}}{\partial {{\bf{x}}}_{c}}$ results in a 1/λ factor from the Fourier transformed vector potential (see (A1) in appendix A), while (x − x_c) determines the wave packet spread in real space Δx. Using a Lagrangian formalism within the semiclassical single band approximation [12–14], the EOM are found as

$$\begin{eqnarray}&&\dot{\langle {\bf{x}}\rangle }=\displaystyle \frac{1}{{\hslash }}{\rm{\nabla }}{{ \mathcal E }}_{M}({{\bf{k}}}_{c})-\dot{\langle {\bf{k}}\rangle }\times {\bar{{\boldsymbol{\Omega }}}}_{n},\dot{\langle {\bf{k}}\rangle }=\displaystyle \frac{e}{{\hslash }}({\bf{E}}+\dot{\langle {\bf{x}}\rangle }\times {\bf{B}}),\end{eqnarray} \tag{ 4 }$$

where ${\bar{{\boldsymbol{\Omega }}}}_{n}={\rm{\nabla }}\times {\bar{{\boldsymbol{ \mathcal A }}}}_{n}$ is the Berry curvature and ${\bar{{\boldsymbol{ \mathcal A }}}}_{n}={\rm{i}}\langle {u}_{n}({{\bf{k}}}_{c})| {\rm{\nabla }}| {u}_{n}({{\bf{k}}}_{c})\rangle$ is the Berry connection with Bloch functions u_n (k_c) for the nth band. The energy ${{ \mathcal E }}_{M}({{\bf{k}}}_{c})={ \mathcal E }({{\bf{k}}}_{c})-{\bar{{\bf{M}}}}_{n}\cdot {\bf{B}}$ is the lattice band dispersion with a contribution from a reciprocal space magnetic moment ${\bar{{\bf{M}}}}_{n}=-\tfrac{{\rm{i}}{e}}{2{\hslash }}\langle {\rm{\nabla }}{u}_{n}({{\bf{k}}}_{c})| \left({ \mathcal E }({{\bf{k}}}_{c})-{H}_{0}({\bf{p}}+{\hslash }{{\bf{k}}}_{c},{\bf{x}})\right)\,\times | {\rm{\nabla }}{u}_{n}({{\bf{k}}}_{c})\rangle$ 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 $\dot{\langle {\bf{x}}\rangle }$. 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 [15] which uses a projection procedure of operators onto a wave packet subspace consisting of a single energy band. The resulting operators from this projection lead to noncanonical commutation relations giving rise to the Berry phase contributions in $\dot{\langle {\bf{x}}\rangle }$. Another approach [16], 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 [15, 16, 29].

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 semiclassical expansion has been offered in [30]. Also, the majority of the works have focused on a single band projection with the exception of [5, 31] 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 packets with different dispersions must be resolved, especially for situations when the Fermi level is in close proximity to crossing points between bands.

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⁰(α(t)) + H¹(α(t)) where the unperturbed part H⁰ 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¹(α(t)) with perturbation parameter , (ii) α(t) is turned on and thereafter varies slowly in time, and (iii) the eigenfunctions of H⁰(α(t)) are ${\psi }_{n}^{0}(\alpha (t))$ with the initial eigenstates of the wave packet confined to a subspace whose eigenstates 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 ${\psi }_{n}^{0}\in 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

$$\begin{eqnarray}&&\dot{\langle O\rangle }=\displaystyle \frac{{\rm{i}}}{{\hslash }}\langle \left[{{ \mathcal P }}_{D}(H),{{ \mathcal P }}_{D}(O)\right]\rangle +\left\langle \displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{{ \mathcal P }}_{D}(O)\right\rangle ,\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{{ \mathcal P }}_{D}(O)=\displaystyle \sum _{m,n\in {I}_{D}}| {\psi }_{m}^{0}\rangle \langle {\psi }_{m}^{0}| O| {\psi }_{n}^{0}\rangle \langle {\psi }_{n}^{0}| ,\end{eqnarray} \tag{ 5b }$$

$$\begin{eqnarray}&&{\dot{c}}_{n}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}{c}_{n}\langle {\psi }_{n}^{0}| H| {\psi }_{n}^{0}\rangle -\displaystyle \sum _{m\in {I}_{D}}{c}_{m}\langle {\psi }_{n}^{0}| {\dot{\psi }}_{m}^{0}\rangle ,\end{eqnarray} \tag{ 5c }$$

where ${{ \mathcal P }}_{D}$ is the projection operator onto D. Here c_n are the coefficients for the expansion of the wave packet $\psi ={\sum }_{n\in {I}_{D}}{c}_{n}{\psi }_{n}^{0}$ 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 appendix A.

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 $\delta {{ \mathcal E }}_{{bw}}$, which is separated by a large energy gap ${\rm{\Delta }}{ \mathcal E }$ from the rest of the band structure ($\delta {{ \mathcal E }}_{{bw}}\ll {\rm{\Delta }}{ \mathcal E }$), as schematically illustrated in figure 2.

Zoom In Zoom Out Reset image size

As explained in [14] and appendix A, the series expansion of the electromagnetic Hamiltonian in (3) can be considered as having an unperturbed term coming from the bare crystal potential H⁰ = H_c(p − eA_c , x ; x_c) and a perturbation ${H}^{1}=\tfrac{1}{2}\tfrac{\partial {H}_{c}}{\partial {{\bf{x}}}_{c}}\cdot ({\bf{x}}-{{\bf{x}}}_{c})+{\rm{h}}.{\rm{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 $\delta {{ \mathcal E }}_{{bw}}$ of states near a Fermi level tuned to a crossing point is well separated from the rest of the band structure ($\delta {{ \mathcal E }}_{{bw}}\ll {\rm{\Delta }}{ \mathcal E }$), leading to condition (iii). In figure 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.

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 ${\psi }_{n}^{0}({\bf{q}})={{\rm{e}}}^{{\rm{i}}{\bf{q}}\cdot {\bf{x}}}{u}_{n}^{0}({\bf{q}}-\frac{e}{\hslash }{{\bf{A}}}_{c})$. For clarity and concreteness, we take the static electromagnetic circular gauge with ${\bf{A}}=\tfrac{1}{2}\left({\bf{B}}\times {\bf{x}}\right)$ and $\varphi =-{\bf{E}}\cdot {\bf{x}}$ and time-independent fields (although the general theory applies to fields slowly varying in time). After some calculations with extensive details given in appendix B, we find

$$\begin{eqnarray}\begin{array}{rcl}\dot{\langle {\bf{x}}\rangle } & = & \langle {\bf{v}}\rangle -\displaystyle \frac{e}{2{{\hslash }}^{2}}| {z}_{{n}_{1}}{| }^{2}\left\{\left[\left({{\bf{x}}}_{c}-{{\boldsymbol{\nu }}}_{{\gamma }_{{n}_{1}}}\right)\times {\bf{B}}\right]\cdot {\rm{\nabla }}\right\}{{\bf{v}}}_{{n}_{1}}-\displaystyle \frac{e}{{\hslash }}{z}_{{n}_{1}}^{* }{z}_{{n}_{2}}\left\{\displaystyle \frac{1}{4}\left[\left({\bf{B}}\times {\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{2}}\right)\cdot {\rm{\nabla }}\right]\left({{\bf{v}}}_{{n}_{1}}+{{\bf{v}}}_{{n}_{2}}\right)\right.\\ & & +\displaystyle \frac{{\rm{i}}}{4}\left[{\bf{B}}\cdot \left(\left({{\boldsymbol{\nu }}}_{{\gamma }_{{n}_{1}}}+{{\boldsymbol{\nu }}}_{{\gamma }_{{n}_{2}}}-{{\bf{x}}}_{c}\right)\times {\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{2}}\right)\right]\left({{\bf{v}}}_{{n}_{1}}-{{\bf{v}}}_{{n}_{2}}\right)+\displaystyle \frac{{\rm{i}}}{4}\left[{\bf{B}}\cdot \left(\left({{\boldsymbol{\nu }}}_{{\gamma }_{{n}_{1}}}+{{\boldsymbol{\nu }}}_{{\gamma }_{{n}_{2}}}-{{\bf{x}}}_{c}\right)\times \left({{\bf{v}}}_{{n}_{1}}-{{\bf{v}}}_{{n}_{2}}\right)\right)\right]{\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{2}}\\ & & +\displaystyle \frac{1}{4}\left[\left({{\bf{v}}}_{{n}_{1}}+{{\bf{v}}}_{{n}_{2}}\right)\times {\bf{B}}\right]\times {\bar{{\boldsymbol{\Omega }}}}_{{n}_{1}{n}_{2}}-\displaystyle \frac{1}{4}\left({{\bf{v}}}_{{n}_{3}}\times {\bf{B}}\right)\times {\rm{i}}{\bar{{\bf{Q}}}}_{{n}_{1}{n}_{3}{n}_{2}}+\left[\displaystyle \frac{1}{2}\left(\langle {\bf{v}}\rangle \times {\bf{B}}\right)+{\bf{E}}\right]\times \left({\bar{{\boldsymbol{\Omega }}}}_{{n}_{1}{n}_{2}}-{\rm{i}}{\bar{{\bf{Q}}}}_{{n}_{1}{n}_{3}{n}_{2}}\right)\\ & & \left.+\displaystyle \frac{1}{e}{\rm{\nabla }}\left({\bf{B}}\cdot {\bar{{\bf{M}}}}_{{n}_{1}{n}_{2}}\right)+{\rm{i}}{\bf{B}}\cdot \left(\displaystyle \frac{1}{4}{\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{3}}\times {{\bf{v}}}_{{n}_{1}}-\displaystyle \frac{1}{e}{\bar{{\bf{M}}}}_{{n}_{1}{n}_{3}}\right){\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{3}{n}_{2}}-{\rm{i}}{\bf{B}}\cdot \left(\displaystyle \frac{1}{4}{\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{3}{n}_{2}}\times {{\bf{v}}}_{{n}_{2}}-\displaystyle \frac{1}{e}{\bar{{\bf{M}}}}_{{n}_{3}{n}_{2}}\right){\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{3}}\right\},\end{array}\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&\dot{\langle {\bf{k}}\rangle }=\displaystyle \frac{e}{{\hslash }}\left({\bf{E}}+\langle {\bf{v}}\rangle \times {\bf{B}}\right),\end{eqnarray} \tag{ 6b }$$

$$\begin{eqnarray}\begin{array}{rcl}{\rm{i}}{\hslash }{\dot{z}}_{n} & = & \left\{{{ \mathcal E }}_{n}+e\left[\displaystyle \frac{1}{2}\left({{\bf{v}}}_{n}\times {\bf{B}}\right)+{\bf{E}}\right]\cdot {{\bf{x}}}_{c}\right\}{z}_{n}-\displaystyle \sum _{m\in {I}_{D}}\left\{\displaystyle \frac{e}{2}\left[\left(\displaystyle \frac{1}{2}\left({{\bf{v}}}_{n}+{{\bf{v}}}_{m}\right)+\dot{\langle {\bf{x}}\rangle }\right)\times {\bf{B}}+2{\bf{E}}\right]\right.\\ & & \cdot \,{\bar{{ \mathcal A }}}_{{nm}}+{\bf{B}}\left.\cdot {\bar{{\bf{M}}}}_{{nm}}\Space{0ex}{3.2ex}{0ex},\}\right\}{z}_{m}-\displaystyle \frac{{\rm{i}}{e}}{4}\left[\left({{\bf{v}}}_{n}-\langle {\bf{v}}\rangle \right)\times {\bf{B}}\right]\cdot {\rm{\nabla }}{z}_{n},\end{array}\end{eqnarray} \tag{ 6c }$$

where summation over all bands in the multiband wave packet for each band index n₁, n₂, and n₃ in each term in $\dot{\langle {\bf{x}}\rangle }$ and $\dot{\langle {\bf{k}}\rangle }$ is implied. We have defined the interband Berry connection ${\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{2}}={\rm{i}}\langle {u}_{{n}_{1}}({{\bf{k}}}_{c})| {\rm{\nabla }}| {u}_{{n}_{2}}({{\bf{k}}}_{c})\rangle$ and interband Berry curvature ${\bar{{\boldsymbol{\Omega }}}}_{{n}_{1}{n}_{2}}={\rm{\nabla }}\times {\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{2}}$, as well as an interband magnetization ${\bar{{\bf{M}}}}_{{n}_{1}{n}_{2}}\,=-\tfrac{{\rm{i}}{e}}{4{\hslash }}\left\langle {\rm{\nabla }}{u}_{{n}_{1}}({{\bf{k}}}_{c})\right|\left({{ \mathcal E}}_{{n}_{1}}({{\bf{k}}}_{c})+{{ \mathcal E }}_{{n}_{2}}({{\bf{k}}}_{c})-2{H}_{0}({\bf{p}}+{\hslash }{{\bf{k}}}_{c},{\bf{x}})\right)$ $\times \left.| {\rm{\nabla }}{u}_{{n}_{2}}({{\bf{k}}}_{c})\right\rangle$.The non-Abelian coupling between interband Berry connections is denoted as ${\bar{{\bf{Q}}}}_{{n}_{1}{n}_{2}{n}_{3}}={\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{2}}\times {\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{2}{n}_{3}}$. We note that ${{\bf{x}}}_{c}$ and k_c are expectation values calculated with the wave function $\psi ={\sum }_{n\in {I}_{D}}\int {{\rm{d}}}^{3}{\bf{k}}{c}_{n}({\bf{k}}){\psi }_{n}^{0}({\bf{k}})$ for the wave packet. The coefficients c_n are found from first solving for the z_n in (6) and then using (1). Also, $\langle {\bf{v}}\rangle ={\sum }_{n\in {I}_{D}}{a}_{n}^{2}{{\bf{v}}}_{n}$ depends only on the wave packet distribution across bands and their dispersions, ${{\bf{v}}}_{n}=\tfrac{1}{{\hslash }}{\rm{\nabla }}{{ \mathcal E }}_{n},{{\boldsymbol{\nu }}}_{{\gamma }_{n}}={\rm{\nabla }}{\gamma }_{n}$, and ∇ is understood as differentiation with respect to k_c everywhere. Since ${\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{2}}$,${\bar{{\boldsymbol{\Omega }}}}_{{n}_{1}{n}_{2}}$, and ${\bar{{\bf{M}}}}_{{n}_{1}{n}_{2}}$ are Hermitian , it is noted that (6a) is real.

A comparison with the single band EOM in (4) shows that the results for $\dot{\langle {\bf{x}}\rangle }$ 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 single band case, these dispersion terms reduce to the ${\rm{\nabla }}{ \mathcal E }({{\bf{k}}}_{c})$ contribution in $\dot{\langle {\bf{x}}\rangle }$ in (4). All terms in the third line and the terms containing ${{\bf{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 single band $\dot{\langle {\bf{k}}\rangle }\times {\bar{{\boldsymbol{\Omega }}}}_{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 the ${\bar{{\bf{M}}}}_{n}\cdot {\bf{B}}$ contribution to $\dot{\langle {\bf{x}}\rangle }$ in (4). Surprisingly, the only difference between the $\dot{\langle {\bf{k}}\rangle }$ in (6) and that in (4) is the $\langle {\bf{v}}\rangle$ that appears in (6) instead of the full expression for $\dot{\langle {\bf{x}}\rangle }$. The reason that $\dot{\langle {\bf{k}}\rangle }$ in (6) is similar to $\dot{\langle {\bf{k}}\rangle }$ in (4) is because k is diagonal in the zeroth order Bloch basis. Physically, the real space wave packet does not 'see' the distribution across bands in reciprocal space in regards to the momentum and we thus retain a Lorentz-like force expression for $\dot{\langle {\bf{k}}\rangle }$.

We examine (6) in some specific cases. Firstly, the single band wave packet EOM reduce to $\dot{\langle {\bf{x}}\rangle }\,=\tfrac{1}{{\hslash }}{\rm{\nabla }}{{ \mathcal E }}_{M}({{\bf{k}}}_{c})-\dot{\langle {\bf{k}}\rangle }\times \bar{{\boldsymbol{\Omega }}}$ and $\dot{\langle {\bf{k}}\rangle }=\tfrac{e}{{\hslash }}({\bf{E}}+\langle {\bf{v}}\rangle \times {\bf{B}})$, which is the same as (4), except for the ${\bf{v}}$ that appears in $\dot{{\bf{k}}}$ instead of $\dot{{\bf{x}}}$. We also consider the case of degenerate bands, i.e. ${{ \mathcal E }}_{n}={ \mathcal E }$ for all N bands, for which (6) reduces to

$$\begin{eqnarray}&&\dot{\langle {\bf{x}}\rangle }=\langle {\bf{v}}\rangle -\displaystyle \frac{1}{{\hslash }}{z}_{{n}_{1}}^{* }{z}_{{n}_{2}}\left[{\rm{\nabla }}\left({\bf{B}}\cdot {\bar{{\bf{M}}}}_{{n}_{1}{n}_{2}}\right)-{\rm{i}}{\bf{B}}\cdot {\bar{{\bf{M}}}}_{{n}_{1}{n}_{3}}{\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{3}{n}_{2}}+{\rm{i}}{\bf{B}}\cdot {\bar{{\bf{M}}}}_{{n}_{3}{n}_{2}}{\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{3}}-\dot{\langle {\bf{k}}\rangle }\times \left({\bar{{\boldsymbol{\Omega }}}}_{{n}_{1}{n}_{2}}-{\rm{i}}{\bar{{\bf{Q}}}}_{{n}_{1}{n}_{3}{n}_{2}}\right)\right],\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&\dot{\langle {\bf{k}}\rangle }=\displaystyle \frac{e}{{\hslash }}\left({\bf{E}}+\langle {\bf{v}}\rangle \times {\bf{B}}\right),\end{eqnarray} \tag{ 7b }$$

$$\begin{eqnarray}&&{\rm{i}}{\hslash }{\dot{z}}_{n}=\left\{{ \mathcal E }+e\left[\displaystyle \frac{1}{2{\hslash }}\left({\rm{\nabla }}{ \mathcal E }\times {\bf{B}}\right)+{\bf{E}}\right]\cdot {{\bf{x}}}_{c}\right\}{z}_{n}-\displaystyle \sum _{m\in {I}_{D}}\left\{\displaystyle \frac{e}{2}\left[\left(\displaystyle \frac{1}{{\hslash }}{\rm{\nabla }}{ \mathcal E }+\dot{\langle {\bf{x}}\rangle }\right)\times {\bf{B}}+2{\bf{E}}\right]\cdot {\bar{{ \mathcal A }}}_{{nm}}+{\bf{B}}\cdot {\bar{{\bf{M}}}}_{{nm}}\right\}{z}_{m}.\end{eqnarray} \tag{ 7c }$$

All notation and conventions for summation over band indices and ∇ differentiation are the same as in (6). Comparing the results in (7a) with the general expression in (6a) shows that the only terms that do not vanish in $\dot{\langle {\bf{x}}\rangle }$ are those coming from the coupling between the magnetic field and interband magnetizations. Also, the anomalous velocity now contains only terms coupling $\dot{\langle {\bf{k}}\rangle }$ to ${\bar{{\boldsymbol{\Omega }}}}_{{n}_{1}{n}_{2}}$ and ${\bar{{\bf{Q}}}}_{{n}_{1}{n}_{3}{n}_{2}}$, which are manifestations of the non-canonical commutation relations of the projected position operators. We also note that given the global degeneracy of the energy bands, the results in (7) are SU(N) gauge invariant similar to the results obtained in [5] as well as those in [15] when applied to the specific electromagnetic gauge taken here. Since our approach applies to systems with diverse energy dispersions that may not contain globally degenerate bands, an SU(N) gauge invariance is not expected in (6).

The general EOM in (6) can also be applied to a pair of linear bands ${{ \mathcal E }}_{\pm }({\bf{k}})=\pm {vk}$, where v is the isotropic Fermi velocity (schematics in figure 2). Such a band structure captures the vicinity of a band crossing in a Weyl semimetal or graphene [18, 19] and may not have any gauge symmetry. In this case, the opposite velocities lead to spreading of the packet wave function (1) into two different band components, thus the question about the collapse of the packet in real and reciprocal spaces needs to be examined. We note that as long as the semiclassical approximation is maintained (also required by semiclassical transport theory), such that the distance between the two components is much smaller than the fields wavelength, the wave packet description remains valid. Also, in a transport calculation, which includes scattering (not considered here), the time scale for spreading of the packet in real space must be larger than the scattering time. The collapse of the packet in reciprocal space is currently under investigation as noted earlier, although it has been shown that even for multiband wave packets with crossing points in the band structure, the collapse times are relatively long [22–25, 32]. Although this issue needs further separate investigations, the EOM for the pair of linear bands can be obtained as

$$\begin{eqnarray}\begin{array}{rcl}\dot{\langle {\bf{x}}\rangle } & = & \langle {\bf{v}}\rangle -\displaystyle \frac{1}{{\hslash }}| {z}_{{n}_{1}}{| }^{2}{\rm{\nabla }}\left({\bf{B}}\cdot {\bar{{\bf{M}}}}_{{n}_{1}{n}_{1}}\right)-\displaystyle \frac{e}{{{\hslash }}^{2}}{z}_{{n}_{1}}^{* }{z}_{{n}_{2}}\left\{\displaystyle \frac{{\rm{i}}}{4}\left[{\bf{B}}\cdot \left(\left({{\boldsymbol{\nu }}}_{{\gamma }_{{n}_{1}}}+{{\boldsymbol{\nu }}}_{{\gamma }_{{n}_{2}}}-{{\bf{x}}}_{c}\right)\times \left({{\bf{v}}}_{{n}_{1}}-{{\bf{v}}}_{{n}_{2}}\right)\right)\right]{\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{2}}\right.\\ & & \left.+\displaystyle \frac{{\rm{i}}}{4}\left[{\bf{B}}\cdot \left(\left({{\boldsymbol{\nu }}}_{{\gamma }_{{n}_{1}}}+{{\boldsymbol{\nu }}}_{{\gamma }_{{n}_{2}}}-{{\bf{x}}}_{c}\right)\times {\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{2}}\right)\right]\left({{\bf{v}}}_{{n}_{1}}-{{\bf{v}}}_{{n}_{2}}\right)+\displaystyle \frac{1}{4}\left[\left({{\bf{v}}}_{{n}_{1}}+{{\bf{v}}}_{{n}_{2}}\right)\times {\bf{B}}\right]\times {\bar{{\boldsymbol{\Omega }}}}_{{n}_{1}{n}_{2}}\right\},\end{array}\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&\dot{\langle {\bf{k}}\rangle }=\displaystyle \frac{e}{{\hslash }}\left({\bf{E}}+\langle {\bf{v}}\rangle \times {\bf{B}}\right),\end{eqnarray} \tag{ 8b }$$

$$\begin{eqnarray}\begin{array}{rcl}{\rm{i}}{\hslash }{\dot{z}}_{n} & = & \left[{{ \mathcal E }}_{n}+e\left(\displaystyle \frac{1}{2{\hslash }}\left({\rm{\nabla }}{{ \mathcal E }}_{n}\times {\bf{B}}\right)+{\bf{E}}\right)\cdot {{\bf{x}}}_{c}-\displaystyle \frac{e}{2{\hslash }}\left({\rm{\nabla }}{{ \mathcal E }}_{n}\times {\bf{B}}\right)\cdot {\bar{{ \mathcal A }}}_{{nn}}\right]{z}_{n}\\ & & -\displaystyle \sum _{m\in \pm }\left[\displaystyle \frac{e}{2}\left(\dot{\langle {\bf{x}}\rangle }\times {\bf{B}}+2{\bf{E}}\right)\cdot {\bar{{ \mathcal A }}}_{{nm}}+{\bf{B}}\cdot {\bar{{\bf{M}}}}_{{nm}}\right]{z}_{m}\\ & & -\displaystyle \frac{{\rm{i}}{e}}{4}\left[\left({{\bf{v}}}_{n}-\langle {\bf{v}}\rangle \right)\times {\bf{B}}\right]\cdot {\rm{\nabla }}{z}_{n},\end{array}\end{eqnarray} \tag{ 8c }$$

where {n₁, n₂} = ± and $\langle {\bf{v}}\rangle =\tfrac{v}{{\hslash }}\left({a}_{+}^{2}-{a}_{-}^{2}\right)\hat{k}$. 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 ${{\bf{v}}}_{{n}_{1}}-{{\bf{v}}}_{{n}_{2}}$ 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. There is also an explicit dependence on a_±, γ_± in the exponents and in the ${{\boldsymbol{\nu }}}_{{\gamma }_{\pm }}$, thus to obtain more specific results, detailed knowledge of these properties is needed, which is beyond the scope of this paper.

In this study, the time dynamics of the position and gauge invariant crystal momentum subject to electromagnetic fields in a periodic environment is considered using a general Hamiltonian constructed for a multiband wave packet. One of the main contributions here is showing that when the 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 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 packet affect transport however, one must utilize a multiband generalization of the Boltzmann equation [21, 33] 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 are thankful to Professor Qian Niu for discussions. LW acknowledges financial support from the US Department of Energy under grant No. DE-FG02-06ER46297. CT acknowledges financial support from the Deutsche Forschungsgemeinschaft, in part through Research Training Group GRK 1621 and Collaborative Research Center SFB 1143, project A04. TS also acknowledges support from the ERASMUS+ Programme of the European Union and the hospitality of TU Dresden.

Here we demonstrate in detail how the 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.

A.1. General model

A.2. Perturbation

It is noted that the Hamiltonian in the derivation of the EOM of operators in the main text involves the perturbation ${H}^{1}=\tfrac{1}{2}\tfrac{\partial {H}_{c}}{\partial {{\bf{x}}}_{c}}\cdot ({\bf{x}}-{{\bf{x}}}_{c})+{\rm{h}}{\rm{.}}{\rm{c}}.$, where the perturbation parameter is $\epsilon =\tfrac{{\rm{\Delta }}x}{\lambda }$ with Δx as the wave packet spread in real space and λ as the wavelength of the Fourier components of the vector potential. It is clear that Δx comes from the operator (x − x_c). To show that $\tfrac{\partial {H}_{c}}{\partial {{\bf{x}}}_{c}}$ is proportional to $\tfrac{1}{\lambda }$, we first simply observe that $\tfrac{\partial {H}_{c}}{\partial {{\bf{x}}}_{c}}={\sum }_{j}\tfrac{\partial {H}_{c}}{\partial {A}_{c,j}}\tfrac{\partial {A}_{c,j}}{\partial {{\bf{x}}}_{c}}$. Now we write A_c as a Fourier transform

$$\begin{eqnarray}&&{{\bf{A}}}_{{c}}=\int \tilde{{\bf{A}}}({\bf{q}}){{\rm{e}}}^{-{\rm{i}}{\bf{q}}\cdot {{\bf{x}}}_{{\bf{c}}}}{{\rm{d}}}^{3}{\bf{q}}.\end{eqnarray} \tag{ A1 }$$

Here $q=\tfrac{2\pi }{\lambda }$ 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, $\tilde{{\bf{A}}}({\bf{q}})$ is non-negligible only in the vicinity of q ≈ 0. Additionally, one finds for the spatial derivative $\tfrac{\partial {A}_{c,j}}{\partial {x}_{c,k}}={\rm{i}}\int {q}_{k}{\tilde{A}}_{j}({\bf{q}}){{\rm{e}}}^{-{\rm{i}}{\bf{q}}\cdot {{\bf{x}}}_{{\bf{c}}}}{{\rm{d}}}^{3}{\bf{q}}\propto \tfrac{1}{\lambda }$.

A.3. 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 (equation (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)

$$\begin{eqnarray}&&\psi (t)=\displaystyle \sum _{n\in {I}_{{ \mathcal H }}}{c}_{n}(t){\psi }_{n}(t),\end{eqnarray} \tag{ A2 }$$

where the sum is over the index set ${I}_{{ \mathcal H }}$ for a complete set of instantaneous eigenstates in the Hilbert space ${ \mathcal H }$.

For a time-dependent Hamiltonian, one can write the Schrödinger equation for the expansion coefficients

$$\begin{eqnarray}&&{\dot{c}}_{n}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}{c}_{n}{E}_{n}-\displaystyle \sum _{m\in {I}_{{ \mathcal H }}}{c}_{m}\langle {\psi }_{n}| {\dot{\psi }}_{m}\rangle .\end{eqnarray} \tag{ A3 }$$

The instantaneous eigenstates ψ_m are perturbatively expanded as

$$\begin{eqnarray}&&{\psi }_{m}={\psi }_{m}^{0}+\displaystyle \sum _{n\in {I}_{{ \mathcal H }}}\left({c}_{m,n}^{1}+{c}_{m,n}^{2}+...\right){\psi }_{n}^{0}.\end{eqnarray} \tag{ A4 }$$

One notes that each coefficient ${c}_{m,n}^{j}$ for j > 0 contains factors of the form $\tfrac{\langle {\psi }_{k}^{0}| {H}^{1}| {\psi }_{l}^{0}\rangle }{{E}_{l}^{0}-{E}_{k}^{0}}$ for some states ${\psi }_{k}^{0},{\psi }_{l}^{0}$ and energies ${E}_{k}^{0},{E}_{l}^{0}$. Since states in D are separated from all other states by a large energy gap, such terms are negligible when either ${\psi }_{k}^{0}\notin D$ and ${\psi }_{l}^{0}\in D$ or vice versa, and thus (using the index set I_D for the set of basis eigenstates that span D)

$$\begin{eqnarray}{\psi }_{m}=\left\{\begin{array}{ll}{\psi }_{m}^{0}+\displaystyle \sum _{n\notin D}\left({c}_{m,n}^{1}+{c}_{m,n}^{2}+...\right){\psi }_{n}^{0}, & m\notin {I}_{D}\\ {\psi }_{m}^{0}, & m\in {I}_{D}.\end{array}\right.\end{eqnarray} \tag{ A5 }$$

The above result implies that ψ_m is approximated as simply ${\psi }_{m}^{0}$ for $m\in {I}_{D}$ and is entirely contained outside of D when $m\notin {I}_{D}$. Next, the kth order of the expansion coefficients from (A3) are considered and can be cast as

$$\begin{eqnarray}&&{\dot{c}}_{n}^{k}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}\displaystyle \sum _{j=0}^{k}{c}_{n}^{k-j}{E}_{n}^{j}-\displaystyle \sum _{j=0}^{k-i}\displaystyle \sum _{i=0}^{k}\displaystyle \sum _{m\in {I}_{{ \mathcal H }}}{c}_{m}^{k-i-j}\langle {\psi }_{n}^{j}| {\dot{\psi }}_{m}^{i}\rangle .\end{eqnarray} \tag{ A6 }$$

It is realized that in the zeroth order term, ${\dot{c}}_{n}^{0}=-\tfrac{{\rm{i}}}{{\hslash }}{c}_{n}^{0}{E}_{n}^{0}-{\sum }_{m\in {I}_{{ \mathcal H }}}{c}_{m}^{0}\langle {\psi }_{n}^{0}| {\dot{\psi }}_{m}^{0}\rangle$, one can make use of the adiabatic approximation in the term containing $\langle {\psi }_{n}^{0}| {\dot{\psi }}_{m}^{0}\rangle$, giving

$$\begin{eqnarray}{\dot{c}}_{n}^{0}=\left\{\begin{array}{ll}-{c}_{n}^{0}\left(\tfrac{{\rm{i}}}{{\hslash }}{E}_{n}^{0}+\langle {\psi }_{n}^{0}| {\dot{\psi }}_{n}^{0}\rangle \right), & n\notin {I}_{D}\\ -\tfrac{{\rm{i}}}{{\hslash }}{c}_{n}^{0}{E}_{n}^{0}-\displaystyle \sum _{m\in {I}_{D}}{c}_{m}^{0}\langle {\psi }_{n}^{0}| {\dot{\psi }}_{m}^{0}\rangle , & n\in {I}_{D}.\end{array}\right.\end{eqnarray} \tag{ A7 }$$

From the above relation it is realized that since there is no perturbation at t = 0, then

$$\begin{eqnarray}&&{c}_{n}^{0}=0,n\notin {I}_{D}.\end{eqnarray} \tag{ A8 }$$

Consequently, there is no zeroth order contribution to the coefficients c_n for $n\notin {I}_{D}$ as expected from adiabatic theory. To zeroth order, the particle stays in the degenerate space D.

We now claim that ${c}_{n}\approx 0$ to all orders for $n\notin {I}_{D}$, which can be shown by induction. Suppose that ${c}_{n}^{j}=0$ for all $j\lt k$ when $n\notin {I}_{D}$. Thus using (A6), for $n\notin {I}_{D}$

$$\begin{eqnarray}&&{\dot{c}}_{n}^{k}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}{c}_{n}^{0}{E}_{n}^{0}-\displaystyle \sum _{m\notin {I}_{D}}{c}_{m}^{k}\langle {\psi }_{n}^{0}| {\dot{\psi }}_{m}^{0}\rangle -\displaystyle \sum _{j=0}^{k-i}\displaystyle \sum _{i=0}^{k}\displaystyle \sum _{m\in {I}_{D}}{c}_{m}^{k-i-j}\langle {\psi }_{n}^{j}| {\dot{\psi }}_{m}^{i}\rangle .\end{eqnarray} \tag{ A9 }$$

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 (A5), the terms in the second sum are negligible. Therefore, the non-negligible terms give

$$\begin{eqnarray}&&{\dot{c}}_{n}^{k}=-{c}_{n}^{k}\left(\displaystyle \frac{{\rm{i}}}{{\hslash }}{E}_{n}^{0}+\langle {\psi }_{n}^{0}| {\dot{\psi }}_{n}^{0}\rangle \right)\to {c}_{n}^{k}={c}_{n}^{k}(t=0){{\rm{e}}}^{-\tfrac{{\rm{i}}}{{\hslash }}{\int }_{0}^{t}{E}_{n}^{0}(t^{\prime} ){\rm{d}}{t}^{\prime} -{\int }_{0}^{t}\langle {\psi }_{n}^{0}(t^{\prime} )| {\dot{\psi }}_{n}^{0}(t^{\prime} )\rangle {\rm{d}}{t}^{\prime} },n\notin {I}_{D}.\end{eqnarray} \tag{ A10 }$$

However, there is no is perturbation at t = 0 and thus

$$\begin{eqnarray}&&{c}_{n}^{k}=0,n\notin {I}_{D}.\end{eqnarray} \tag{ A11 }$$

Consequently, there is no kth order contribution to the coefficients c_n for $n\notin {I}_{D}$. Coupled to the fact that ${c}_{n}^{0}\approx 0$ for $n\notin {I}_{D}$, induction shows that c_n = 0 to all orders for $n\notin {I}_{D}$. Therefore, (A3) transforms to

$$\begin{eqnarray}&&{\dot{c}}_{n}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}{c}_{n}{E}_{n}-\displaystyle \sum _{m\in {I}_{D}}{c}_{m}\langle {\psi }_{n}| {\dot{\psi }}_{m}\rangle ,\end{eqnarray} \tag{ A12 }$$

where the summation now runs over I_D instead of ${I}_{{ \mathcal H }}$. We are interested in (A12) when $n\in {I}_{D}$. By using the large energy gap condition and (A5), (A12), it is found that

$$\begin{eqnarray}&&{\dot{c}}_{n}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}{c}_{n}{E}_{n}-\displaystyle \sum _{m\in {I}_{D}}{c}_{m}\langle {\psi }_{n}^{0}| {\dot{\psi }}_{m}^{0}\rangle ,n\in {I}_{D}.\end{eqnarray} \tag{ A13 }$$

One further notes that ${\psi }_{n}={\psi }_{n}^{0}$ for n ∈ I_D from (A5), thus the instantaneous energies can be written as

$$\begin{eqnarray}&&H{\psi }_{n}^{0}={E}_{n}{\psi }_{n}^{0},n\in {I}_{D}.\end{eqnarray} \tag{ A14 }$$

We arrive at the important result for the expansion coefficients in (A13)

$$\begin{eqnarray}&&{\dot{c}}_{n}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}{c}_{n}\langle {\psi }_{n}^{0}| H| {\psi }_{n}^{0}\rangle -\displaystyle \sum _{m\in {I}_{D}}{c}_{m}\langle {\psi }_{n}^{0}| {\dot{\psi }}_{m}^{0}\rangle ,n\in {I}_{D}.\end{eqnarray} \tag{ A15 }$$

A.4. Dynamics of an observable

We now turn to the dynamics of an observable O

$$\begin{eqnarray}&&\dot{\langle O\rangle }=\displaystyle \sum _{m,n\in {I}_{{ \mathcal H }}}\left({\dot{c}}_{m}^{* }{c}_{n}+{c}_{m}^{* }{\dot{c}}_{n}\right)\langle {\psi }_{m}| O| {\psi }_{n}\rangle +{c}_{m}^{* }{c}_{n}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\langle {\psi }_{m}| O| {\psi }_{n}\rangle .\end{eqnarray} \tag{ A16 }$$

Due to the approximation that c_n = 0 to all orders for $n\notin {I}_{D}$ and (A5), (A16) can be approximated as

$$\begin{eqnarray}&&\dot{\langle O\rangle }=\displaystyle \sum _{m,n\in {I}_{D}}\left({\dot{c}}_{m}^{* }{c}_{n}+{c}_{m}^{* }{\dot{c}}_{n}\right)\langle {\psi }_{m}^{0}| O| {\psi }_{n}^{0}\rangle +{c}_{m}^{* }{c}_{n}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\langle {\psi }_{m}^{0}| O| {\psi }_{n}^{0}\rangle .\end{eqnarray} \tag{ A17 }$$

If we use the results in (A15), we find

$$\begin{eqnarray}\begin{array}{rcl}\dot{\langle O\rangle } & = & \displaystyle \frac{{\rm{i}}}{{\hslash }}\displaystyle \sum _{m,n\in {I}_{D}}{c}_{m}^{* }{c}_{n}\left(\langle {\psi }_{m}^{0}| H| {\psi }_{m}^{0}\rangle \langle {\psi }_{m}^{0}| O| {\psi }_{n}^{0}\rangle -\langle {\psi }_{m}^{0}| O| {\psi }_{n}^{0}\rangle \langle {\psi }_{n}^{0}| H| {\psi }_{n}^{0}\rangle \right)\\ & & +\displaystyle \sum _{m,n\in {I}_{D}}{c}_{m}^{* }{c}_{n}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\langle {\psi }_{m}^{0}| O| {\psi }_{n}^{0}\rangle -\displaystyle \sum _{k,m,n\in {I}_{D}}\left({c}_{k}\langle {\psi }_{n}^{0}| {\dot{\psi }}_{k}^{0}\rangle +{c}_{k}^{* }\langle {\dot{\psi }}_{k}^{0}| {\psi }_{m}^{0}\rangle \right)\langle {\psi }_{m}^{0}| O| {\psi }_{n}^{0}\rangle .\end{array}\end{eqnarray} \tag{ A18 }$$

Since $H={H}^{0}+{H}^{1}$ is diagonal in D using the basis $\left\{{\psi }_{n}^{0}\right\}$, we can rewrite (A18) as

$$\begin{eqnarray}\begin{array}{rcl}\dot{\langle O\rangle } & = & \displaystyle \frac{{\rm{i}}}{{\hslash }}\displaystyle \sum _{k,m,n\in {I}_{D}}{c}_{m}^{* }{c}_{n}\left(\langle {\psi }_{m}^{0}| H| {\psi }_{k}^{0}\rangle \langle {\psi }_{k}^{0}| O| {\psi }_{n}^{0}\rangle -\langle {\psi }_{m}^{0}| O| {\psi }_{k}^{0}\rangle \langle {\psi }_{k}^{0}| H| {\psi }_{n}^{0}\rangle \right)\\ & & +\displaystyle \sum _{m,n\in {I}_{D}}{c}_{m}^{* }{c}_{n}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\langle {\psi }_{m}^{0}| O| {\psi }_{n}^{0}\rangle +\displaystyle \sum _{k,m,n\in {I}_{D}}\left({c}_{k}\langle {\dot{\psi }}_{n}^{0}| {\psi }_{k}^{0}\rangle +{c}_{k}^{* }\langle {\psi }_{k}^{0}| {\dot{\psi }}_{m}^{0}\rangle \right)\langle {\psi }_{m}^{0}| O| {\psi }_{n}^{0}\rangle .\end{array}\end{eqnarray} \tag{ A19 }$$

Here we also moved the time derivatives on the states in the last two terms in (A18).

Finally, (A19) can be cast as

$$\begin{eqnarray}&&\dot{\langle O\rangle }=\displaystyle \frac{{\rm{i}}}{{\hslash }}\left[{{ \mathcal P }}_{D}(H),{{ \mathcal P }}_{D}(O)\right]+\left\langle \displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{{ \mathcal P }}_{D}(O)\right\rangle ,\end{eqnarray} \tag{ A20a }$$

$$\begin{eqnarray}&&{{ \mathcal P }}_{D}(O)=\displaystyle \sum _{m,n\in {I}_{D}}| {\psi }_{m}^{0}\rangle \langle {\psi }_{m}^{0}| O| {\psi }_{n}^{0}\rangle \langle {\psi }_{n}^{0}| .\end{eqnarray} \tag{ A20b }$$

These are the projected EOM with $\psi ={\sum }_{n\in {I}_{D}}{c}_{n}(t){\psi }_{n}^{0}$ and c_n (t) are the solutions to (A15). The essence of the projection of operators is that if a particle starts out in the subspace D, then the EOM take a familiar Heisenberg form with operators projected onto D. It is this projection of operators that gives rise to (5) and (6) in the main text, leading to altered commutation relations for the operators time dynamics.

A.5. Dynamics of the coefficients

We wish to write (A15) for the case of a Bloch wave packet in D, for which $\psi ={\sum }_{n\in {I}_{D}}\int {{\rm{d}}}^{3}{\bf{q}}{c}_{n}({\bf{q}}){\psi }_{n}({\bf{q}})$, with ${\psi }_{n}^{0}({\bf{q}})={{\rm{e}}}^{{\rm{i}}{\bf{q}}\cdot {\bf{x}}}{u}_{n}^{0}({\bf{q}}-\frac{e}{\hslash }{{\bf{A}}}_{c})$ and ${c}_{n}({\bf{q}})=\sqrt{\rho ({\bf{q}}-{{\bf{q}}}_{c})}{z}_{n}({\bf{q}})$. With this form for the wave packet, (A15) becomes

$$\begin{eqnarray}&&{\dot{c}}_{n}({\bf{k}})=-\displaystyle \frac{{\rm{i}}}{{\hslash }}\displaystyle \sum _{m\in {I}_{D}}\int {{\rm{d}}}^{3}{\bf{q}}{c}_{m}({\bf{q}})\langle {\psi }_{n}^{0}({\bf{k}})| H| {\psi }_{m}^{0}({\bf{q}})\rangle -\displaystyle \sum _{m\in {I}_{D}}{c}_{m}({\bf{k}})\langle {\psi }_{n}^{0}({\bf{k}})| {\dot{\psi }}_{m}^{0}({\bf{k}})\rangle ,n\in {I}_{D}.\end{eqnarray} \tag{ A21 }$$

In addition to the diagonal unperturbed elements of the Hamiltonian, using (B4) leads to

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{m\in {I}_{D}}\displaystyle \int {{\rm{d}}}^{3}{\bf{q}}{c}_{m}({\bf{q}})\langle {\psi }_{n}^{0}({\bf{k}})| H| {\psi }_{m}^{0}({\bf{q}})\rangle & = & {\rm{i}}{{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{{ \mathcal E }}_{n}\cdot {{\rm{\nabla }}}_{{\bf{k}}}{c}_{n}+\displaystyle \frac{{\rm{i}}}{2}{c}_{n}{{\rm{\nabla }}}_{{\bf{k}}}\cdot {{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{{ \mathcal E }}_{n}+{c}_{n}\left({{ \mathcal E }}_{n}-{{\bf{x}}}_{c}\cdot {{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{{ \mathcal E }}_{n}\right)\\ & & +{c}_{m}\left\{{{ \mathcal E }}_{{mn}}^{M}+\displaystyle \frac{{\rm{i}}}{2}{{\rm{\nabla }}}_{{{\bf{x}}}_{c}}\left({{ \mathcal E }}_{m}+{{ \mathcal E }}_{n}\right)\cdot {{\boldsymbol{ \mathcal A }}}_{{mn}}\right\}.\end{array}\end{eqnarray} \tag{ A22 }$$

Finally, by multiplying both sides of (A21) by $\sqrt{\rho }$ and integrating over ${\bf{k}}$, one finds

$$\begin{eqnarray}\begin{array}{rcl}{\rm{i}}{\hslash }{\dot{z}}_{n}({{\bf{k}}}_{c}) & = & \left({{ \mathcal E }}_{n}-{{\bf{x}}}_{c}\cdot {{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{{ \mathcal E }}_{n}\right){z}_{n}+\displaystyle \frac{{\rm{i}}}{2}\left({{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{{ \mathcal E }}_{n}-\displaystyle \sum _{m\in D}| {z}_{m}{| }^{2}{{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{{ \mathcal E }}_{m}\right){{\rm{\nabla }}}_{{{\bf{k}}}_{c}}{z}_{n}\\ & & +\displaystyle \sum _{m\in {I}_{D}}\left\{{{ \mathcal E }}_{{mn}}^{M}+\displaystyle \frac{{\rm{i}}}{2}{{\rm{\nabla }}}_{{{\bf{x}}}_{c}}\left({{ \mathcal E }}_{m}+{{ \mathcal E }}_{n}\right)\cdot {{\boldsymbol{ \mathcal A }}}_{{mn}}-{\rm{i}}{\hslash }\langle {u}_{n}^{0}| {\dot{u}}_{m}^{0}\rangle \right\}{z}_{m}.\end{array}\end{eqnarray} \tag{ A23 }$$

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 ${{\bf{k}}}_{c}$.

A.6. Near degeneracy

The above results require the subspace D to be degenerate. Next we consider the near degenerate case, where there are states outside of D whose energy eigenvalues are close to those in D. Schematically, this situation can be given in figure A1, where the degenerate subspace D and the near degenerate subspace $\tilde{D}$ are shown. The eigenstates of the unperturbed Hamiltonian H⁰ are ${\psi }_{n}^{0}\in D$ and ${\psi }_{m}^{0}\in \tilde{D}$. Since H⁰ is nearly degenerate in $\tilde{D}$, then $| {E}_{n}^{0}-{E}_{m}^{0}| \lt \langle {\psi }_{m}^{0}| {H}^{1}| {\psi }_{n}^{0}\rangle$ meaning that $| {E}_{n}^{0}-{E}_{m}^{0}|$ is small with respect to H¹.

Zoom In Zoom Out Reset image size

The full Hamiltonian is represented as $H={H}^{0}+{H}^{1}={\bar{H}}^{0}+{\rm{\Delta }}\bar{H}+{H}^{1}$. Here ${\bar{H}}^{0}$ has a degenerate spectrum, ${\bar{H}}^{0}{\psi }_{m}^{0}={\bar{E}}^{0}{\psi }_{m}^{0}$ for ${\psi }_{m}^{0}\in \tilde{D}$, where ${\bar{E}}^{0}$ is an average of the energies for $\tilde{D}$. Also, ${\bar{H}}^{0}{\psi }_{m}^{0}={H}^{0}{\psi }_{m}^{0}\,={E}_{m}^{0}{\psi }_{m}^{0}$ for ${\psi }_{m}^{0}\notin \tilde{D}$. We then have, $[{\bar{H}}^{0},{H}^{0}]=0$, and these operators are diagonal with matrix elements $\langle {\psi }_{m}^{0}| {\bar{H}}^{0}| {\psi }_{n}^{0}\rangle ={\bar{E}}^{0}{\delta }_{{mn}}$ for ${\psi }_{m}^{0},{\psi }_{n}^{0}\in D$ and $\langle {\psi }_{m}^{0}| {\bar{H}}^{0}| {\psi }_{n}^{0}\rangle ={E}_{n}^{0}{\delta }_{{mn}}$ otherwise. The diagonal Hermitian operator ${\rm{\Delta }}\bar{H}={H}^{0}-{\bar{H}}^{0}$, with matrix elements $\langle {\psi }_{m}^{0}| {\rm{\Delta }}\bar{H}| {\psi }_{n}^{0}\rangle =\left({E}_{n}^{0}-{\bar{E}}^{0}\right){\delta }_{{mn}}$. Essentially, the full Hamiltonian H, containing a degenerate in $\tilde{D}$ unperturbed ${\bar{H}}^{0}$ and ${\rm{\Delta }}\bar{H}+{H}^{1}$, is treated as a perturbation. However, for the degenerate ${\bar{H}}^{0}$, we can make use of the projected EOM from (A20) giving

$$\begin{eqnarray}&&\dot{\langle O\rangle }=\displaystyle \frac{{\rm{i}}}{{\hslash }}\left[{{ \mathcal P }}_{D}(H-{\bar{H}}^{0}),{{ \mathcal P }}_{D}(O)\right]+\left\langle \displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{{ \mathcal P }}_{D}(O)\right\rangle .\end{eqnarray} \tag{ A24 }$$

However, ${{ \mathcal P }}_{D}({\bar{H}}^{0})$ is proportional to the identity matrix in D and hence commutes with any other operator in D, so one obtains

$$\begin{eqnarray}&&\dot{\langle O\rangle }=\displaystyle \frac{{\rm{i}}}{{\hslash }}\left[{{ \mathcal P }}_{D}(H),{{ \mathcal P }}_{D}(O)\right]+\left\langle \displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{{ \mathcal P }}_{D}(O)\right\rangle .\end{eqnarray} \tag{ A25 }$$

So, even in cases of near degeneracy, the projected EOM hold.

The projected EOM in (5) and (6) of the main text within first order of the perturbation H¹ can now be used to obtain the time dynamics of the position x and the gauge invariant wave vector ${\bf{k}}={\bf{q}}-\tfrac{e}{{\hslash }}{{\bf{A}}}_{c}$.

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 ${{ \mathcal P }}_{D}({\bf{x}})$ and ${{ \mathcal P }}_{D}({\bf{k}})$, are found as

$$\begin{eqnarray}&&\langle {\psi }_{{n}_{1}}({\bf{k}})| {\bf{x}}| {\psi }_{{n}_{2}}({\bf{k}}^{\prime} )\rangle ={\delta }_{{n}_{1}{n}_{2}}\left({\rm{i}}{{\rm{\nabla }}}_{{\bf{k}}}+{{\boldsymbol{ \mathcal A }}}_{{n}_{1}{n}_{2}}\left({\bf{k}}\right)\right)\delta \left({\bf{k}}-{\bf{k}}^{\prime} \right),\end{eqnarray} \tag{ B1a }$$

$$\begin{eqnarray}&&\langle {\psi }_{{n}_{1}}({\bf{k}})| {\bf{k}}| {\psi }_{{n}_{2}}({\bf{k}}^{\prime} )\rangle ={\delta }_{{n}_{1}{n}_{2}}{\bf{k}}\delta \left({\bf{k}}-{\bf{k}}^{\prime} \right),\end{eqnarray} \tag{ B1b }$$

where ${{\boldsymbol{ \mathcal A }}}_{{n}_{1}{n}_{2}}({\bf{k}})={\rm{i}}\left\langle {u}_{{n}_{1}}({\bf{k}})\right|{\rm{\nabla }}| {u}_{{n}_{2}}({\bf{k}})\rangle$. One also needs to consider ${{ \mathcal P }}_{D}({H}^{1})$ where ${H}^{1}=\tfrac{1}{2}\left(\tfrac{\partial {H}_{c}}{\partial {{\bf{x}}}_{c}}\right.\cdot ({\bf{x}}-{{\bf{x}}}_{c})\,+({\bf{x}}-{{\bf{x}}}_{c})\cdot \left.\tfrac{\partial {H}_{c}}{\partial {{\bf{x}}}_{c}}\right)$. We find that

$$\begin{eqnarray}&&\left\langle {\psi }_{{n}_{1}}({\bf{k}})\right|\displaystyle \frac{\partial {H}_{c}}{\partial {{\bf{x}}}_{c}}\left|{\psi }_{{n}_{2}}({\bf{k}}^{\prime} )\right\rangle =\delta \left({\bf{k}}-{\bf{k}}^{\prime} \right)\left[{\delta }_{{n}_{1}{n}_{2}}{{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{{ \mathcal E }}_{{n}_{1}}\left({\bf{k}}\right)-\left({{ \mathcal E }}_{{n}_{1}}\left({\bf{k}}\right)-{{ \mathcal E }}_{{n}_{2}}\left({\bf{k}}^{\prime} \right)\right)\langle {u}_{{n}_{1}}\left({\bf{k}}\right)| {{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{u}_{{n}_{2}}\left({\bf{k}}^{\prime} \right)\rangle \right].\end{eqnarray} \tag{ B2 }$$

With these results, we can then obtain

$$\begin{eqnarray}&&\begin{array}{c}\left. \langle {\psi }_{{n}_{1}}({\bf{k}})\right|\displaystyle \frac{{\rm{\partial }}{H}_{c}}{{\rm{\partial }}{{\bf{x}}}_{c}}\cdot {\bf{x}}\left|{\psi }_{{n}_{2}}({{\bf{k}}}^{{\rm{{\prime} }}})\right. \rangle ={{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{{ \mathcal E }}_{{n}_{1}}\cdot \langle {\psi }_{{n}_{1}}({\bf{k}})| {\bf{x}}| {\psi }_{{n}_{2}}({{\bf{k}}}^{{\rm{{\prime} }}}) \rangle +{\rm{i}}{{\rm{\nabla }}}_{{{\bf{k}}}^{{\rm{{\prime} }}}}\cdot \left[\left({{ \mathcal E }}_{{n}_{1}}\left({\bf{k}}\right)-{{ \mathcal E }}_{{n}_{2}}\left({{\bf{k}}}^{{\rm{{\prime} }}}\right)\right) \langle {u}_{{n}_{1}}\left({\bf{k}}\right)| {{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{u}_{{n}_{2}}\left({{\bf{k}}}^{{\rm{{\prime} }}}\right) \rangle \delta ({\bf{k}}-{{\bf{k}}}^{{\rm{{\prime} }}})\right]\\ \,+\,{\rm{i}} \langle {{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{u}_{{n}_{1}}\left({\bf{k}}\right)\left|\cdot \left({{ \mathcal E }}_{{n}_{1}}\left({\bf{k}}\right)-{H}_{c}\left({\bf{k}}\right)\right)\right|{{\rm{\nabla }}}_{{{\bf{k}}}^{{\rm{{\prime} }}}}{u}_{{n}_{2}}\left({{\bf{k}}}^{{\rm{{\prime} }}}\right) \rangle \delta \left({\bf{k}}-{{\bf{k}}}^{{\rm{{\prime} }}}\right).\end{array}\end{eqnarray} \tag{ B3 }$$

Finally, the elements of the perturbation H¹ are calculated as

$$\begin{eqnarray}\begin{array}{rcl}\langle {\psi }_{{n}_{1}}({\bf{k}})\left|{H}^{1}\right|{\psi }_{{n}_{2}}({\bf{k}}^{\prime} )\rangle & = & \displaystyle \frac{1}{2}{{\rm{\nabla }}}_{{{\bf{x}}}_{c}}\left({{ \mathcal E }}_{{n}_{1}}\left({\bf{k}}\right)+{{ \mathcal E }}_{{n}_{2}}\left({\bf{k}}^{\prime} \right)\right)\cdot \langle {\psi }_{{n}_{1}}({\bf{k}})| {\bf{x}}| {\psi }_{{n}_{2}}({\bf{k}}^{\prime} )\rangle \\ & & -\displaystyle \frac{{\rm{i}}}{2}\left({{\rm{\nabla }}}_{{\bf{k}}}-{{\rm{\nabla }}}_{{\bf{k}}^{\prime} }\right)\cdot \left[\left({{ \mathcal E }}_{{n}_{1}}\left({\bf{k}}\right)-{{ \mathcal E }}_{{n}_{2}}\left({\bf{k}}^{\prime} \right)\right)\langle {u}_{{n}_{1}}\left({\bf{k}}\right)| {{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{u}_{{n}_{2}}\left({\bf{k}}^{\prime} \right)\rangle \delta \left({\bf{k}}-{\bf{k}}^{\prime} \right)\right]\\ & & +{{ \mathcal E }}_{{n}_{1}{n}_{2}}^{M}\delta \left({\bf{k}}-{\bf{k}}^{\prime} \right)-{{\bf{x}}}_{c}\cdot \left\langle {\psi }_{{n}_{1}}({\bf{k}})\left|\displaystyle \frac{\partial {H}_{c}}{\partial {{\bf{x}}}_{c}}\right|{\psi }_{{n}_{2}}({\bf{k}}^{\prime} )\right\rangle .\end{array}\end{eqnarray} \tag{ B4 }$$

Here the interband magnetization contribution to (B4) is

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal E }}_{{n}_{1}{n}_{2}}^{M} & = & \displaystyle \frac{{\rm{i}}}{2}\left\{\langle {{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{u}_{{n}_{1}}\left({\bf{k}}\right)\left|\cdot \left({{ \mathcal E }}_{{n}_{1}}\left({\bf{k}}\right)-{H}_{c}\left({\bf{k}}\right)\right)\right|{{\rm{\nabla }}}_{{\bf{k}}^{\prime} }{u}_{{n}_{2}}\left({\bf{k}}^{\prime} \right)\rangle \right.\\ & & \left.-\langle {{\rm{\nabla }}}_{{\bf{k}}}{u}_{{n}_{1}}\left({\bf{k}}\right)\left|\cdot \left({{ \mathcal E }}_{{n}_{2}}\left({\bf{k}}^{\prime} \right)-{H}_{c}\left({\bf{k}}\right)\right)\right|{{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{u}_{{n}_{2}}\left({\bf{k}}^{\prime} \right)\rangle \right\}.\end{array}\end{eqnarray} \tag{ B5 }$$

At this point, let us note that if a function f(k) depends on k, then ${{\rm{\nabla }}}_{{{\bf{x}}}_{c}}f=-\tfrac{e}{{\hslash }}{\sum }_{j}{{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{({{\bf{A}}}_{c})}_{j}\tfrac{\partial }{\partial {{\bf{k}}}_{j}}f$ and $\tfrac{{\rm{d}}}{{\rm{d}}{t}}f=-\tfrac{e}{{\hslash }}{\dot{{\bf{A}}}}_{c}\cdot {{\rm{\nabla }}}_{{\bf{k}}}f$. This observation applies to quantities like ${u}_{{n}_{2}}({\bf{k}})$ and ${{ \mathcal E }}_{{n}_{2}}({\bf{k}})$, for example. Therefore, (B4) becomes

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal E }}_{{n}_{1}{n}_{2}}^{M} & = & -\displaystyle \frac{{\rm{i}}{e}}{2{\hslash }}\displaystyle \sum _{j}{{\rm{\nabla }}}_{{{\bf{x}}}_{c}}{\left({{\bf{A}}}_{c}\right)}_{j}\left\{\langle {{\rm{\nabla }}}_{{{\bf{k}}}_{j}}{u}_{{n}_{1}}\left({\bf{k}}\right)\left|\cdot \left({{ \mathcal E }}_{{n}_{1}}\left({\bf{k}}\right)-{H}_{c}\left({\bf{k}}\right)\right)\right|{{\rm{\nabla }}}_{{\bf{k}}^{\prime} }{u}_{{n}_{2}}\left({\bf{k}}^{\prime} \right)\rangle \right.\\ & & \left.-\langle {{\rm{\nabla }}}_{{\bf{k}}}{u}_{{n}_{1}}\left({\bf{k}}\right)\left|\cdot \left({{ \mathcal E }}_{{n}_{2}}\left({\bf{k}}^{\prime} \right)-{H}_{c}\left({\bf{k}}\right)\right)\right|{{\rm{\nabla }}}_{{{\bf{k}}}_{j}^{{\prime} }}{u}_{{n}_{2}}\left({\bf{k}}^{\prime} \right)\rangle \right\}.\end{array}\end{eqnarray} \tag{ B6 }$$

In the electromagnetic circular gauge ${\bf{A}}=\tfrac{1}{2}\left({\bf{B}}\times {\bf{x}}\right)$, and when ${\bf{k}}={\bf{k}}^{\prime}$, the above reduces to

$$\begin{eqnarray}&&{{ \mathcal E }}_{{n}_{1}{n}_{2}}^{M}={\bf{B}}\cdot {{\bf{M}}}_{{n}_{1}{n}_{2}}({\bf{k}})=-\displaystyle \frac{{\rm{i}}{e}}{4{\hslash }}{\bf{B}}\cdot \langle {\rm{\nabla }}{u}_{{n}_{1}}({\bf{k}})| \left({{ \mathcal E }}_{{n}_{1}}+{{ \mathcal E }}_{{n}_{2}}-2{H}_{c}({\bf{k}})\right)\times | {\rm{\nabla }}{u}_{{n}_{2}}({\bf{k}})\rangle ,\end{eqnarray} \tag{ B7 }$$

where we have defined the interband magnetization ${{\bf{M}}}_{{n}_{1}{n}_{2}}({\bf{k}})=-\tfrac{{\rm{i}}{e}}{4{\hslash }}\langle {\rm{\nabla }}{u}_{{n}_{1}}({\bf{k}})| \left({{ \mathcal E }}_{{n}_{1}}+{{ \mathcal E }}_{{n}_{2}}-2{H}_{c}({\bf{k}})\right)\,\times | {\rm{\nabla }}{u}_{{n}_{2}}({\bf{k}})\rangle$, similar to what is found in the main text after (6). The other interband quantities in (6) of the main text are obtained by considering the following commutator, which can be calculated using (B1)

$$\begin{eqnarray}&&\left[{{ \mathcal P }}_{D}({{\bf{x}}}_{i}),{{ \mathcal P }}_{D}({{\bf{x}}}_{j})\right]={\rm{i}}\displaystyle \sum _{{n}_{1},{n}_{2}\in D}\int {{\rm{d}}}^{3}{\bf{k}}\left\{{\left({{\rm{\Omega }}}_{{n}_{1}{n}_{2}}\right)}_{{ij}}({\bf{k}})-{\rm{i}}\displaystyle \sum _{{n}_{3}\in {I}_{D}}{\left({Q}_{{n}_{1}{n}_{3}{n}_{2}}\right)}_{{ij}}({\bf{k}})\right\}\left|{\psi }_{{n}_{1}}({\bf{k}})\rangle \langle {\psi }_{{n}_{2}}({\bf{k}})\right|,\end{eqnarray} \tag{ B8 }$$

where

$$\begin{eqnarray}&&{\left({{\rm{\Omega }}}_{{n}_{1}{n}_{2}}\right)}_{{ij}}({\bf{k}})=\displaystyle \frac{\partial }{\partial {k}_{j}}{\left({{\boldsymbol{ \mathcal A }}}_{{n}_{1}{n}_{2}}\right)}_{i}({\bf{k}})-\displaystyle \frac{\partial }{\partial {k}_{i}}{\left({{\boldsymbol{ \mathcal A }}}_{{n}_{1}{n}_{2}}\right)}_{j}({\bf{k}})={\rm{i}}\left(\displaystyle \frac{\partial }{\partial {k}_{j}}\left\langle {u}_{{n}_{1}}({\bf{k}})\right|\displaystyle \frac{\partial }{\partial {k}_{i}}{u}_{{n}_{2}}({\bf{k}})\rangle -\displaystyle \frac{\partial }{\partial {k}_{i}}\left\langle {u}_{{n}_{1}}({\bf{k}})\right|\displaystyle \frac{\partial }{\partial {k}_{j}}{u}_{{n}_{2}}({\bf{k}})\rangle \right),\end{eqnarray} \tag{ B9a }$$

$$\begin{eqnarray}\begin{array}{r}\begin{array}{c}\begin{array}{ccl}{\left({Q}_{{n}_{1}{n}_{3}{n}_{2}}\right)}_{ij}({\bf{k}}) & = & {\left({{\boldsymbol{ \mathcal A }}}_{{n}_{1}{n}_{3}}\right)}_{i}({\bf{k}}){\left({{\boldsymbol{ \mathcal A }}}_{{n}_{3}{n}_{2}}\right)}_{j}({\bf{k}})-{\left({{\boldsymbol{ \mathcal A }}}_{{n}_{1}{n}_{3}}\right)}_{j}({\bf{k}}){\left({{\boldsymbol{ \mathcal A }}}_{{n}_{3}{n}_{2}}\right)}_{i}({\bf{k}})\\ & = & -\left(\left. \langle {u}_{{n}_{1}}\left({\bf{k}}\right)\right|\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}{k}_{i}}{u}_{{n}_{3}}({\bf{k}})\left.\left. \rangle \langle {u}_{{n}_{3}}({\bf{k}}\right)\right|\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}{k}_{j}}{u}_{{n}_{2}}({\bf{k}}) \rangle -\left. \langle {u}_{{n}_{1}}\left({\bf{k}}\right)\right|\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}{k}_{j}}{u}_{{n}_{3}}({\bf{k}})\left.\left. \rangle \langle {u}_{{n}_{3}}({\bf{k}}\right)\right|\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}{k}_{i}}{u}_{{n}_{2}}({\bf{k}}) \rangle \right).\end{array}\end{array}\end{array}\end{eqnarray} \tag{ B9b }$$

These can be recast into a more familiar form found after (6) of the main text by defining their vector components as ${\left({{\boldsymbol{\Omega }}}_{{n}_{1}{n}_{2}}\right)}_{i}\left({\bf{k}}\right)=\tfrac{1}{2}{\epsilon }_{{ijk}}{\left({{\rm{\Omega }}}_{{n}_{1}{n}_{2}}\right)}_{{jk}}\left({\bf{k}}\right)$ and ${\left({{\bf{Q}}}_{{n}_{1}{n}_{3}{n}_{2}}\right)}_{i}({\bf{k}})=\tfrac{1}{2}{\epsilon }_{{ijl}}{\left({Q}_{{n}_{1}{n}_{3}{n}_{2}}\right)}_{{jl}}({\bf{k}})$. In vector form, these become

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{\Omega }}}_{{n}_{1}{n}_{2}}({\bf{k}}) & = & {\rm{\nabla }}\times {{\boldsymbol{ \mathcal A }}}_{{n}_{1}{n}_{2}}({\bf{k}}),\\ {{\bf{Q}}}_{{n}_{1}{n}_{3}{n}_{2}}({\bf{k}}) & = & {{\boldsymbol{ \mathcal A }}}_{{n}_{1}{n}_{3}}({\bf{k}})\times {{\boldsymbol{ \mathcal A }}}_{{n}_{3}{n}_{2}}({\bf{k}}).\end{array}\end{eqnarray} \tag{ B10 }$$

The last terms in (5) and (6) of the main text involve the time derivatives of the projected operators. Using the above results, it is found that

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{{ \mathcal P }}_{D}({\bf{x}}) & = & {\rm{i}}\displaystyle \frac{e}{{\hslash }}\displaystyle \sum _{j}{\left({\dot{{\bf{A}}}}_{c}\right)}_{j}\left[{{ \mathcal P }}_{D}({{\bf{x}}}_{j}),{{ \mathcal P }}_{D}({\bf{x}})\right],\\ \displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}{{ \mathcal P }}_{D}({\bf{q}}) & = & 0.\end{array}\end{eqnarray} \tag{ B11 }$$

The results in (B4) and (B11) then give the time dynamics in (5) and (6) of the main text for the position x and gauge invariant wave vector k upon taking expectation values using a wave packet distribution $\rho ({\bf{k}})=\delta ({\bf{k}}-{{\bf{k}}}_{c})$.

We finally alert the reader that the interband properties in the main text were written as ${\bar{{\boldsymbol{ \mathcal A }}}}_{{n}_{1}{n}_{2}},{\bar{{\boldsymbol{\Omega }}}}_{{n}_{1}{n}_{2}},{\bar{{\bf{M}}}}_{{n}_{1}{n}_{2}},{\bar{{\bf{Q}}}}_{{n}_{1}{n}_{2}{n}_{3}}$ to denote that these were evaluated at k_c, while in this document these interband properties are functions of the general k.