Dimension transcendence and anomalous charge transport in magnets with moving multiple-$Q$ spin textures

Multiple-$Q$ spin textures, such as magnetic bubble and skyrmion lattices, can be driven into motion by external stimuli. The motion of spin textures affects the electronic states. Here we show that to describe correctly the electronic dynamics, the momentum space needs to be transcended to higher dimensions by including the ancillary dimensions associated with phason modes of the translational motion of the spin textures. The electronic states have non-trivial topology characterized by the first and second Chern numbers in the high dimensional hybrid momentum space. This gives rise to an anomalous electric charge transport due to the motion of spin textures. By deforming the spin textures, a nonlinear response associated with the second Chern number can be induced. The charge transport is derived from the semi-classical equation of motion of electrons that depends on the Berry curvature in the hybrid momentum space. Our results suggest that the motion of multiple-$Q$ spin textures has significant effects on the electronic dynamics and provides a new platform to explore high dimensional topological physics.

The static spin texture has significant effects on electronic wave functions and affects the dynamics of conduction electrons [21,22]. For instance, when electrons interact with a static skyrmion, they experience an effective magnetic field produced by the noncoplanar spin texture. As a consequence, there is a topological Hall effect, which has been observed experimentally [23][24][25]. The noncoplanar spin texture can open an energy gap in the electronic spectrum and stabilize a Chern insulator [22,[26][27][28]. The spin texture can be driven into motion by magnetic or electric field gradients [9][10][11], thermal gradients [12][13][14], and spin-polarized currents [29][30][31]. How the motion of spin textures affects the electronic dynamics is an important question. In ferromagnetic metals, the interaction between conduction electrons and moving spin textures can generate an emergent electric field that results in a spin motive force affecting the electronic dynamics [7,[32][33][34]. Nevertheless, how it works for a gapped system is not well understood.
In this paper, we reveal that the description of the dynam-ics of conduction electrons coupled to the moving multiple-Q spin textures requires a transcendent high dimensional momentum space. The high dimensional momentum space is spanned by the physical dimensions and ancillary dimensions associated with phason modes of the translational motion of spin textures. The electronic states can have nontrivial topology characterized by the first and second Chern numbers in the hybrid momentum space when the energy spectrum is gapped. Due to the nontrivial topology, there is an anomalous charge transport in magnets with moving multiple-Q spin textures. To explain the anomalous charge transport, we develop a general semi-classical theory to study the response of electronic states to the motion and deformation of spin textures. Our results demonstrate that the motion of multiple-Q spin textures has profound consequences on the electronic wave functions and dynamics.
spin texture, r i → r i +vτ, is depicted by the phase φ ν = Q ν ·vτ, where v is velocity and τ is time.
The coupling of conduction electrons to multiple-Q spin textures yields a magnetic superlattice. We denote the crystal momentum of the magnetic superlattice as k. The Hamiltonian in the momentum space is also a periodic function of φ ν besides k, i.e. H(k, φ) with φ = (φ 1 , φ 2 , · · · ). Therefore, when the spin texture is moving, the electronic dynamics is depicted in the effective momentum space spanned bỹ k = k ⊕ φ= (k x , k y , φ 1 , φ 2 , · · · ), a hybridization of crystal momenta and phason modes of the translational motion of spin textures. The motion of the spin textures is slow compared to electronic dynamics. For instance, for a spin texture with period 1 nm moving at a velocity 1 m/s, φ ν changes at a rate of 1 GHz. As a good approximation, the electronic spectrum E n (k) and the corresponding eigenstate |ψ nk evolve adiabatically with φ ν .

III. NONLINEAR RESPONSE
The semi-classical equation of motion for an electronic wave packet with the dispersion E n (k) and center of mass r n followsṙ µ n = where Ω µν n = i∂˜k µ u nk |∂˜k ν u nk − i∂˜k ν u nk |∂˜k µ u nk is the Berry curvature and |u nk is the periodic part of the Bloch state |ψ nk . Eq. (3) is meaningful only when µ = x or y since the center of mass r n is confined in the real space. Here the first term is the group velocity and the last term is the anomalous velocity due to nonzero Berry curvatures [35]. In the absence of external electromagnetic fields,k ν =k ν = 0 for ν = {x, y}. However, the translational motion of spin textures described by φ ν = ω ν τ ensures a nonzerok ν =φ ν = ω ν for ν = {1, 2, · · · }. In addition to the conventional anomalous velocity −k ν Ω µν n (k) [35], there is a new contribution −φ ν Ω µν n (k) that can be obtained in the semi-classical theory (see Appendix A-C). The Berry curvature in the new contribution is defined in the hybrid momentum space consist of the crystal momentum k µ and phason mode φ ν . The semi-classical equation of motion is derived in the leading order. The higher order corrections do not change the form of equation of motion but modify the dispersion relation and Berry curvature due to interband coupling (see Appendix B). Since the higher order corrections do not change our finally results (proved in Appendix D), the modifications of dispersion relation and Berry curvature are not considered here. The corresponding current density due to the translational motion of spin textures is where f (E n − E F ) is the Fermi distribution function and E F is the Fermi energy. When the energy spectrum is gapped in the hybrid momentum space and the Fermi energy E F is in the gap, the first term in Eq. (4) vanishes. The transported charge over one period ∆τ = 2π/|ω ν | at zero temperature is where C µν 1 is the first Chern number defined on thek µkν plane and is independent of the other momenta if the energy gap retains open in the entire hybrid momentum space, and N µ is the number of unit cells in the cross section perpendicular to the µ direction. Therefore the translational motion of the spin textures can result in quantized charge transport in magnetic insulators with C µν 1 0, which is the topological pumping. The response function Eq. (4) becomes nonlinear when the dynamics of φ ν couples withṙ n . This can be achieved by deforming the spin texture, Q ν → Q ν + Q ν . Therefore, φ ν = ω ν τ + Q ν · r n , and its dynamics followṡ that depends on the velocity of wave packet (see Appendix A). Substituting Eq. (7) into Eq. (3) and retaining the terms up to the second order in ω ν and Q νµ , we obtaiṅ where the first two terms correspond to the linear response in Eq. (4). The last term in Eq. (8) originated from the deformation of spin textures induces a nonlinear response where F µνγδ n = Ω µν n Ω γδ n + Ω µγ n Ω δν n + Ω µδ n Ω νγ n (see Appendix D). In this case, the transported charge over ∆τ at zero temperature is where L µ is the size of the cross section perpendicular to the µ direction and is the second Chern number [36] defined on thek µkνkγkδ hypersurface of the hybrid momentum space. Here we make the approximation that F µνγδ n ≈ dφ δ F µνγδ n /2π. The approximation is valid when the Berry curvatures are weakly dispersive  (a (a (a (a (a (a (a (a (a (a) ) ) ) ) ) (a) (a) (a) (a) (c (c (c (c (c (c (c) ) ) ) ) ) ) ) (c) (c) (c) (c) (c) (c) (c)  (d (d) ) ) ) ) ) ) ) ) ) ) ) ) ) (d) (e (e (e (e (e (e (e (e (e (e (e (e (e) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) as a function of φ δ . This is true when the spin texture is incommensurate with the lattice or the period of spin texture is much larger than the lattice constant. In our definition, both C µν 1 and C µνγδ 2 are antisymmetric under the permutation of the indices because Ω µν n = −Ω νµ n . Equation (10) demonstrates that, in addition to the contribution from the first Chern number, the charge transport has the contribution from the second Chern number when the spin texture is deformed.
To substantiate the anomalous charge transport due to the translational motion of the spin textures from the semiclassical analysis above, we study the charge transport in magnets with several typical multiple-Q spin textures obtained from the Monte Carlo simulation (see Appendix E). Here we consider three kinds of spin textures: double-Q and triple-Q collinear magnetic bubble lattices where spins are nearly parallel or antiparallel due to the strong easy axis anisotropy, and triple-Q skyrmion lattice as shown in Figs. 1(a)-1(c), respectively. For conduction electrons hopping on a square lattice (whose lattice constant is set to unity), the Q ν vectors for the double-Q magnetic bubble lattice are Q 1 = (2π/5, 0) and Q 2 = (0, 2π/5), while for the triple-Q magnetic bubble and skyrmion lattices are Q 1 = (−2π/5, 2 3π/5), and Q 3 = (4π/5, 0). For the triple-Q states, because 3 ν=1 φ ν = 0 for a translational motion, only two of φ ν are independent of each other.

IV. TOPOLOGICAL PROPERTY
To identify the topology of electronic states in magnets with the multiple-Q spin textures, we first study their electronic spectra in the hybrid momentum space. For the double-Q magnetic bubble lattice, the translational motion of the spin texture is parameterized by φ ν = ω ν τ. The energy spectrum as a function of φ 2 (with φ 1 = 0) is shown in Figs. 1(d) for the double-Q spin texture with J = 1.5t and B = 0. Here we use the periodic boundary condition along the x direction and the open boundary condition along the y direction. Apparently, there are topological edge states in the bulk energy gaps indicating the system is topologically nontrivial. For the triple-Q spin textures, to make them commensurate with the lattice, we compress the lattice constant along the y direction to 1/ √ 3. In this case, because 3 ν=1 φ ν = 0, we definẽ φ 1 = −(φ 1 + φ 2 )/2 = φ 3 /2 andφ 2 = (φ 1 − φ 2 )/2 that are independent. The change ofφ 1 andφ 2 corresponds to the shift of triple-Q spin textures along the x and y directions, respectively. Thus we can parameterize the translational motion of the triple-Q spin textures byφ ν = ω ν τ. The energy spectra as a function ofφ 2 (withφ 1 = 0) for the triple-Q magnetic bubble and skyrmion lattices with J = 1.5t and B = −1.2t are displayed in Figs. 1(e) and 1(f), respectively. The topological edge states and bulk energy gaps persist for the triple-Q spin textures.
To characterize the nontrivial band topology, we calculate the Chern numbers of the system with an efficient algorithm by using the U(1) link variable [37,38]. For commensurate spin textures, the system has translational symmetry and we can describe it in a high dimensional hybrid momentum space spanned by the generalized crystal momentumk = (k x , k y , φ 1 , φ 2 ) for the double-Q spin texture and k = (k x , k y ,φ 1 ,φ 2 ) for the triple-Q spin textures. Even if the spin texture is incommensurate with the lattice, k x and k y can still be introduced by using the twisted boundary condition [39]. Consequently, the Chern numbers in Eqs. (6) and (11) are well defined on the compact manifold. As an example, we consider the Fermi energy in the bottom bulk energy gaps in Figs. 1(d)-1(f). In the 4D hybrid momentum space, there are six first Chern numbers: C x1 1 = C y2 1 = 2 and C xy 1 = C x2 1 = C y1 1 = C 12 1 = 0, and one second Chern number: C x1y2 2 = −2 for the double magnetic bubble lattice. For the triple-Q magnetic bubble and skyrmion lattices, they have the same Chern numbers: C x1 Interestingly, our results demonstrate an insulator that is trivial in the crystal momentum space, can be topologically nontrivial in the hybrid momentum space.

V. ANOMALOUS CHARGE TRANSPORT
The translational motion of spin textures can induce a quantized charge transport according to Eqs.
time (ħ/t) (g (g (g (g (g (g (g (g (g (g) ) ) ) ) (g) (g) (g) (g) (g) magnetic bubble lattice, the currents along the x and y directions over ten periods are shown in Fig. 2(a). By integrating the currents over time, the average transported charges in one period are q x = 0.01e and q y = 15.95e. According to Eq. (5), q x = eN x ω 2 C x2 /|ω 2 | = 0 and q y = eN y ω 2 C y2 /|ω 2 | = 16e. Apparently, the number obtained by simulation shows a good agreement with that obtained from the semi-classical equation. For the triple-Q magnetic bubble lattice, we have q x = 0.14e and q y = 16.04e from the currents in Fig. 2(b), while for the triple-Q skyrmion lattice, we have q x = 0.00e and q y = 15.96e from the currents in Fig. 2(c). The numerical results for the triple-Q spin textures are also consistent with q x = 0 and q y = 16e obtained from Eq. (5).
The nonlinear response can be realized by deforming the spin texture as shown in Eqs. (7) and (9). For the double-Q magnetic bubble lattice, we deform the spin texture by setting φ 1 = Q 1 · r with Q 1 = (0, 2π/20). This induces a shear strain γ xy = −Q 1y /Q 1x = −0.25 on the spin texture that results in x → x + γ xy y, as shown in Fig. 2(g). The currents due to the translational motion and shear strain of the spin texture are displayed in Fig. 2(d). In this case, the transported charges in one period are q x = −3.92e and q y = 15.89e from the currents in Fig. 2(d). Meanwhile, according to the Eq. (10), we have q x = eL x C x21y = −2, and q y = 16e as before. The results from numerical simulation and semi-classical analysis are consistent with each other. Therefore, we demonstrate that, besides the longitudinal charge transport due to the linear response, there is a transverse charge transport caused by the nonlin-ear response of the deformed spin texture. Moreover, one can probe the Chern numbers experimentally from the quantized charge transport. For the triple-Q spin textures, the shear strain can also be induced by takingφ 1 = Q 1 · r, where we set Q 1 = (0, 2 √ 3π/40) and γ xy = Q 1y /Q 1x = −0.125 √ 3 in Figs. 2(h) and 2(i). The currents of the triple-Q magnetic bubble lattice are shown in Fig. 2(e) and the transported charges in one periodic are q x = −3.43e and q y = 15.60e. The currents of the triple-Q skyrmion lattice are shown in Fig. 2(f) and the transported charges are q x = −3.48e and q y = 15.92e. Equation (10) yields q x = −4e and q y = 16e for the triple-Q spin textures. In this case, the numerical results are less quantized. One possible reason is due to the fact that the spin textures from Monte Carlo simulation are not exactly periodic. As a comparison, we also calculate the current for the double-Q and triple-Q spin texture ansatz described by the linear superposition of cosine functions, and the transported charge for the ansatz are perfectly quantized (see Appendix G).
When the deformed spin texture is commensurate with the lattice, one may wonder whether the anomalous charge transport can be fully characterized by the first Chern numbers of the deformed system. To answer this question, we take the deformed double-Q magnetic bubble lattice studied above as an example. Under the deformation induced by Q 1 = (0, 2π/20), the new unit cell is enlarged to contain 5 × 20 lattice sites (fourth larger than the undeformed unit cell that is 5 × 5). Therefore on the 40 × 40 lattice, N x = 2 and N y = 8. Because the lattice constant along the y direction is increased by 4 times, the period of motion of the spin texture is also in-creased by 4 times (since ω 2 is fixed). By performing exactly the same calculation, we find the relevant first Chern numbers of the deformed system become C x2 1 = −8 and C y2 1 = 8. Plugging these into Eq. (5), we get q x = −16e and q y = 64e in one period of motion. To compare with the previous results, we need to divide the new q x and q y by 4 to count the transported charges in the same period of time as the undeformed case. Thus, the two independent approaches yield the same result. As shown in this example, even a small deformation of spin textures can increase the size of unit cell significantly. With the periodic modulation extended to a larger length scale, there are more energy bands due to the band folding, accompanying complicated band gap closing and opening in this process, where the topological phase transitions are possible. Therefore, without the nonlinear response function associated with the second Chern number, one needs to calculate different first Chern numbers for different deformation to understand the charge transport. This calculation turns out to be impractical for small deformation that results in huge unit cells. In this sense, the nonlinear response theory can make life much easier since one only needs to know the first and second Chern numbers of the undeformed system and these Chern numbers are usually relatively easy to calculate. VI. LINEAR PIEZOELECTRICITY.
As we have shown above, the deformation of spin textures is essential to generate the nonlinear response. The deformation of magnetic skyrmion lattices has been realized in a strained crystal [17] or by applying an electric field [40]. Very recently, the deformation of the moving magnetic skyrmion lattice under electric current flow has also been observed [41]. According to Eq. (9), the current depends linearly on the deformed Q v vector, i.e. j µ ∝ Q γδ . For instance, when the spin texture is deformed by an effective shear strain γ xy ∝ Q 1y in our numerical simulation above, the linear piezoelectric current is j x ∝ γ xy . To verify the linear piezoelectricity, we calculate q x as a function of γ xy and we expect q x ∝ γ xy . The numerical results are displayed in Fig. 3. For the finite system size (40 × 40 lattice sites) considered in our numerical calculations, the shear strain are limited to a few discrete values in order to retain the periodic boundary condition. For much larger systems, the shear strain can be tuned continuously. We show q x as a function of γ xy for the double-Q magnetic bubble lattice in Fig. 3(a), and for the triple-Q magnetic bubble and skyrmion lattices in Fig. 3(b). Here the symbols are from the numerical calculation and the black lines are from the Eq. (10). Apparently, the q x for the double-Q spin texture collapses on the black line and shows a good linear dependence on γ xy . The q x for the triple-Q spin textures are also linear in γ xy but deviate from the black line due to poor quantization of charge transport as reasoned above. Therefore, the linear piezoelectricity is verified numerically and works for moderate strains that preserve the intrinsic topology. Our results suggest a new mechanism of piezoelectricity as a consequence of the nontrivial topology. .

VII. DISCUSSION
The nonlinear response due to the translational motion and deformation of the spin textures is associated with the second Chern number. The quantized charge transport requires the Fermi energy in the bulk energy gap. When this is not the case, Eq. (9) is still valid but the transported charge is not quantized. Therefore, the anomalous charge transport generated by the motion of spin textures occurs both in insulating and metallic magnets. In our system, the electronic dynamics is governed by the nontrivial topology of electronic structure in the hybrid momentum space spanned by the crystal momentum and phason mode of the spin texture. This is different from the emergent electrodynamics induced by the motion of non-collinear spin textures according to the generalized Faradays law [7,32,33] or magnetization dynamics as a reciprocal process of the spin transfer torque [34]. Evidently, the anomalous charge transport in the collinear spin textures and in magnetic insulators can not be captured by the emergent electrodynamics. Our predictions can be tested in B20 chiral magnets [1,2], f -electron materials [42,43] and magnetic heterostructures [44][45][46], where multiple-Q spin textures have been observed in experiments. The quantized charge transport can be possibly realized in a heterostructure with the multiple-Q spin texture formed in a magnetic insulator layer or at the interface through the proximity effect [47][48][49]. In materials, the multiple-Q spin textures are pinned by defects. To drive the spin texture into motion, a large driving force is needed in order to overcome the pinning energy. The pinning-depinning transition provides a clear way to distinguish the anomalous charge transport from other contributions.
In summary, we reveal that the motion of spin textures has significant effects on the dynamics of conduction electrons. The semi-classical equation of motion of the electronic wave packet is described in a high dimensional hybrid momentum space. The hybrid momentum space is constituted of the physical dimensions and ancillary dimensions associated with the phason modes of the translational motion of the spin textures. The electronic states can be topologically nontrivial when the energy spectrum is gapped, and the nontrivial topology is characterized by the first and second Chern numbers in the hybrid momentum space. The nontrivial topology results in an anomalous charge transport. The charge transport is quantized and can be used to extract the Chern numbers experimentally. Our theory presents an entirely new mechanism for electric current generation induced by the motion of multiple-Q spin texture in magnetic insulators. Our results can be extended to 3D magnets with multiple-Q spin textures where even higher dimensional topological physics can be realized.

ACKNOWLEDGMENTS
We would like to thank Joel Moore, Chih-Chun Chien, and Cristian Batista for insightful discussions. Computer resources for numerical calculations were supported by the Institutional Computing Program at LANL. This work was carried out under the auspices of the U.S. DOE NNSA under contract No. 89233218CNA000001 through the LDRD Program, and was supported by the Center for Nonlinear Studies at LANL (Y. S. and S.Z.L). S.H. is supported by the JSPS KAKENHI Grants No. JP18K13488.

Appendix A: Semi-classical equation of motion
Here we consider a very general model of free electrons coupled to a periodic modulation potential described by and we study how the electron dynamics is affected by the motion and deformation of the potential. The potential is the linear superposition of a set of periodic functions V ν (Q ν · r) with the modulating vector Q ν . The perturbation can be induced by asking Here the first term shifts the potential and the second term deforms the potential. Apparently, the multiple-Q spin texture model belongs to the category depicted by Eq. (A1). In the semi-classical approach, the electronic state is described by a wave packet that is highly localized at the center of mass r c and k c in the phase space. Due to the localized nature, when the perturbation is smooth compared to the length scale spread by the wave packet, the approximate Hamiltonian for the wave packet can be obtained by linearizing the perturbation as φ ν (r, τ) ≈ φ ν (r c , τ) [35]. The leading order local Hamiltonian is and the first order correction is H 1 = ∂ r c H c · (r − r c ). Here we haveφ that depends on the velocity of the wave packet. Because H c has the same translation symmetry as H 0 , its eigenstates are still the Bloch states of H 0 as |ψ nk (r) = e ik µr µ |u nk (r) (A5) but with a shifted positionr µ = r µ + φ ν /Q νµ due to φ ν (r c , τ).
The wave packet to the leading order can be constructed from the Bloch states as where the electron dynamics is confined in a single band whose dispersion is E n (k) with n being the band index. The interband coupling is encoded in the higher order corrections and will be studied in the next section. In the following, we denote r n and k n as the center of mass of the wave packet |W n in the phase space. To be self-consistent, we have that indicates w n (k, τ) is highly localized at k n , and where we have used the integration by parts and the fact that ψ nk (r)| e ik γr γ |∂ k µ u nk (r) Here N is the number of unit cells (UC), s 0 is the area of a UC, and we parameterizer = R + ξ where R is the lattice vector and ξ is within the UC centered at the the origin. Therefore dr = R UC dξ since u nk (R + ξ) = u nk (ξ).
The equation of motion can be obtained from the La-grangian L = W n | i ∂ τ − H |W n . For the first term, we have where we use the fact that − ∂ τ arg[w n (k n , τ)] = − d dt arg[w n (k n , τ)] + k n µ ∂ k n µ arg[w n (k n , τ)] and the total time derivative is dropped because it just contributes a constant to the action and does not affect the equation of motion.
To the leading order, the second term of the Lagrangian is (A11) Now we combine Eq. (A8) and Eq. (A10) together that yields the final expression of the Lagrangian According to the Euler-Lagrange equation, the leading order semi-classical equation of the motion iṡ where the Berry curvatures Ω µγ n,kk = i ∂ k n µ u nk n |∂ k n γ u nk n − ∂ k n γ u nk n |∂ k n µ u nk n , (A14) Ω µν n,kφ = i ∂ k n µ u nk n |∂ φ ν u nk n − ∂ φ ν u nk n |∂ k n µ u nk n , (A15) are respectively defined in the crystal momentum space and the hybrid momentum space consist of the crystal momentum k n µ and phason mode φ ν . In this sense, φ ν does play the same role as the crystal momentum and the space dimension is transcendent accordingly. The dynamics of φ ν follows Eq. (A4) which is part of the equation of motion. The new contribution −φ ν Ω µν n,kφ to the velocity of wave packet is due to the translational motion of the periodically modulated potential.

Appendix B: Higher order corrections to the equation of motion
To include the higher order corrections to the equation of motion, we consider the wave packet with the first order corrections. According to the perturbation theory, the eigenstate of H c + H 1 up to the first order is where |u (1) nk (r) is the first order correction to the eigenstate and Thus the wave packet should be constructed from the corrected eigenstates as The self-consistent conditions of the wave packet are still k n µ = W n | k µ |W n and up to the first order. Here the first term is the same as Eq. (A8) and the second term is where we use the fact that u nk |u (1) nk = 0. Therefore the center of mass is r µ n = −∂ k n µ arg[w n (k n , τ)] − φ ν /Q νµ + i u nk n |∂ k n µ u nk n + i u nk n |∂ k n µ u (1) nk n + i u (1) nk n |∂ k n µ u nk n , and the first order corrections come from the interband coupling.
For the Lagrangian L = W n | i ∂ τ − H |W n , the first term now becomes where the leading order term W (0) n | i ∂ τ |W (0) n is the same as Eq. (A10) and the first order correction is The dispersion relation is corrected to the first order as Now the Lagrangian becomes L(r n ,ṙ n , k n ,k n ) = k n µ ∂ k n µ arg[w n (k n , τ)] − k n µφν /Q νµ + i φ ν u nk n |∂ φ ν u k n + i φ ν u nk n |∂ φ ν u (1) nk n + i φ ν u (1) nk n |∂ φ ν u nk n = k n µ −r µ n − φ ν /Q νµ + i u nk n |∂ k n µ u nk n + i u nk n |∂ k n µ u (1) nk n + i u (1) nk n |∂ k n µ u nk n − k n µφν /Q νµ + i φ ν u nk n |∂ φ ν u k n + i φ ν u nk n |∂ φ ν u (1) nk n + i φ ν u (1) nk n |∂ φ ν u nk n where the Berry curvatures with the higher order corrections areΩ µγ n,kk =Ω µγ n,kk + i ∂ k n µ u nk n |∂ k n γ u (1) nk n − ∂ k n γ u nk n |∂ k n µ u (1) Apparently, the equation of motion with the higher order corrections are still in the same form as the leading order equation of motion derived in the section above. The modified dispersion relation and Berry curvature with first order corrections are induced by the interband coupling.
Appendix C: Unify the dummy index In the notation above, the index ν = {1, 2, · · · } is for the ancillary dimensions associated with phason modes of the translational motion of the periodic modulation, while µ, γ = {x, y} are for the physical dimensions. To unify the dummy indices, we introduce the hybrid momentum space spanned by the generalized momentumk (C1) In the current notation, the equation of motion becomeṡ where the superscript ofk n is dropped hereafter and in the main text for simplification. Now µ, ν, γ = {x, y, 1, 2, · · · } are equivalent. However, Eq. (C2) is meaningful only for µ = x or y since the center of mass r n is confined in the real space.
Here the unified Berry curvature is Ω µν n = Ω µν n + Ω (1)µν n , Ω µν n = i ∂˜k µ u nk |∂˜k ν u nk − ∂˜k ν u nk |∂˜k µ u nk , where Ω µν n is the leading order Berry curvature and Ω (1)µν n is the first order correction to the Berry curvature. The dynamics of the generalized momentum followṡ (C4) The deformation of the periodic modulation potential coupleṡ k ν withṙ γ n that results in the nonlinear response of the system to the deformation. The nonlinear response function will be derived in the next section.
where D n (r,k) is the phase space density of states (PSDOS). In the absence of deformation, the phase space is homogeneous and the PSDOS is a constant as D n = 1/(2π) d where d is the dimension of the system. However, due to the deformation, the phase space is inhomogeneous and we need to consider the modification of PSDOS. The basic idea is to study the time evolution of phase space volume element ∆V = ∆r∆k that follows (1/∆V)d∆V/dt = ∇ r ·ṙ + ∇˜k ·k [50]. The equation can be solved [51] and yields the modified PSDOS D n (r,k) = 1 ∆V = 1 (2π) 2 where ε µνγδ and ε ηραβ are Levi-Civita tensor with ε xy12 = ε xy12 = 1. Plugging Eq. (C4) into Eq. (C2) and retaining the terms up to the second order in ω ν and Q νµ yieldṡ which is in the same form of Eq. (8) but with modified dispersion relation and Berry curvature due to the first order corrections. In the following, we will show that our finally results are not changed by the first order corrections which can be dropped. Now we substitute Eq. (D2) and Eq. (D3) into Eq. (D1) and retain the terms up the second order in ω ν and Q νµ that gives (D4) Here we collect different terms with similar orders and forms into four groups enclosed by brackets. The integral of the first group just gives Eq. (4) where higher order corrections are dropped. Since we focus on the Fermi energy in the bulk energy gap that remains open upon weak perturbation, the modified group velocity ∂ k µẼ n / and the first order correction to the Berry curvature Ω (1)µν n vanish under the integration over hybrid Brillouin zone [52]. Therefore, the transported charge described by Eq. (5) is unaffected by the higher order corrections. Notably, the second group in Eq. (D4) is unaffected by the higher order corrections. Because Ω µν n = −Ω νµ n , Q µν = −Q νµ , and the dummy indices can be freely replaced, we have for the second group.
Next, we are going to prove the integral of the last two groups are zero. For the third group, we use the same trick as above and we get (D6) For a fully occupied band at zero temperature, the integral of Eq. (D6) yields where we use the integration by parts and the fact that both E n andΩ µν n are periodic over the Brillouin zone. The three terms in the bracket equal to zero according to the Bianchi identity [53]. Therefore, the third group in Eq. (D4) does not contribute to the current density. The last group vanishes due to the antisymmetric property of Ω µν n and Q µν . To show the three terms are indeed canceled out, we take j y as an example and we collect all the terms associated with ∂E n ∂k y Q 1y Q 2x . Here the first two terms give ∂E n ∂k y Q 1y Q 2x Ω x1 n Ω y2 n − Ω y1 n Ω x2 n (D8) and the last term yields where we use the property Ω µν n Q γδ = Ω νµ n Q δγ , Ω µν n Ω γδ n = Ω νµ n Ω δγ n , and Q µν Q γδ = Q νµ Q δγ . Obviously, Eq. (D8) and Eq. (D9) are opposite to each other, and the last group has no contribution to the current density.

Appendix E: Monte Carlo simulation
We present details of the model and numerical simulations for the spin textures. We here consider two models in order to describe three multiple-Q states: one is the effective spin model derived from the Kondo lattice model for the double-Q and triple-Q collinear spin textures and the other is the frustrated spin model for the triple-Q skyrmion spin texture.
The effective spin model is given by [54,55] where S q is the Fourier transform of S i and Q ν [ν = 1, 2 (1, 2, 3) for the square (triangular) lattice] are the equivalent wave vectors for the multiple peaks in the bare susceptibility of itinerant electrons. The bilinear interactionJ = 1 represents the effective Ruderman-Kittel-Kasuya-Yosida interaction and the positive biquadratic interaction K = NK with N the number of sites represents the four-spin interactions leading to the multiple-Q instability. α is introduced as the anisotropy parameter for exchange coupling where α = 0 represents the Ising-type coupling. The third term represents the easy-axis anisotropy (A > 0) for localized spins. The simulations for the effective spin model are carried out with the classical Monte Carlo simulations at low temperatures. Our simulations are carried out with the standard Metropolis local updates. The results are obtained for systems with N = 96 × 96 sites under periodic boundary conditions. In each simulation, we perform simulated annealing to find the low-energy spin configuration by gradually reducing the temperature with the rate T n+1 = aT n , where T n is the temperature in the nth step. We set the initial temperature T 0 = 1.0 − 10.0. We take a = 0.99995-0.99999 and the final temperature is typically taken at T = 0.001-0.01. We also start the simulations from the multiple-Q spin configurations, such as the double-Q and triple-Q collinear spin textures. We determine the magnetic phase by comparing their energies of the obtained magnetic patterns in simulations. The double-Q collinear spin texture in Fig. 1(a) is obtained at T = 0.01, K = 0.5, α = 1.0, A = 0.5, and |Q| = 2π/6 on the square lattice, while the triple-Q collinear spin texture in Fig. 1(b) is obtained at T = 0.001, K = 1.0, α = 0.0, A = 0.0, and |Q| = 2π/6 on the triangular lattice.
The frustrated spin model on the triangular lattice is given by [56] where we consider the ferromagnetic nearest-neighbor coupling J 1 = −1 and antiferromagnetic third nearest-neighbor coupling J 3 = 0.5 in the first term, which gives the ordering vector |Q| = 2π/5. The second and third terms represent the Zeeman coupling to an external magnetic field and the easyaxis anisotropy, respectively. The optimal magnetic phases in the frustrated spin model are obtained from the Monte Carlo simulations based on the Metropolis algorithm. The lattices have N = 100 × 100 spins and the periodic boundary conditions. In the simulations, the 10 5 − 10 7 Monte Carlo sweeps measurements are performed after equilibration. The triple-Q skyrmion spin texture in Fig. 1(c) is obtained at T = 0.01, H = 0.27, and A = 0.5.

Appendix F: Time-dependent current
In the lattice model, the current through a bond connecting lattice sites i and j is where σ =↑ or ↓ for spin degree of freedom and |ψ n (τ) is the time evolution of the eigenstate |ψ n that is obtained by exact diagonalization of the Hamiltonian Eq. (1) at τ = 0. The time evolution of eigenstates follows the time-dependent Schrödinger equation i ∂ τ |ψ n (τ) = H(τ) |ψ n (τ) , which is solved numerically by the Runge-Kutta method. The total current I x (I y ) is the sum of all the bond currents through a cross section perpendicular to the x (y) direction.
Appendix G: Anomalous charge transport for the spin texture ansatz To compare with the numerical results for the spin textures from Monte Carlo simulation, we study the anomalous charge transport for the spin texture ansatz described by The spin texture ansatz is exactly periodic and hence well quantized anomalous charge transport is expected. Here we consider a double-Q collinear spin texture with Q 1 = (2π/a, 0) and Q 2 = (0, 2π/a), and a triple-Q collinear spin texture with Q 1 = (−2π/a, 2 3π/a), and Q 3 = (4π/a, 0). We first study the electronic spectrum to identify the topological property of the system. For conduction electrons hopping on a square lattice (whose lattice constant is set to unity), we fix a = 5 for both the double-Q and triple-Q spin textures depicted by Eq. (G1). While for the triple-Q spin texture, to make it commensurate with the lattice, we compress the lattice constant along the y direction to 1/ √ 3 such that each supercell also contains 5 × 5 lattice sites. The energy spectra as a function of φ 1 for the double-Q and triple-Q spin textures when J = 1.5t and B = 0 in Eq. (1) are shown in Figs. 4(a) and 4(b), respectively. Here we use the open boundary condition along the x direction and the periodic boundary condition along the y direction. Apparently, there are topological edge states in the bulk energy gaps indicating the system is topologically nontrivial. To characterize the nontrivial topology, we calculate the Chern numbers of the system. Here we focus on the Fermi energy in the bottom bulk energy gaps in Fig. 4. In the 4D hybrid momentum space of the double-Q spin texture, there are six first Chern numbers: C x1 1 = C y2 1 = 2 and C xy 1 = C x2 1 = C y1 1 = C 12 1 = 0, and one second Chern number: C x1y2 2 = −2. In the 5D hybrid momentum space of the triple-Q spin texture, there are ten first Chern numbers: C x1 1 = C x2 1 = −C y1 1 = C y2 1 = −2 and C xy 1 = C x3 1 = C y3 1 = C 12 1 = C 13 1 = C 23 1 = 0, and five second Chern numbers: We then study the anomalous charge transport induced by the translational motion of the spin textures, which is parameterized by φ ν = ω ν τ. Here we consider the spin texture moving in the y direction same as that in the main text. For the double-Q spin texture, we have ω 1 = 0 and ω 2 = 2πt/100 , while for the triple-Q spin texture, we have ω 1 = −ω 2 = 2πt/100 and ω 3 = 0. Thus the period of motion of both the double-Q and triple-Q spin textures is ∆τ = 100 /t. We calculate the currents due to the motion of spin textures on a 40 × 40 lattice (containing 8 × 8 supercells). For the double-Q spin texture, the currents over ten periods are shown in 5(a). By integrating the currents over time, the average transported charges in one period are quantized as q x = 0.00e and q y = 15.99e. According to Eq. (5), the transported charges are q x = eN x ω 2 C x2 1 /|ω 2 | = 0e and q y = eN y ω 2 C y2 1 /|ω 2 | = 16e. For the triple-Q spin texture, the currents over ten periods are shown in Fig. 5(c) that gives q x = 0.00e and q y = 31.97e. The transported charges are q x = eN x ω 1 C x1 1 /|ω 1 | + eN x ω 2 C x2 1 /|ω 2 | = 0e and q y = eN y ω 1 C y1 1 /|ω 1 | + eN y ω 2 C y2 1 /|ω 2 | = 32e according to Eq. (5). Apparently, the numerical results show a very good agreement with Eq. (5).
To trigger the nonlinear response, one needs to deform the spin texture. For the double-Q spin texture, we apply a shear deformation of the spin texture by deforming the Q 1 vector to Q 1 = (2π/5, 2π/20). By shifting the spin texture in the same way mentioned above, we obtain q x = −4.00e and q y = 16.00e from the currents in Fig. 5(b). The numerical results are also consistent with q x = eL x C x21y 2 ω 2 Q 1y /2π|ω 2 | = −4e and q y = 16e according to Eq. (10). For the triple-Q spin texture, we apply a shear deformation of the spin texture by deforming the Q ν vectors to Q 1 = (−2π/5, 2 3π/10). In this case, the currents over ten periods are shown in Fig. 5(d) that gives q x = −7.98e and q y = 31.96e. Eq. (5) yields q x = eL x C x12y 2 ω 1 Q 2y /2π|ω 1 | + eL x C x21y 2 ω 2 Q 1y /2π|ω 2 | = −8e and q y = 32e as before that agree with the numerical results. The numerical results of anomalous charge transport for the spin texture ansatz are perfectly quantized and exhibit very good agreement with the response function.