Charge Conservation Beyond Uniformity: Spatially Inhomogeneous Electromagnetic Response in Periodic Solids

Nonlinear electromagnetic response functions have reemerged as a crucial tool for studying quantum materials. Most attention has been paid to responses to spatially uniform electric fields, relevant to optical experiments in conventional materials. However, magnetic and magnetoelectric phenomena are naturally connected by responses to spatially varying electric fields due to Maxwell's equations. Furthermore, in moir\'{e} materials, characteristic lattice scales are much longer, allowing spatial variation of optical electric fields to potentially have a measurable effect in experiments. To address these issues, we develop a formalism for computing linear and nonlinear responses to spatially inhomogeneous electromagnetic fields. Starting with the continuity equation, we derive an expression for the current operator that is manifestly conserved and model independent. Crucially, our formalism makes no assumptions on the form of the microscopic Hamiltonian and so is applicable to model Hamiltonians derived from tight-binding or ab initio calculations. We then develop a diagrammatic Kubo formalism for computing the wavevector dependence of linear and nonlinear conductivities, using Ward identities to fix the value of the diamagnetic current order-by-order in the vector potential. We apply our formula to compute the magnitude of the Kerr effect at oblique incidence for a model of a moir\'{e} Chern insulator and demonstrate the experimental relevance of spatially inhomogeneous fields in these systems. We further show how our formalism allows us to compute the magnetic multipole moments and magnetic susceptibility in insulators. We next use our formalism to compute the second-order transverse response to spatially varying transverse electric fields in our moir\'{e} Chern insulator model, with an eye towards the next generation of experiments in these systems.

Hall Response when q = (q 0 , 0, 0 ) 16 The study of optical properties of materials has long been one of the guiding themes in condensed matter physics.From the ability of X-ray crystallography to determine crystal structure to spectroscopic probes of band structure and collective dynamics, optical response experiments have served as a vital tool to learn about the structure of solid state systems.On the theoretical side, sum rules for optical response functions provide some of the few exact and experimentally relevant results in quantum many body physics [1][2][3][4].Combined with first principles calculations, the Kubo formula for linear and nonlinear optical responses has allowed for an understanding of the electronic properties of materials [5][6][7][8][9][10].More recently, it has come to be appreciated that optical response functions such as the linear and nonlinear Hall conductivity [11][12][13][14], the photogalvanic effect [15][16][17], the magnetoelectric polarizability [18][19][20][21][22], and nonlinear magnetoresistance [23][24][25][26][27] all receive large topological contributions, and can be used to both diagnose and functionalize topological materials.
The fundamental objects of interest for optical responses are the (nonlinear) conductivity tensors of a material, which relate the measured current density ⟨j(r, t)⟩ at position r and time t that flow in a material in response to an electric field E(r, t).Typically, the electric field E(r, t) is introduced via a voltage source or an optical probe.The typical objects of interest are the experimentally-measurable (nonlinear) conductivities σ µ,ν1,...,νn (n) (r, r 1 , . . ., r n , t, t 1 , . . ., t n ), (1) which determine the measured current via ⟨j µ (r, t)⟩ = n i dr i dt i σ µ,ν1,...,νn (r, r 1 , . . ., r n , t, t 1 , . . ., t n )E ν1 (r Here and throughout we use Greek indices µ, ν = 1, 2, 3 to index spatial components of vectors, and the Einstein summation convention is used for repeated indices.Since we work entirely in flat space time, there is no functional distinction between upper and lower indices; we will try to choose the index arrangement that makes expressions easiest to parse.Eq. ( 2) parametrizes the current density order-by-order in powers of the external electric field.At order n = 1 we recover the familiar generalized Ohm's law with σ µν 1 (r, r 1 , t, t 1 ) = σ µν (r, r 1 , t − t 1 ) being the linear conductivity tensor.Note that in writing Eq. (2) we have not assumed that our system is translationallyinvariant.
Because the speed of light is large compared with typical velocity scales in a solid state material, and because the typical energy scales of interest are in the microwave, terahertz, and optical regimes, the typical wavelength of the probe electric field (on the order of thousands of angstroms for visible light all the way up to the order of meters for AC fields) is significantly larger than the typical lattice spacing of crystals (on the order of nanometers).This means that for most experimentally-relevant situations, we can ignore the spatial variation of the electric field [28].In this approximation, we can introduce the uniform optical conductivities σ µ,ν1,...,νn n (t, t 1 , . . ., t n ) which depend on time and not space, and which deter-mine the spatially averaged current density ⟨j(t)⟩ via ⟨j µ (t)⟩ = n i dt i σ µ,ν1,...,νn n (t, t 1 , . . ., t n ) i E νi (t i ).
The uniform optical conductivities can be computed using generalizations of the Kubo formula, in terms of correlation functions of the electronic velocity operator v. Recent theoretical work [29][30][31][32][33] has developed streamlined and efficient analytical machinery for computing σ µ,ν1,...,νn n (t, t 1 , . . ., t n ) for a variety of interesting materials, both for finite systems, continuum crystalline systems, and crystalline systems in the tight-binding approximation.For noninteracting electron systems, recent work has shown how the Kubo formula for the uniform optical conductivities can be interpreted in terms of the quantum geometry of electronic wavefunctions in the systems [34][35][36].This theoretical work has been complemented by advances in nonlinear spectroscopy, which have seen signatures of intriguing topological and bandgeometric effects in quantum materials [37,38].
That said, there are limitations to focusing on the uniform optical conductivity.First, the uniform optical conductivity captures response to spatially uniform electric fields, but does not capture response to magnetic fields.
To understand this, we can see from Faraday's law, that a spatially uniform and time-varying magnetic field generates a spatially inhomogeneous (transverse) electric field.Magnetic and magnetoelectric responses are thus encoded in the gradient expansion of σ µ,ν1,...,νn (n) (r, r 1 , . . ., r n , t, t 1 , . . .t n ).Such a gradient expansion has been carried out to low orders to obtain, for instance, the Streda formula relating the Hall conductivity to the magnetization density [39,40], the Kubo formula for the magnetoelectric polarizability [41], and optical gyrotropy [42].Similarly, a gradient expansion of the second order (n = 2) longitudinal optical conductivity yielded diagrammatic Kubo formulas for electric quadrupole transitions in periodic crystals [43].Lastly, the study of electron hydrodynamics has brought renewed attention to the connection between the viscosity tensor of an Galilean-invariant electron fluid and the (linear) optical conductivity expanded to second order in gradients [44][45][46][47][48][49][50].Attempts to generalize this relationship to periodic solids have yielded unintuitive results [51].
Second, the explosion of interest in moiré materials [52][53][54][55][56][57][58][59][60][61] gives cause to reevaluate the focus on the uniform optical conductivity.In systems such as twisted multilayer graphene, twisted transition metal dichalcogenides, and lattice-mismatched van der Waals interfaces, the effective moiré lattice constant can be several orders of magnitude larger than the atomic spacing in individual layers.In such a system, the wavelength of optical frequency light can be an appreciable fraction of the moiré lattice spacing.Thus, we can expect that optical properties such as Kerr and Faraday rotation, as well as nonlinear shift current and second harmonic generation may obtain measurable corrections from the spatial inhomogeneity of the applied optical electric field.Recent advances in ultrafast spectroscopy have put the study of nonlinear optical effects in moiré and van der Waals materials within reach of modern experiments [50,62,63].
Thus, in order to study both generalized magnetoelectric responses and optical properties of moiré materials, it would be desirable to have a complete theory of spatially inhomogeneous electromagnetic response applicable to periodic systems.From the Kubo formula, we know that σ µ,ν1,...,νn (n) (r, r 1 , . . ., r n , t, t 1 , . . ., t n ) can be computed from a retarded correlation function of n current operators j(r).When the full microscopic Hamiltonian H 0 for a system of interest is known, the operator j(r) can be identified via minimal coupling, by replacing the momentum of each particle p i by the covariant combination p i − eA(x i ), where A(x i ) is the electromagnetic vector potential as a function of the position operator x i of particle i, and e is the charge of particle i [64,65].Denoting by H A the Hamiltonian we obtain after minimal coupling, we can define the current operator as where δ/δA(r) denotes a variational derivative with respect to the vector potential.This allows us to write from which the nonlinear conductivities can be extracted via standard response theory upon identifying (in this gauge [64,65]) the electric field E = −∂A/∂t.Crucially, the current defined via Eq.( 5) manifestly obeys the continuity equation (charge conservation equation) where ρ(r) is the charge density operator.Eq. ( 7) reflects the conservation of electric charge and is equivalent to gauge invariance.While this gives an in-principle complete procedure for computing σ µ,ν1,...,νn (n) (r, r 1 , . . ., r n , t, t 1 , . . ., t n ), it has one main drawback: in most situations in solid state physics we do not have access to the full microscopic Hamiltonian, but instead only an effective low-energy model for the degrees of freedom of interest.This can arise due to an approximate treatment of relativistic effects of core electrons in a solid (as in, e.g., a DFT calculation), or via the use of a tight-binding model that describes only the low-energy subset of states in the full Hilbert space of a material.When only an effective or approximate Hamiltonian is known, the minimal coupling substitution H 0 → H A cannot be carried out uniquely.
One common approach [41,42,[66][67][68] to circumvent this problem in nonrelativistic materials is to assume that the full and unknown microscopic electronic Hamiltonian has the form where m is the electron mass, V is the periodic potential of the ions in the system, ⃗ λ is the approximate spin-orbit potential (usually expressed in terms of gradients of V ), ⃗ σ is a vector of Pauli matrices acting on electron spins, and U is the Coulomb potential.With this assumption, the minimally coupled current operator can be evaluated exactly.In the absence of the external electromagnetic field, it is given by j min (r) ≡ − δH A δA(r) A=0 = q e iq•r j q,min (9) where the Fourier components j q,min can be expressed in two equivalent and often used forms = e i e −iq•xi/2 v i e −iq•xi/2 = j mid,q , where v i = i [H, x i ] is the velocity operator for particle i, and e is the charge of the particles (which for convenience we assume to be the same for all particles throughout this work).The customary approach is to then take either Eq.(10) or Eq.(11) as the definition of the current operator in the low-energy effective model, and use that to compute optical response functions.
There are two problems with this approach.The first is that, although Eqs.(10) and (11) are equivalent for the microscopic Hamiltonian Eq. ( 8), they will in general give different results when applied to any effective model.This raises the question of which, if either, of them should be used to compute optical response functions.Additionally, neither Eq. (10) nor Eq. ( 11) coincide with the minimally coupled current j A (r) when relativistic corrections to the kinetic energy are incorporated, which may be important for core electrons in solid state systems.Even more severe, we will show below that generally neither Eq. (10) nor Eq. ( 11) satisfy the continuity equation Eq. ( 7) even when A = 0.When used in the context of effective low-energy models for even nonrelativistic systems, this means that Eqs.(10) and (11) represent nonnumber-conserving approximations to the minimally coupled current.Thus neither Eq. (10) nor Eq. ( 11) provide a satisfactory starting point for the computation of optical response functions in low-energy models of quantum materials.
In what follows, we will develop a theory of spatially inhomogeneous linear and nonlinear optical response that manifestly obeys the continuity equation Eq. ( 7) when applied to arbitrary effective models of condensed matter systems.First in Section II, we show how to construct the position-dependent current operator for a generic Hamiltonian that manifestly respects the continuity equation.Crucially, our construction makes no assumptions on the form of the Hamiltonian, and is applicable to nonrelativistic, semi-relativistic, and approximately tightbinding Hamiltonians.We show that our current agrees with Eqs.(10) and (11) to linear order in the wavevector q, implying that our theory gives the same linear response to uniform electric and magnetic fields as in Refs.[18,41,42,66].We will also show how to define the diamagnetic current such that the continuity equation is satisfied order-by-order in the electromagnetic vector potential.
Next in Section III, we set out the Feynman diagrammatic rules for evaluating the spatially inhomogeneous conductivity in Fourier space as a function of frequency and wavevector (where the uniform optical conductivities correspond to the limit of zero wavevector).Focusing on the linear conductivity, we derive a modified f-sum rule relating the diamagnetic conductivity to the density-density response function.In Section V, we explore various observable phenomena that arises from our wavevector-dependent response theory at the linear response level, applied to non-interacting electron systems in a periodic potential.As a proof of concept, we examine the wavevector-dependence of the Hall conductivity in a time-reversal breaking Weyl semimetal model first.We then show how our formulation of the wavevector dependent linear conductivity can be applied to compute the Kerr rotation and ellipticity in a model of a moiré-Chern insulator.We compare our predictions with an analogous calculation using the current operators in Eqs.(10) and (11), showing that our gauge-invariant formulation predicts a measurably different Kerr response on the order of arcseconds to arcminutes for typical moiré length and energy scales.We round off this section be examining the magnetic properties of insulators.We analyze the relationship between magnetic susceptibility and transverse conductivity, present expressions for the magnetic quadrupole moment, and give a new derivation of the Streda formula for the Hall conductivity.
Then in Section VI, we move on from phenomena derived from linear conductivity, and evaluate the diagrammatic Kubo formula for the second-order electromagnetic conductivity for noninteracting electrons in a periodic potential.We then focus on the component of the secondorder conductivity that is second harmonic generation in frequency, and self-focusing in wavevector, which is relevant to future ultrafast optical experiments.We show how to compute this component of the conductivity for our model of a moiré-Chern insulator.Finally, we conclude our paper by summarizing the key results and suggesting new directions in Section VII.

II. THE CURRENT OPERATOR
Motivated by recent efforts to obtain a complete spatially-inhomogeneous electromagnetic response theory [43,51,66,68], our goal in this section will be to derive the Fourier components of the conserved current operator in a model-independent fashion.Measurable phenomena such as conductivities, optical response, and magnetization can only be accurately modeled in terms of a charge-conserving current operator.As we will show, the conventionally defined current densities Eq. ( 11) and (10) fail to satisfy the continuity equation (7) for generic Hamiltonians.
To find the conserved current density, let us consider a general Hamiltonian for a (possibly interacting) system of electrons.Here, i = 1 . . .N indexes the particles of the system, p i is the momentum operator for particle i, and x i is the position operator for particle i.Also, T (p i ) is the kinetic energy operator for the i-th particle.Although this kinetic operator is usually quadratic in momenta, it can extend to semi-relativistic systems.For example, the kinetic energy could take the form which is relevant for all-electron ab initio modeling of spin-orbit coupled materials [69].Our results will hold for arbitrary T (p i ).The single particle potential V (x i , p i ) includes both the momentum-independent external potential as well as momentum-dependent external potential terms such as the spin-orbit potential.Finally, U (x i − x j ) is the interaction potential, which we take to be momentum-independent for simplicity (we expect this is a good assumption for most models of interest; however, for long-range interacting models of Mott insulators that have recently attracted attention, this assumption must be checked [70]).We assume that the potential V has discrete translation symmetry such that where R is a Bravais lattice vector.The global U (1) symmetry of H 0 implies that the density operator satisfies the continuity equation (7), where all particles have charge e.We can rewrite the continuity equation in the Heisenberg picture as where we work in natural units where ℏ = c = 1 unless stated otherwise.In the presence of an external electromagnetic field with scalar potential A 0 and vector potential A, we have that the minimally coupled Hamiltonian becomes We can we take the functional derivative with respect to the scalar potential to arrive at where ρ A (r) indicates the density operator in the presence of the electromagnetic field ρ(r) is defined in Eq. ( 15), and we have made use of the identity Furthermore, we have that where by evaluating the commutator we find Eq. ( 21) allows us to define the current in terms of the the Hamiltonian H A minimally coupled to the external electromagnetic field.In particular, we can find the unperturbed current operator j min (r) from Eq. ( 16) as However, in order to utilize Eq. ( 22), we need to know the explicit microscopic form of the Hamiltonian in Eq. (12).In many cases of interest, we do not have this information available.For instance, in modeling electronic systems we often only have access to Wannier-based tightbinding models, or other effective Hamiltonians, that are projected into a subset of bands of interest.For an effective Hamiltonian, it is not possible to implement the minimal coupling of Eq. ( 17) in a model independent fashion.One possible way to circumvent this complication is to make use of a Peierls substitution to couple a Wannierized effective model to the vector potential, provided the Wannier orbitals are sufficiently well localized and the vector potential varies sufficiently slowly [4].However, defining the current via a Peierls substitution requires making choices for how to evaluate path integrals of the vector potential between Wannier centers that are not necessarily natural.It also requires assuming that the vector potential varies slowly on the length scale of the Wannier function localization length, which may not be appropriate for all models (especially for models of moiré or topologically nontrivial systems).Finally, the Peierls substitution cannot be used to define the current when a Wannierized tight-binding model has not been explicitly constructed.
In order to circumvent these issues, we will derive an alternative expression for the current j(r) that is manifestly conserved and that does not require detailed knowledge of the microscopic Hamiltonian.To do so, let us first recast the continuity equation, Eq. ( 16), in momentum space.We introduce the Fourier-transformed density In Fourier space, the continuity equation reads where j q = dr e −iq•r j(r).
Rather than attempt to evaluate the commutator in Eq. ( 24) directly, we will take a more general approach.
In the Heisenberg picture, we see that ρ q in Eq. ( 23) depends on time implicitly through the operators x i .Defining the single-particle velocity operators we can write Thus, using the definition of the derivative, the rate of change of the density operator is To simplify this further, we make use of a general result of Karplus and Schwinger [71], who showed that for any operators A and B, We prove this result in App. A. Inserting Eq. ( 29) with A = −iq • x i and B = −iq • v i into Eq.( 28), we find that Comparing with the continuity equation in Eq. ( 24), we can identify Eq. ( 31) has several desirable features.Most importantly, it is manifestly conserved.Second, for nonrelativistic Hamiltonians (where v(p, x) is a linear function of p), a short calculation (see App. B) shows that j q coincides with the minimally coupled current.Note that since the Hamiltonian H 0 commutes with the projection operator onto a set of low-energy bands, the projected low-energy density operator satisfies Eq. ( 16) with the conserved current given by the low-energy projection of Eq. (31).Thus, for nonrelativistic systems with microscopic Hamiltonians of the form of Eq. ( 8), Eq. ( 31) is the approximation to the minimally coupled current in Eq. ( 5) that is manifestly conserved in any low-energy approximation to the Hamiltonian.However, we emphasize that a choice was made in going from Eq. (30) to Eq. ( 31): we can add any operator δj q that is purely transverse (q • δj = 0) to Eq. ( 31) without changing the continuity equation (16).Different choices of δj correspond to different definitions for the (orbital) magnetization current, which do not contribute to transport and are not constrained by gauge invariance.As such, it is in principle not possible to fix δj unambiguously from effective Hamiltonians alone.Even in the context of the microscopic Hamiltonian nonrelativistic approximations to the Dirac equation, nonminimal coupling to the electromagnetic field can modify the definition of the magnetization current.Faced with these challenges, we will take the natural choice δj = 0-corresponding to our definition Eq. ( 31) for the current operator-as our working, model-independent definition of the current.
As we will now show, Eq. ( 31) resolves many of the problems of the non-conserved choices Eqs.(10) and (11) for the current.In Sec.II A we will show that both the longitudinal and transverse components of Eqs.(10) and (11) differ from those of Eq. ( 31) at second order in wavevector.Then, in Sec.II B we will compute the matrix elements of the conserved current Eq.(31) in the basis of Bloch eigenstates, allowing us to apply our formalism to compute electromagnetic responses.Finally, in Sec.II C, we will show that the formula of Karplus and Schwinger can be used to define the generalized wavevector-dependent diamagnetic contributions to the current operator in the presence of a nonzero electromagnetic field.

A. The Non-Conserved Current Operators
Eq. ( 31) is our first main result.Let us compare it quantitatively with the (non-conserved) current operators Eqs. ( 10) and (11), which have been used to compute spatially inhomogeneous electromagnetic response functions to low order in the wavevector [51,66].Recall that Eq. (10) proposes the definition for the current operator.We will refer to this as the "trapezoidal" current, since it is the trapezoidal approximation the integration over λ in our conserved current Eq.(31).Similarly, Eq. ( 11) proposes the commonly-used [51,66,68] definition We will refer to this as the "midpoint" current, since it is the midpoint approximation the integration over λ in our conserved current Eq.(31).
To see how the definitions of current operator from Eqs. (32) and (33) compare with our Eq.(31), it is helpful to write v i = v i (p i , x i ) to make explicit the p i -and x idependence of v i .Inserting this into Eq.( 31) and using the fact that we find Similarly, using Eq. ( 34) to simplify Eq. ( 32) yields Lastly, using Eq. ( 34) to simplify Eq. (33) yields We thus generically have that j q ̸ = jq ̸ = j mid,q .Note that if v i is a linear function of momentum (which occurs when the kinetic energy includes no relativistic corrections and the spin-orbit potential is linear in momentum), then we can carry out the integral in Eq. ( 35) to find that j q = jq = j mid,q .However, for Hamiltonians with relativistically-corrected kinetic energy and complicated spin-orbit potentials, this will not be the case.This is relevant not just for ab-initio calculations of charge transport in heavy elements, but also for effective models where integrating out high energy degrees of freedom can lead to a renormalization of the effective Hamiltonian away from the standard nonrelativistic form.
We can compute the difference j q − jq and j q − j mid,q order-by-order in q.We find that We see from Eqs. (38) and (39) that j q , jq , and j mid,q all coincide to linear order in q, with the first discrepancy at order |q| 2 .This suggests that while calculations of uniform magnetoelectric response functions (which requires knowing j q to linear order in q, as in Refs.[41,42,66]) can use jq or j mid,q in place of the conserved Eq. ( 31), the calculations of the O(|q| 2 ) corrections to the conductivity in Ref. [51,68] should be reexamined.We will show below that while j q , jq , and j mid,q all yield the same Hall conductivity to order O(|q| 2 ), they yield different predictions for the longitudinal conductivity at the same order.Even further, we see from Eq. ( 38) that q • j q − jq ̸ = 0 and q • (j q − j mid,q ) ̸ = 0, so that only j q and not jq nor j mid,q satisfies the continuity equation Eq. ( 16).Generically, we have q • j q − q • j min,q = 0, (40) In App.B we give a direct analysis of how the conserved j q differs from non-conserved jq and j mid,q for a concrete model of a semi-relativistic free-electron system.We show explicitly that once the free electron is energetic enough for quartic or higher-order corrections to the kinetic energy to become appreciable, the conserved j q and non-conserved jq and j mid,q will disagree at order |q| 2 .
Although we have shown how the conservedj q of Eq. ( 31) differs from the nonconserved jq and j mid,q , in the current form our results only apply to continuum first-quantized Hamiltonians of the form of Eq. (12).To apply our same approach to formulate the current density operator for a wider variety of condensed matter problems, we will in Sec.II B calculate the second-quantized form of the current operator Eq. ( 31) by examining its matrix elements in a basis of single-particle Bloch eigenstates.

B. Matrix Elements of the Current Operator in
the Bloch Basis with Second-Quantization In this section, we will evaluate the matrix elements of the conserved current operator Eq. ( 31)-as well as the non-conserved current operators Eqs. ( 32) and ( 33)-in a basis of single-particle Bloch eigenstates.Consider a noninteracting Hamiltonian of the form of Eq. ( 12) with U = 0 with discrete translation symmetry as in Eq. ( 14).From Bloch's theorem we can introduce single-particle eigenstates where N is the nominally infinite number of unit cells in the system.We denote the eigenstates with kets |ψ nk ⟩ with inner product It will also be convenient to introduce kets |u nk ⟩ to represent the cell-periodic functions u nk (r).The inner product of cell-periodic kets is given by where cell denotes an integration over a single unit cell.We will work in the periodic gauge such that for any reciprocal lattice vector G.Note that the periodic gauge always exists, even for topologically-nontrivial systems, as we discuss further at the end of this subsection.The periodic gauge constraint implies [72] |u nk+G ⟩ = e −iG•x |u nk ⟩ .
Let us begin by computing the matrix elements of the density operator ρ q = e i e −iq•xi between two Bloch states.We have Introducing creation and annihilation operators c † nk and c nk that create or annihilate electrons in the state |ψ nk ⟩, Eq. ( 48) implies that we can write the density operator in second-quantized notation as Next, we can examine the second quantized form of our conserved current j q from Eq. (31).Note first that Eq. (50) indicates that j q |ψ nk ⟩ has crystal momentum k − q.Since Bloch states with different crystal momentum are orthogonal, we deduce that the only nonvanishing single-particle matrix elements of j q are where we have defined Using Eq. ( 51), we conclude that the current j q has the second-quantized representation By the same logic, we find the second-quantized representation of the non-conserved "trapezoid" current jq from Eq. ( 32) is Finally, for the non-conserved "midpoint" current j mid,q in Eq. ( 33) we find the second-quantized representation Eq. ( 53) is the main result of this section, and allows us to reformulate our conserved current operator Eq. ( 31) in terms of second-quantized Bloch orbitals.Eq. ( 53) can be used to compute the matrix elements of the conserved current for effective models that include only a subset of the energy bands of interest in a solid, in which the effective kinetic energy need not be a quadratic function of momentum.Additionally, Eq. ( 53) manifestly satisfies the continuity equation, and so gives a proper starting point for the evaluation of linear and nonlinear electromagnetic response coefficients.
Although we supposed that the wavefunctions |u nk ⟩ were the (cell periodic parts of the) energy eigenstates for the non-interacting part of a Hamiltonian, we can follow the logic leading to Eq. ( 53) for any orthonormal Bloch-like basis set.In App.C for example, we derive the second-quantized representation of the current Eq. ( 31) in a basis of ultra-localized tight-binding orbitals which will be useful for calculations in model systems.Additionally, provided that the interaction energy depends only on (unprojected) density operators, Eq. (53) will still give the second-quantized current operator.
Note that for systems with nonvanishing Chern numbers, the periodic gauge constraint for |u nk ⟩ in Eq. (46) requires that the phase of |u nk ⟩ be non-smooth in the Brillouin zone [73,74].This presents no difficulty for our formalism, which will not involve derivatives or Taylor expansions of the wavefunctions themselves with respect to k, but only derivatives of the manifestly-smooth operator H k .Additionally, when we work in a basis of tightbinding orbitals as in App.C, we can include a sufficiently large number of bands to guarantee that the basis Wannier orbitals have smooth Fourier transforms.Finally, we note that in any numerical computation, the wave functions |u⟩ nk are typically evaluated on a discrete grid of points in the Brillouin zone; we can then choose a periodic gauge such that any nonremovable singularities due to nonvanishing Chern number occur at points that are not included in the grid.
Having now expressed the current operator Eq. ( 31) in the second-quantized Bloch basis, we can now turn our attention to how the current density is modified in the presence of a background electromagnetic field.This will allow us to generalize the usual diamagnetic current to nonzero wavevector for systems with generic Hamiltonians.We will see that these generalized diamagnetic cur-rent operators-which are essential to maintain gaugeinvariance of response functions-appear naturally in our formalism for j q defined in Eq. (31).

C. Diamagnetic Current Operators
In order to compute electromagnetic response functions, we need to know the current j A (r) order-by-order in the vector potential A. Although we have primarily focused on j q,min in Eq. ( 22), for a general system, we can write or in Fourier space While knowledge of j µν q,−q ′ requires knowledge of the minimally coupled Hamiltonian in Eq. ( 17), the continuity equation at Eq. ( 16) places constraints on the longitudinal components of j µν q,−q ′ .In particular, note first that Inserting Eq. ( 58) into the continuity equation ( 16), we find that where we have used from Eq. ( 23) that ρ q is independent of the vector potential.Eq. (59) shows that the longitudinal component of the diamagnetic current j µν q,−q ′ can be expressed entirely in terms of the unperturbed current j ν −q ′ ,min .Furthermore, since the longitudinal component of j ν −q ′ ,min is equal to the current in Eq. ( 31), we find that the longitudinal component q µ q ′ ν j µν q,−q ′ of the diamagnetic current is entirely expressible in terms of the velocity operator.
Eq. ( 59) is equivalent to the Ward identity from quantum field theory [75].We can iterate the logic leading to Eq. ( 59) in order to obtain the longitudinal part of all higher-order current vertices q µ1 q ′ µ2 q ′′ µ3 . . .j µ1µ2µ3... q,q ′ ,q ′′ ... .Since the operators j µ1µ2µ3... q,q ′ ,q ′′ ... appear as vertex functions in calculations of higher-order electromagnetic responses [29,76], we see that our expression for the current in Eq. (31) in terms of the velocity operator can be used to generate the longitudinal components of these higher-order vertices.
Using Eqs. ( 58) and ( 59), we can recast this second order operator in terms of an integration over an auxiliary variable λ via the Karplus-Schwinger relation.Namely, we show in App.D, that for any single-particle operator Applying Eq. ( 60) to our Ward identity Eq. ( 59) we find where we have used the fact that both j min,q and j q of Eq. ( 31) are conserved.The advantage to using the Karplus-Schwinger relation is being able to strip off the q and q ′ , thus allowing us to define j µν q,−q ′ that satisfies the Ward identity for generic systems.
We can rewrite Eq. ( 61) in terms of the secondquantized creation and annihilation operators for Bloch eigenstates following the formalism of Sec.II B in order to define In a similar manner, we can iterate this procedure to derive expressions for the N -th order current operator which satisfies the generalized Ward identity Applying the Karplus-Schwinger relation Eq. ( 60) iteratively, and using our expressions Eqs.(31) and ( 53) we find that we can write the n-th order diamagnetic current operator in second-quantized notation as Note, importantly, that Eq. ( 65) is explicitly symmetric under the exchange of any pair of indices, since partial derivatives commute.
For later convenience, we will define the N -th order velocity vertex v ν N (N ) (k, q, −q 1 , • • • , −q N ) as the operator appearing in the matrix element of Eq. ( 65), i.e., Therefore, applying Eq. ( 66) to Eq. ( 65), we can simply write Note that since the "trapezoid" and "midpoint" currents [Eqs.(10) and (11)] do not satisfy the continuity equation, there is no general procedure for determining their corresponding diamagnetic currents.Nevertheless, for the sake of comparison, we will define "trapezoid" and "midpoint" diamagnetic current operators by approximating each auxiliary integral in Eq. (65) with the trapezoid or midpoint approximation, respectively.
Given Eq. ( 65) for the generalized diamagnetic current operators as a function of wavevector(s), we will now develop a diagrammatic formalism to compute linear and nonlinear conductivities as a function of frequency and wavevector.Our construction of Eq. ( 65) and Eq. ( 53) from generalized Ward identities ensures that these conductivities will give a current that respects the continuity equation.

III. FEYNMAN DIAGRAMMATIC METHOD AND SUM RULES
Using the conserved current operator, Eq. ( 31) will now develop the formalism for computing the spatiallyinhomogeneous linear and nonlinear conductivities defined in Eq. ( 2).We will build off the work of Ref. [29], which presented a simple yet powerful framework for calculating the spatially uniform (q = 0) nonlinear conductivities.Continuing off of this framework, we look to expand this framework to accommodate spa-tially inhomogeneous fields and currents.This framework utilizes diagrammatic perturbation theory in terms of non-interacting Matsubara Green's functions for electrons, with interactions with the external electromagnetic field given by the (noninteracting) vertex functions governed by Eq. (66).We will generalize this diagrammatic method of Ref. [29,33] to allow for electromagnetic fields and currents with nonzero wavevectors.As a consequence, we will be able to capture not just response to electric fields, but-via Faraday's law Eq.(4)-response to magnetic fields as well.We begin in Sec.III A by establishing the foundation for diagrammatic evaluation of the Kubo formula.Then, in Sec.III B we will give the rules for setting up and evaluating our Feynman diagrams.Finally, in Secs.IV-VI we will apply our formalism to analyze the linear and second order response as a function of wavevector for 3D Weyl and 2D moire systems.

A. Feynman Diagram Setup
To write the Feynman diagram rules for evaluating the average current, we must first relate the conductivities to the generating function for correlations.Since the Feynman diagrams are a direct interpretation of perturbative expansions of the generating function, this allows us to write the conductivities as a sum of diagrams, which may then be translated to a mathematical statement involving Green's functions and interaction vertices [29,33,[77][78][79][80][81][82][83].
First consider the generic Hamiltonian that has been coupled to a vector potential, as posed in Eq. (17).Notice that the minimally coupled Hamiltonian can be expanded in powers of the vector potential as: In the last term of Eq. ( 68), we replaced the functional derivatives of the coupled Hamiltonian with the generalized current term described in Eq. ( 63).We will now specify to studying systems of electrons, where we take e = −e, with e being the (positive) elementary charge.With this definition of the minimally coupled Hamiltonian, we can define the expectation values of the current operator to all orders in the vector potential [29] where j µ A (r) is defined via Eq.( 56) in terms of our conserved current, Eq. ( 53), and the diamagnetic currents, Eq. 65.Also, T t is the time-ordering operator.
In Eq. ( 69), two levels of expansion may be implemented in orders of the vector potential: an expansion in j µ A (r) and an expansion in H A .We will assume the vector potential can be written as a superposition of plane waves, To obtain the average current, we can expand the lefthand-side of Eq. ( 69) in terms of a Fourier-transformed response coefficient and the corresponding fields, as described in Eq. ( 2), taking the Fourier transform of every time and position argument.
Note that for systems of electrons in a periodic potential, Umklapp processes allow for the generation of currents that oscillate in space with wavevectors that differ from the wavevectors of the applied field by a reciprocal lattice vector.Concretely, expanding Eq. ( 69) in Fourier space [or, equivalently, taking the Fourier transform of Eq. ( 2)] yields (q, q 1 , . . . where G are the reciprocal lattice vectors.Since the spatial scale of variation of a current with wavevector G ̸ = 0 is on the order of (fractions of) an angstrom, these Umklapp terms are typically not experimentally interesting.Therefore, we will focus the remaining work on computing the G = 0 component of the electromagnetic response.We emphasize, however, that our diagrammatic calculation can be easily generalized to compute the G ̸ = 0 Umklapp components of the current as well.
By comparing Eq. ( 69) with Eq. ( 71), we may extract the n-th order conductivity by using the following process: we first expand the average current to n-th order in the vector potential.Next, we Fourier transform with respect to time, and make use of the definition of the electric field in the A 0 = 0 gauge.Note that this produces a timeordered, rather than a causal response function.To recover the usual causal response function that can be measured in experiment, we can analytically continue all frequencies into the complex plane, ω → ω + iη for a small infinitesimal η, following the procedure outlined in Refs.[29,65,80].In additional to preserving causality, a small positive finite η can be used to approximate the electron self-energy due to impurity scattering [29,65].
Practically speaking, we will carry out the perturbative calculation using the imaginary-time Matsubara formalism, using Matsubara frequencies iν n = (2n+1)iπ β for fermions and iω n = 2niπ β for bosons, where β = (k B T ) −1 [78].We will take the T → 0 limit at the end of all computations, which has the effect of turning Fermi-Dirac distributions, n F , into Heaviside step functions.Moreover, the T → 0 limit will still give a good approximation for low temperature transport coefficients in insulators.Since the Matsubara frequencies become continuous in this limit, we will suppress the subscript n [29,78,81,82].
With this as our starting point, we can now outline the diagrammatic rules that allow us to express Eq. ( 71) in terms of Matsubara Green's functions and the vertex functions of Eq. ( 65).

B. List of Feynman Diagram Rules
In this section, we will describe the rules for writing down and evaluating the Feynman diagrams for the wavevector-and frequency-dependent conductivity tensors.These diagrams rely on the generalized current operator derived in Section II C. Our formalism will closely follow Ref. [29], contrasting only in the inclusion of the wavevector dependence of the external fields and the cur-rent operator.
First let us define the symbology that will compose our diagrams: 1.The free fermion propagator (Matsubara Green's function), G(k), is denoted by .
2. The perturbing vector potential field is denoted by , and carries with it a (Matsubara) frequency and a wavevector.

3.
Where the fermion and photon propagator meet, we place a vertex in the α direction.When corresponding to the output current operator it is denoted by an outgoing photon vertex, , and when corresponding to a perturbing field it is denoted by an incoming vertex, .
We have introduced a four vector notation for Matsubara frequencies, with k = (iν, k) representing the fermionic Matsubara frequencies and momenta, and q = (iω, q) representing the bosonic Matsubara frequencies and momenta.We will also use the shorthand dk ≡ dkdν 1 to denote the integration measure, and q 1 + q 2 ≡ q 12 to denote the componentwise sum of two four vectors.With these components defined, we may write diagrams describing the contributions to ⟨j µ q ⟩ order-by-order in perturbation theory.For concreteness, we will specialize at this point to a free-electron system.As such, the Feynman rules below require us to consider only diagrams with a single fermion loop.We note, however, that our diagrammatic rules can be extended in the presence of interactions, such as from phonons [80], by including additional interaction vertices.To simplify the bookkeeping arising from the negative signs in Eq. ( 65) and from the negative sign of the electron charge, we will formulate our diagrammatic rules in terms of the velocity vertex functions in Eq. (66).
The diagrammatic rules are: 1. Every loop has an output vertex, and most have input vertices.Diagrams contributing the the n-th order response have n + 1 photon lines.All photon lines are incoming at input vertices.At output vertices, exactly one photon line is outgoing.
2. Every loop denotes an integral over both k space as well as a Matsubara sum (which can be made into an integral) in frequency ν.
3. Every closed loop conserves momentum and energy [84].By convention, incoming photon lines carry momentum out of the loop, and outgoing photon lines carry momentum into the loop.
4. Each incoming field vector also contributes a factor of ei ℏiωα , where ω α is the Matsubara frequency for the α-th photon line.

To avoid double counting, only topologically unique
diagrams should be considered.In particular, if the exchange of two or more four-vector and index labels on photon lines (not including the photon corresponding to the output current) does not change a diagram, then the diagram should be divided by the appropriate multiplicity factor [29].In practice, this means that N -th order input diamagnetic current vertices are accompanied by a factor of 1/N !, and N -th order output diamagnetic vertices are accompanied by a factor of 1/(N − 1)!We will see an application of this rule in Sec.VI.
6.A vertex with N photon lines corresponds a factor of (a matrix elements of) v ν N (N ) (k, q 1 , . . ., q N ) defined in Eq. ( 65), where k is the fermion momenta going into the vertex.
After evaluating a particular set of diagrams, we can finally obtain a causal response function by analytically continuing each Matsubara frequency back to a real frequency, iω → ω + iη, for η a positive infinitesimal.

IV. LINEAR RESPONSE
In this section, we will apply the Feynman diagrammatic rules from Sec. III B to calculate the linear conductivity as a function of frequency and wavevector for a general noninteracting periodic solid.First, in Sec.IV A we will use our diagrammatic method to recover the Kubo formula for the conductivity, where we comment on the relationship between our approach and that of Refs.[51,66,68].Then, in Sec.IV B we will show how our formalism relates to the generalized f-sum rule [4,85].
The linear conductivity is given by the sum of the two diagrams in Fig. 1.Using our Feynman rules, this becomes where G n (k) = (ν − ϵ nk ) −1 is the Matsubara Green's function in the energy eigenbasis |ψ nk ⟩ of Eq. ( 43), ϵ nk is the corresponding energy, and we introduce the notation v µν1,...ν N (N +1),n1n2 to denote the matrix elements of the velocity vertex Eq. ( 66) in the basis of Bloch eigenstates |u nk ⟩.
The first term in Eq. ( 73), corresponding to the dia-gram in Fig. 1(a), is the diamagnetic conductivity.We can evaluate Eq. ( 73) and analytically continue to real frequency to find where n F is the Fermi-Dirac distribution and ω + = ω + iη.In the limit q → 0 this reproduces the generalized diamagnetic conductivity of Ref. [29].Here, however, we see that when H k contains non-quadratic momentum dependence-as it will for any semirelativistic system or for any system modeled by a low-energy effective Hamiltonian-the diamagnetic conductivity acquires a nontrivial q dependence.The diamagnetic conductivity serves to regularize the conductivity in the limit ω → 0. We can gain some intuition for this by noting that Fig. 1(a) gives the diamagnetic conductivity in terms of the average of the diamagnetic current Eq. ( 62).The generalized Ward identity from Eq. ( 59) constrains the diamagnetic conductivity.
The bubble diagram in Fig. 1(b) gives the paramagnetic conductivity, corresponding to the second term in Eq. (73).By evaluating the Matsubara integrals in Eq. ( 73) and analytically continuing back to real frequencies, we find that the paramagnetic conductivity is given by where ) is the retarded current-current correlation function in the time domain.In the frequency domain, we can use our expression Eq. ( 53) for the current operator to find where v α ab (k, q) is the charge-conserving velocity, given by ) as indicated in Eq. ( 65).
Our improved definition for the conserved current operator enters into the matrix elements F ab αβ (k, q), while the ratio of Fermi functions to energy differences is universal and arises from the evaluation of the Green's function integral in Eq. ( 73).We can compare our result for the paramagnetic conductivity to that of Refs.[51,68] which uses the midpoint current from Eq. ( 33) instead of the conserved current.This leads to the appearance of a modified matrix element in place of F ab αβ (k, q), with vα k ≡ ∂ Ĥk ∂kα .The subscript "mid" is assigned to this quantity to remind the reader that it was derived using the midpoint definition of the current operator.
Similarly, if the (non-conserved) trapezoid current j from Eq. ( 32) were used in place of j to derive the response, then we would find for the matrix element We emphasize that since Eq. ( 77) is computed using the conserved current in Eq. ( 53), it gives the physicallymeaningful conductivity.We will use Eqs.( 79) and (80) in Sec.V to show how the trapezoid and midpoint currents give quantitatively different predictions for the conductivity as compared to Eq. ( 77).
In particular, Ref. [51] used the (non-conserved) Eq. ( 79) to evaluate the q-dependent Hall response for a two-band tight-binding model, where we expect the use of the conserved current Eq. ( 53) to be important to obtain reliable results.We now argue that, since the current operators Eqs. ( 31), (32), and (33) agree to linear order in q, they yield identical calculations of the paramagnetic conductivity Eq. ( 76) to order |q| 2 .To see this, let us rewrite the current-current correlator in Eq. ( 76) as where we have introduced To proceed, note that for crystalline systems, we know from conservation of (crystal) momentum that the ground state operator of any operator with nonzero momentum (modulo reciprocal lattice vectors) must vanish.We can apply this observation to conclude that the first two terms on the right hand side of Eq. ( 81) vanish unless q = 0. Thus, we have Finally, Taylor expanding Eq. ( 83) for small q shows that the O(|q| 2 ) contribution to the paramagnetic conductivity is determined by the O(|q|) term in δj α q .Since Eqs.(31), (32), and (33) agree to linear order in q, they yield the same prediction for the paramagnetic (and hence the Hall) conductivity.
Note, however, that the diamagnetic vertex Eq. ( 75) derived from the conserved current in Eq. ( 31) differs from the naive q-independent diamagnetic current derived from the nonrelativistic Hamiltonian in Eq. ( 8) and used in Refs.[66][67][68] to calculate the symmetric part of the conductivity from tight-binding models.In particular, since the conserved two-photon vertex j µν −q,q , defined in Eq. ( 62), is an even and symmetric function of q, there will be O(|q| 2 ) contributions to the (symmetric) diamagnetic conductivity for tight binding models that are necessary to ensure charge conservation.

B. The f-Sum Rule and the Diamagnetic Conductivity
In this section, we will show how our definition for the conserved current j q allows us to derive a generalized f-sum rule valid for any effective model.While the derivations in this section do not make use of the explicit form of the current operator, they highlight constraints on the conductivity imposed by charge conservation that are satisfied only when the conserved current operator is used.We begin with a short review of density-density response and the derivation of the f-sum rule [86].Recall that the f-sum rule constrains the spectral weight of the density-density response function We start by defining the spectral density which satisfies the Kramers-Kronig relation The spectral density χ ′′ (ω, q) can also be expressed as the imaginary part of χ(ω, q), and can be directly measured through absorption spectroscopy.
From Eq. ( 86) we can deduce an expansion for the large-ω asymptotics of χ(ω, q) in terms of moments of χ ′′ (ω, q).By Taylor expanding the denominator in the integral in Eq. ( 86) we obtain the asymptotic expansion where is (1/π times) the n-th frequency moment of χ ′′ (ω, q).Note that since the (unprojected) density operators commute, i.e. [ρ q , ρ q ′ ] = 0, we have where in going from Eq. ( 90) to Eq. ( 91) we exchanged the order of the ω and t integrals and used This means that the n = 0 term in the asymptotic expansion Eq. ( 87) vanishes.Using Eqs. ( 87) and ( 93) we thus have asymptotically to leading order as ω → ∞, We can go further and evaluate χ ′′ 1 (q) using the continuity equation to arrive at a general form of the f-sum rule.In particular, inserting the definition Eq. ( 85) into the definition Eq. ( 88) of χ ′′ 1 (q) and integrating by parts we find where we have suppressed the time arguments in the equal-time commutator in the last line.Using the continuity equation, Eq. ( 24), along with the Karplus-Schwinger relationship to simplify commutators of the density operator, we find where we have used Eq. ( 75) to express the sum rule in terms of the average of the two-photon velocity vertex, which is the coefficient of the diamagnetic conductivity σ µν dia (q, ω).For Hamiltonians with quadratic momentum dependence, the diamagnetic conductivity is independent of q and this reproduces the usual f-sum rule.Note that since Eq. ( 96) involves an integral over all frequencies, the sum rule of Eq. ( 97) in the form we have written it relates to the full (unprojected) density operator and diamagnetic conductivity.We defer the exploration of restricted sum rules over a limited frequency range to future work (though see the recent Ref. [87] for progress along these lines).
Our derivation shows how the f-sum rule generalizes to nonzero wavevector for general systems.In particular, we have shown in Sec.II C that the two-photon velocity vertex, j µν q ′ ,−q , which uses the conserved current j q is generally given by Eq. (61).We have also shown how the correlation-correlation operator through the spectral density function relates back to conductivity in Eq. ( 97), and, in App.F, we explore this connection even further through the plasmon dispersion.
For general Hamiltonians, such as those arising in tight-binding approximations to solid state systems, the q dependence of j µν q ′ ,−q will lead to a nontrivial qdependent diamagnetic conductivity, and hence a modification to the f-sum rule according to Eq. ( 97).In Ref. [4] the deviation between the generalized f-sum rule [in the form of Eq. ( 96)] and the electron filling was taken as a measure of goodness of fit for tight-binding parameters.
Here, we take a different point of view: the deviation between the generalized f-sum rule and the electron filling is a consequence of gauge invariance and quantifies a combination of semirelativistic effects as well as information about the truncation of the Hilbert space in an effective low-energy model.Since the spectral density at nonzero q is measurable via absorption spectroscopy, Eq. ( 97) gives an experimental probe of the two-photon velocity vertex, and hence of the conserved current j q [85].

V. APPLICATIONS OF THE LINEAR RESPONSE
Having now derived the complete expression for the qdependent linear conductivity in Eq. ( 73), we will now use it to analyze electromagnetic response in insulators and semimetals.We will begin in Sec.V A by computing the frequency-and wavevector-dependent Hall conductivity in a model of a Weyl semimetal.We will see that for large wavevectors the use of the conserved current Eq. ( 53) in the Kubo formula yields quantitatively different predictions for the Hall response as compared to the trapezoid or midpoint approximations prevalent in the literature.This analysis is applicable to studying the response of the system to spatially-modulated AC electromagnetic fields, such as can be applied using standing waves or gate potentials.Next, in Sec.V B we turn our attention to moiré materials, where the wavevector dependence of σ µν (ω, q) has implications even for optical response due to the large effective lattice constant.We will compute the Kerr angle and ellipticity at oblique incidence as a function of frequency for a model of a moiré Chern insulator in two dimensions, showing that spatial inhomogeneous electric fields can lead to experimentally relevant modifications to the Kerr effect.Finally, in Sec.V C we use our formalism for the conserved current to analyze the magnetic moment and magnetic susceptibility of insulators.

A. Linear Hall Effect in a Weyl Semimetal
As a proof of principle, we will apply our definition for the conserved current Eq. ( 53) to compute the wavevector-dependent Hall conductivity for a toy model of a time-reversal symmetry breaking, inversion symmetric Weyl semimetal with two Weyl points first presented in Ref. [88].We will also compare our predictions with analogous calculations of the Hall conductivity using the (non-conserved) midpoint and trapezoid currents via Eqs.(79).and (80 For future convenience, we define the wavevectordependent anomalous Hall conductivity σ xy anom (q) as We will use our formalism to compute how this topological Hall response is modified in the presence of inhomogeneous electric fields.At this stage, we do not assume that our fields are optical-that is, we do not require that ω = c|q|.As such, our analysis is primarily applicable to the response of the Weyl semimetal to spatially dispersive, slowly varying AC electromagnetic fields.Note that the diamagnetic conductivity Eq. ( 74) is explicitly symmetric, owing to the symmetry of the diamagnetic current vertex in Eq. ( 59).Therefore, when we go to calculate the Hall conductivity (which is antisymmetric), it will not contribute.We can thus focus on the paramagnetic conductivity Eqs. ( 76) and ( 77).

The Weyl Semimetal Model
-π - We will consider a model for a time-reversal symmetry breaking, inversion symmetric Weyl semimetal with tight-binding Hamiltonian given by [92][93][94][95][96][97]: where σ i are the 2 × 2 Pauli matrices with σ 0 defined to be the identity matrix.All values of k are measured in reduced coordinates, i.e. k i ∈ [−π, π] for i = x, y, z.Here k 0 quantifies the separation of the Weyl nodes along the ẑ direction.When |m| > 2t z , this model has a gapped spectrum everywhere in the Brillouin zone except at two Weyl points with (reduced) coordinates (0, 0, ±k 0 ); we will focus on this regime for our analysis.The energy parameters t x , t y , and t z determine the velocities close to the Weyl node in each principal direction.The quantities γ and γ z are tilting parameters that allow for type-II Weyl semimetals to occur.[92].A plot of the spectrum of Eq. ( 101) with k x = k y = 0 is shown in Fig. 2.
Using this model along with our Kubo formula Eq. ( 77), we will examine how the Hall conductivity σ xy H (ω, q) deviates from its topological value as a function of both ω and q.
2. Hall Response when q = (q0, 0, 0) Using Eq. ( 77), we first compute the Hall response in this Weyl semimetal model when the electric field is parallel to the ŷ direction and the wavevector is oriented in the x direction with a magnitude q 0 .The Hall current then flows along the x direction, parallel to the wavevector.Note that in this case, the electric field is transverse, while the measured current is longitudinal.We examine this response because it illustrates the differences between the Hall conductivity calculated with the conserved current and the Hall conductivity calculated with the non-conserved currents commonly used in the literature.
We first illustrate the difference in the three definitions as outlined in Equations ( 77), (79), and (80): the conserved current, the (non-conserved) midpoint definition, and the (non-conserved) trapezoidal definition respectively.
Each of these three different definitions, yielding different predictions for the Hall conductivity, is shown in Fig. 3. Increasing the wavevector generally leads to a more pronounced difference among the three predictions, which makes sense since the midpoint and trapezoid currents are approximations to the λ integration in the definition Eq. ( 53) of the conserved current; these integral Black curves correspond to the conductivity computed with the conserved current jq using Eq. ( 77).For comparison, we also show the Hall conductivity computed with the non-conserved trapezoid current jq (blue) and the non-conserved midpoint current (pink).We see that at large wavevectors, the non-conserved currents give quantitatively different predictions for the Hall conductivity as compared to the conserved current.All plots are generated with natural units (i.e.all energy parameters in terms of tz and e = 1).
approximations become less exact as the wavevector increases.This is further demonstrated in Fig. 4(c) where we can see that the discrepancies between the anomalous Hall conductivity σ xy anom (q) [defined in Eq. ( 100)], computed using each of the three methods increase as q increases.We note also that due to the integration over λ in the definition of the conserved current Eq. ( 53), the anomalous Hall conductivity computed using the conserved current decays as 1/|q| 2 for asymptotically large q.Hence, the continuity equation has an influence on the Hall response that becomes more apparent with increasing wavevector.On the contrary, the non-conserved trapezoidal and midpoint definitions of the current fail to capture this feature, and instead predict an anomalous Hall conductivity that is periodic in q x with period 2π and 4π, respectively.This point is illustrated in Fig. 4(c).
We also find, as expected, that the three current operators predict the same Hall conductivity in the limit q → 0 for all frequencies.Furthermore, from Fig. 4(c) we see that for small |q|, deviations between the conserved and nonconserved predictions for the anomalous Hall conductivity vanish faster than |q| 2 .We see that for small q, the anomalous Hall conductivity behaves as This is consistent with the results of Ref. [51] who found that to quadratic order, the non-conserved midpoint ap-proximation Eq. ( 79) predicts a negative O(q 2 ) correction to the anomalous Hall conductivity arising from bandgeometric effects (recall that to quadratic order the midpoint current and the conserved current yield the same Hall conductivity).
Focusing now on the physically-meaningful conserved current, we show in Figs.4(a) and (b) the real and imaginary part of the Hall conductivity as a function of frequency for five different wavevectors in the x direction.Strikingly, we see that the Hall conductivity at q = (π, 0, 0) is identically zero, i.e. σ xy H (ω, πx) = 0.This is true for computations based on both the conserved and non-conserved currents, and arises due to the simplicity of the model, and can be most clearly understood in terms of the trapezoid method: since v x k depends only on k x and includes only nearest-neighbor hopping, we have Fxy (k, πx) = 0 identically.In a more complicated model with longer-range hopping, we would expect σ xy H (ω, πx) ̸ = 0 generically.Examining Figs.4(a-b), we see that as the magnitude of q = qx is increased, the magnitude of the Hall conductivity decreases.We also see peaks in the Hall response at frequencies corresponding to indirect (nonzero momentum transfer) transitions between occupied and unoccupied states induced by the external electromagnetic field.53) and the Kubo formula Eq. ( 77).(c) shows the real part of the anomalous Hall conductivity [as defined in Eq. ( 100)] as a function of the wavevector changing in the x direction.The parameters m = 4tz, tx = ty = tz, γ = 0, γz = 0, and µ = 0 were used for these calculations.
3. Hall Response when q = (0, 0, q0) Now we take the orientation of the incoming field to be in the ẑ direction, which is parallel to the Weyl node separation.Unlike when q was in the x orientation in the Sec.V A 2, the three different definitions of the current operators yield identical predictions for the Hall conductivity in this case.This is because for our toy model Eq. ( 101), only the matrix elements Eqs. ( 77), (79), and (80) change in the Kubo formula for the conductivity.Since our model has no diagonal hopping, and since the Hall conductivity depends only on the x and y components of the current, we have that F ab xy (k, q) = F ab mid;xy (k, q) = F ab xy (k, q) = F ab xy (k, 0).To change this, we could include diagonal hoppings of the form sin(k i ) → sin(k i )f (k z ) or cos(k i ) → cos(k i )f (k z ) for i ∈ {x, y} and f (k z ) being a general periodic function.Mixing the k x or k y dependence with k z gives rise to the possibility of a Hall conductance where each of the three mentioned methods would then yield different results.
In Figs.5(a-b) we show the real and imaginary part of the Hall conductivity σ xy H (q = q z ẑ, ω) for five values of q z .We point out two relevant wavevector values: q z = π/2 (the momentum transfer of the transition from one Weyl node to the highest energy eigenstate) and q z = π (momentum transfer of the transition transition between Weyl nodes at zero frequency).The q z = π/2 transition is shown in Fig. 2 with the green dotted arrow and the q z = π is illustrated in that same figure with a yellow dotted line.At q z = π/2, we would expect a large Hall response, since the density of states at the band maximum is large.To examine the Hall effect in this orientation, we inspect Fig. 5(a-b) which plots the Hall coefficient as a function of frequency for several wavevectors.At the value q z = π/2, we notice the real part of the Hall effect achieves its maximal value at ω ∼ 2.63t z .
The other interesting value of q z is at π.At q z = π a zero-frequency electric field can excite transitions between the two Weyl nodes.This is, an electron can directly hop from one Weyl node to the other without having to absorb energy.This type of excitation at q = πẑ puts the electron in strongest differential of Berry curvature possible, hopping from one topologically charged Fermi point to the Fermi point with opposite charge.We find that at low frequencies, this leads to a Hall conductivity σ xy H (q = πẑ, ω) with real part that is nearly constant over a wide frequency range, before decaying to zero at high frequencies.We see this in Fig. 5(a).In Fig. 5(b), see that the imaginary part of σ xy H (q = πẑ, ω) the conductivity that is nearly linear at lower frequencies and then decays to zero at high frequencies.
We next consider the anomalous Hall conductivity σ xy anom (q = q z ẑ) as a function of the wavevector.We plot the anomalous Hall conductivity as a function of the q = q z ẑ, in Fig. 5(c).We see that the anomalous Hall conductivity for this model is periodic as a function of q z due to the absence of diagonal hopping.We also see that the Hall conductivity decreases quadratically at small q z , consistent with previous investigations into O(q 2 ) corrections to the Hall conductivity [44,45,98].
Lastly, to study the effect of a finite Fermi surface on the Hall conductivity, we study σ xy H (ω, q = q z ẑ) for our Weyl semimetal as a function of chemical potential.When the Fermi surface volume is finite, we expect a diverse range of intra-and inter-Fermi surface indirect transitions to contribute to the q-dependence of the Hall conductivity.
First We show this in Fig. 6(a), where we plot the anomalous Hall conductivity σ xy anom (q = q z ẑ) as a function of q z for five different values of the chemical potential.We see that when µ = 2t z , the anomalous Hall conductivity is positive at q z = 0 but decreases dramatically as a function of q z .For larger q z the anomalous Hall conductivity reverses sign, and achieves a large negative value at q z = π; the Hall response at q z = π and µ = 2t z is almost an order of magnitude larger than the topological anomalous Hall response σ xy anom (q = 0).To investigate this further, we compute the anomalous Hall conductivity σ xy anom (q = q z ẑ) at ω = 0 and at fixed q z as a function of µ, shown in Fig. 6(b) for q z = 0, and we see that σ xy anom (q z = 0) decreases quadratically as µ → 0, consistent with previous results [99][100][101][102].Importantly, the anomalous Hall response at q z = 0 is proportional to the Berry curvature at the Fermi surface leading to this quadratic dependence on the chemical potential [99][100][101][102].Furthermore, we see in Fig. 6(b) that σ xy anom (q = 0) is always positive for any value of µ.Therefore, since we know from previous literature [36,93] that this contribution is proportional to the Berry curvature, we can establish both the positivity and quadratic behavior as indicators of a Berry curvature-dominated influence at the Fermi surface.
However, as noted previously, σ xy anom (q = πẑ) demonstrates very different behavior as a function of µ, as shown in Fig. 6(c).We see that the anomalous response is positive for µ = 0, and then, as |µ| grows, eventually switches sign and achieves a large negative value at the Lifshitz transition.The behavior of σ xy anom (q = πẑ) near the µ = 2t z arises due to the divergence of the joint density of states (essentially the velocity-independent universal factor in our Kubo formula Eq. ( 77)) at the Lifshitz transition.

B. Kerr Rotation in Moiré Materials
While the wavevector dependence of the optical conductivity can usually be ignored in most solids, the large length scales present in moiré lattice systems mean that spatial inhomogeneities in optical electromagnetic fields may have an appreciable effect on transport.In this section, we will focus on applying our formalism to compute the wavevector dependence of the Kerr effect in 2D systems.Recall that the Kerr effect describes the change in relative angle and ellipticity from the light reflected from a material [22,103].and is a direct probe of time-reversal symmetry breaking [104].As such, the polar Kerr effect has been used as a direct probe of time-reversal symmetry breaking in unconventional superconductors and charge-ordered systems [105][106][107][108][109][110].Additionally, reflectivity studies on 2D materials allow us to probe the frequency-and wavevector-dependence of response functions off the light cone, since for a fixed frequency, the magnitude of the in-plane wavevector can be varied by changing the incidence angle.Due to the intrinsic timereversal symmetry breaking of states in twisted graphene and TMD materials, a study of the Kerr effect in these systems is a natural avenue for future experiments [62].
To this end, we will compute the magnitude of the Kerr effect in a model of a moiré-Chern insulator.First, in Sec.V B 1 we introduce a toy model for a Chern insulator in a moiré system.Then in Sec.V B 2 we compute the Kerr angle and ellipticity for the model, using Eq.(73).Unlike in our analysis of the Hall effect in Sec.V A, here both the paramagnetic and diamagnetic conductivities will play a role.We will compare our results with approximate calculations using the non-conserved trapezoid [Eq.(32)] and midpoint [Eq.(33)] definitions of the current prevalent in the literature, showing that they yield quantitatively distinct predictions for the Kerr angle and ellipticity that could be distinguished in experiment.

Moiré Haldane Chern Insulator Model
For our toy model of a moiré Chern insulator, we start with the spinless Haldane model on a honeycomb lattice [111][112][113][114].The Bravais lattice vectors connecting next-nearest-neighbor honeycomb lattice sites (in the same sublattice) are, in Cartesian coordinates, where a is the moiré lattice constant.The nearest neighbor vectors connecting honeycomb lattice sites in opposite sublattices can similarly be written as We construct a tight-binding model with nearest and next-nearest-neighbor hoppings, as well as a staggered on-site potential and an orbital magnetic flux.We can write the Hamiltonian as H = Here M is the staggered on-site potential, t is the nearestneighbor hopping amplitude, and t 2 is the next-nearest neighbor hopping amplitude.The parameter µ is the chemical potential.Also, the last term in Eq. ( 109) is modified from Haldane's original derivation by the extra term, 3t 2 cos(ϕ).We include this extra term for the convenience of putting the zero of energy in the band gap.Finally, ϕ is the time-reversal symmetry breaking magnetic flux per plaquette.
In this model, we will take |t 2 /t| < 1/3, which forces the bands to not overlap.Additionally, to be in a topologically nontrivial state with Chern number C = ±1, our choice of parameters must satisfy |M/t 2 | < 3 √ 3| sin(ϕ)|.The phase diagram for this model is shown in Fig. 7.
Furthermore, we note this model contains two gapped Dirac points at points K and K ′ , with coordinates

The Momentum-Dependent Kerr Effect of the Haldane
Chern Insulator We will now consider the influence of an optical electric field on the moiré Chern insulator, focusing on the Kerr effect.When the moiré lattice constant is large, the in-plane wavevector |q| ∼ ω/c may be an appreciable fraction of the inverse lattice spacing 1/a, at frequencies near the topological gaps.To model this, we will choose the moiré length and energy scales of our model to be comparable to those of twisted bilayer graphene reported in the literature [62,[115][116][117][118].In particular, we take the lattice constant a = 200 Å, which is achievable at small twist angles in graphene.Additionally, we take t 2 = 1THz, t = 4THz, and M ≈ 4.45THz, such that the relevant energy scale for optical excitations is on the order of 10THz.Additionally, by differentiating the Hamiltonian H defined in Eqs. ( 109)- (112) with respect to k, we find that the matrix elements of the velocity operator are all on the order of ta ∼ 2.6 × 10 −4 c or smaller, where c is the speed of light; this is comparable to the Fermi velocity in a typical metal.
To compute the Kerr response, we consider the ex- perimental geometry shown in Fig. 9.We consider an infinite sample oriented in the xy plane, encapsulated within a dielectric substrate with large index of refraction n R = √ ϵ R = 38 [119].We illuminate the system with linearly polarized light at oblique incidence θ i = 7π/16, E I = E 0 xe −iωt+iq•r .Since the in-plane component of the electric field is continuous across the interface, this means that the in-plane component q y of the wavevector relevant for scattering is, in dimensionless units q y a = an R sin(θ I )ω/(c) ∼ (0.0025THz −1 )ω (113) for our choice of moiré lattice constant.This is several orders of magnitude larger than the typical scale q y a 0 ∼ (10 −6 THz −1 )ω for crystalline systems with typical lattice constants a 0 ∼ 2 Å.Taking the incident photon frequency to be on the order of ∼ 10 THz, which should be reasonable with recent breakthroughs in terahertz spectroscopy experiments, we see q in the moiré system is comparable enough to the lattice spacing to make finite wavevector corrections to optical response large enough to be experimentally measurable; although the effects of nonzero q are still small, they are experimentally non-negligible in moiré systems.Although we use a very large n R substrate appropriate to metamaterials, we see from Eq. ( 113) that a lower dielectric constant can be used provided the moiré lattice constant is similarly increased (e.g., by decreasing the twist angle).Additionally, the effect of the wavevector correction will become more noticeable as the angle of incidence becomes closer to being oblique.
We can now compute the Kerr angle (rotation of the plane of polarization of the reflected wave) and ellipticity (imaginary part of the Kerr angle) for light reflected off our model of a moiré Chern insulator.Since our sample is two-dimensional, its primary influence on the propagation of light is through the boundary conditions that enter Maxwell's equations.In particular, the conductivity tensor σ µν (ω, q = q y ) determines the surface current at the sample, which determines the discontinuity in the magnetic field across the sample and thus the reflection coefficient.We give a full derivation of the Kerr angle and ellipticity in terms of the conductivity tensor in App.E. We can thus use our Kubo formula Eq. ( 73) in terms of the conserved current operator in Eq. ( 53) to compute the Kerr angle and ellipticity.Since the derivation for the Kerr angle and ellipticity involves both the longitudinal and Hall components of the conductivity tensor, we need to include both the diamagnetic and paramagnetic parts of the conductivity, as presented in Section IV A.
In Fig. 10(a) we show the Kerr angle, θ K , and Kerr ellipticity, ϵ K , as functions of frequency, using the Kubo formula in Eq. (73).For comparison, we also include the analogous calculation using the non-conserved trape-  109)-( 112) measured in natural units of inverse lattice constants.The dotted green line denotes the momentum transfer between the K ′ and K points.An excitation between these two points can be excited by adjusting the wavevector of the incoming electromagnetic field.zoid and midpoint definitions of the current operator.As noted at the end of Sec.II, to use the non-conserved currents we introduced an approximate diamagnetic current in order to compute the corresponding q-dependent diamagnetic conductivity.For the parameter values considered here, the topological gaps in the band structure occur at ω = 1.5 THz and another at ω ≈ 19.28 THz, which are indicated by vertical dotted lines in Fig. 10(a).As expected, we see that the Kerr angle rises rapidly (and the ellipticity falls rapidly) near the smaller topological gap, and the Kerr angle nearly vanishes above the larger topological gap.Note also that while the ellipticity and Kerr angles also grow at very low frequencies (ω < 0.5THz), the intensity of the reflected light goes rapidly to zero at low frequencies.We also indicate with a vertical dotted line the frequency ω = 25.60 THz, which corresponds to the highest possible excitation energy in this system.
To highlight the non-negligiblity of including wavevector-dependent corrections to the optical response in moiré systems, in Fig. 10(b) we show the difference between the q-dependent conserved current calculation of the Kerr angle and ellipticity and the Kerr angle and ellipticity computed using the common q → 0 uniform approximation to the conductivity tensor.This shows the calculational error arising from neglecting the wavevector dependence in the Kerr effect.We see that the difference varies on the order of 10 3 arcseconds, which shows that including spatial inhomogeneity when modeling light scattering off the material surface has quantitative effects that cannot be ignored in order to accurately describe the Kerr effect in moiré systems.Moreover, this magnitude of difference is detectable using existing experimental techniques.Thus, this further demonstrates that optical response in moiré materials is sensitive to spatial inhomogeneities in electromagnetic fields.
At the scale of Fig. 10(a), it is difficult to distinguish between the predicted value of θ K and ϵ K from the conserved and non-conserved currents.According to our analysis of the dimensionless in-plane wavevector in Eq. ( 113), we expect the difference in prediction between the three definitions of the current operator to be small but measurable.To demonstrate this, in Fig. 10(c) and  (d) we show the difference between θ K and ϵ K calculated with the conserved current and the midpoint [(c)] and trapezoid [(d)] approximations to the current.We see that the differences are on the order of 1 arcsecond for the midpoint current, to as large as almost 1250 arcseconds (or about 20 arcminutes) for the trapezoid current, which are both within reach of experimental detection.Thus, we expect that the conserved current in Eq. ( 53) will provide a better fit to Kerr effect experiments, especially for terahertz frequencies in moiré lattices.
Next, we turn our attention to the dependence of the Kerr angle on the magnetic flux at fixed frequency.In our tight-binding model, the parameter ϕ controls the magnetic flux through each plaquette; varying ϕ at fixed M allows us to tune between the topological phases with Chern numbers −1, 0, and +1 as indicated in Fig. 7 We will examine signatures of the topological phase transition on the Kerr angle and ellipticity, restricting our attention only to predictions using the conserved current of Eq. (53).
In general, we expect the θ K and ϵ K to both be odd functions of ϕ, since all of θ K , ϵ K , and ϕ are odd under time-reversal.Additionally, since ϕ = 0 and ϕ = π are time-reversal invariant values of the magnetic flux per plaquette, we expect the Kerr angle and ellipticity to vanish for these value of ϕ.We see this reflected in Fig. 11 where we show the Kerr angles and Kerr ellipticities as functions of ϕ for various values of ω.Fig. (a-c) show θ K and ϵ K for ω = 1.5THz (the frequency corresponding to the topological gap), 14THz (an intermediate frequency scale in the model), and 25.6THz (the highest possible excitation energy in the model), respectively.We show with vertical dashed lines the values of ϕ corresponding to the topological phase boundaries in Fig. 7.We see that when the incident photon frequency is on resonance with the topological gap [Fig.11(a)], the Kerr angle and ellipticity peak near the value of ϕ corresponding to the topological phase transition.On the other hand, when   -5  the photon frequency is off-resonance [Fig.11(b,c)] the evolution of the Kerr angle and ellipticity as a function of ϕ becomes more smooth; for very large frequencies [Fig.11(c)], we see that the Kerr angle is approximately zero for all ϕ, and the dependence of the ellipticity angle on ϕ resembles a steepened sinusoidal function.
To further explore the sensitivity of the Kerr angle and ellipticity to the topological phase transition in our model, we compute θ K and ϵ K as functions of ϕ at frequency ω(ϕ) = 2|M ± 3 √ 3t 2 sin(ϕ)|, which is on resonance with the smaller of the two topological band gaps for every ϕ (shown in Fig. 7); we show these results in Figs.11(d) and (e).In these figures, we examine the case when the frequency is varied in step with the size of the Dirac points' topological gaps at ω = 2|M ±3 √ 3t 2 sin(ϕ)|.We see that both the Kerr angle and ellipticity on-resonance is large and sharply peaked at the phase boundary, with maximum values Thus, in order to fully understand and quantitatively model optical experiments being performed on moiré materials, it is necessary to account for wavevector dependence of the optical conductivity.In particular, as moiré systems become larger and as fabrication techniques allow for smaller twist angles and longer-wavelength superlattice potentials, the need to understand how q affects the conductivity will become more imperative for explaining the physics.

C. Magnetic Properties of Insulators
While we have derived our expression for the conductivity tensor in terms of response to an electric field, gauge invariance and Maxwell's equations imply that response to electric and magnetic fields are inextricably linked.In particular, Faraday's law Eq.( 4) implies that a transverse electric field is always accompanied by a magnetic field.Rewriting Eq. ( 4) in Fourier space, we have that B µ,q (ω) = 1 ω ϵ µνλ q ν E λ,q (ω), (114) where ϵ µνλ is the totally antisymmetric Levi-Civita symbol.Additionally, the momentum-dependent conductivity tensor σ µν (ω, q) can be used to derive equilibrium susceptibilities.Recall that in the limit ω → 0 for generic, nonzero q, our perturbing fields A q (ω) and A 0q (ω) become time-independent.The electric and magnetic fields corresponding to these potentials are time-independent, bounded, and oscillate in space with wavevector q.Note that in this limit, the static electric and magnetic fields interacting with the material are purely dependent on the wavevector, q, and therefore offer a promising experimental setup to probe the physics discussed in this manuscript.Thus, the Hamiltonian in the presence of a static ω → 0 electromagnetic perturbation at nonzero q has a static ground state.It follows that taking ω → 0 at fixed nonzero q leaves the system in an equilibrium state [120].This means that there can be no charge transport, and so the ground state current is purely a magnetization current [121] where M q is the magnetization density (i.e.magnetic dipole moment per unit volume).Combining Eqs. ( 114) and ( 115), we will be able to use our formalism for the nonuniform conductivity σ µν (ω, q) to calculate magnetic properties of insulators.First, in Sec.V C 1, we will derive a formula for the magnetic susceptibility tensor in insulators.Next, in Sec.V C 2, we will derive expressions for the magnetic quadrupole moment in finite systems.Finally, in Sec.V C 3 we will show that our formalism is consistent with the Streda formula, for which we will provide a new derivation.We note that the results of this section are completely general, relying only on gauge invariance and the assumption that any two-particle interactions are independent of momentum; they apply equally well to interacting and noninteracting systems.

Magnetic Susceptibility in Insulators
We will derive an expression for the magnetic susceptibility starting from the defining relation for the linear conductivity, While the derivations in this section do not make use of the explicit form of the current operator, they highlight constraints on the conductivity imposed by charge conservation that are satisfied only when the conserved current operator is used.We can rewrite Eq. ( 116) in the A 0 = 0 gauge (which we have used for our derivation of the diagrammatic response in Sec.III) as ⟨j µ q (ω)⟩ = iωσ µν (ω, q)A ν,q (ω).
Now let us suppose that our perturbing field A ν,q is timeindependent, which in Fourier space implies that it consists of only an ω = 0 component.Defining we have that the time-independent current response is given by where the current and vector potential are taken at ω = 0. From Eq. ( 115), we know that R µν (q)A ν,q must be expressible as the cross product of q with a vector.This allows us to write for some tensor α ρν (q).Furthermore, gauge invariance restricts the form of α ρν (q); the average current can only depend on A ν,q through the magnetic field B ν,q = iϵ νλρ q λ A ν,q , and so Combining Eqs. ( 115), ( , we have iϵ µνλ q ν M λ,q = −ϵ µλρ q λ χ ρσ (q)ϵ σδν q δ A q,ν = iϵ µλρ q λ χ ρσ (q)B q,σ .
From Eq. ( 124) we can deduce several general features of the conductivity tensor in the ω → 0 limit.First, we see that if the magnetic susceptibility χ µν is finite, then the conductivity tensor is singular in the ω → 0 limit.In particular, Eq. (124) shows that the magnetic susceptibility determines the weight of the 1/ω pole in the conductivity tensor as ω → 0 at nonzero q.This singularity appears only in the transverse component of the conductivity (i.e, the transverse current response to a transverse field).Counterintuitively, this shows that the transverse conductivity can be divergent as ω → 0 even for an insulator.We also deduce for any system with a finite uniform magnetic susceptibility the low-frequency conductivity satisfies lim ω→0 σ µν (ω, q) ∼ i ω ϵ µλρ q λ χ ρσ ϵ σδν q δ + O(q 3 ), (126) and so the singular part of the low-frequency transverse conductivity depends quadratically on q and vanishes as |q| → 0.
To support these conclusions, we use our formalism to compute lim ω→0 σ µν (ω, q) for our model of a moiré Chern insulator described in Eqs. ( 109)- (112).We show the results in Fig. 12.We see that, in accordance with Eq. ( 124), we have that lim ω→0 Im [ωσ xx (ω, q ŷ)] = lim ω→0 Im [ωσ yy (ω, qx)] ∼ −q 2 χ zz (127) for small q.This figure not only supplicates the q dependence but also the singularity that is linear in ω.Notably, this numerical computation does not assume a form of the magnetic susceptibility, since it just faithfully carries out the conductivity calculation; equality of lim ω→0 ωσ xx (q ŷ, ω) and lim ω→0 ωσ yy (qx, ω) in the limit of small q is a reflection of gauge invariance alone.
We can also invert Eq. ( 124) to derive an expression for the uniform magnetic susceptibility χ µν .First, we note that for a system that conserves energy, the uniform magnetic susceptibility must be a symmetric tensor.To derive this, we follow the logic of Ref. [122], which proved an analogous result for the tensor of elastic moduli.We can consider the change in (free) energy for a system as we slowly move the magnetic field through a closed cycle, Since the total energy is a state function, its change over every closed cycle must be zero.Thus, we deduce that χ µν = χ νµ .Taking two derivatives of Eq. ( 124) with respect to q and making use of the antisymmetry of the Levi-Civita symbol, we find after making use of the symmetry of χ µν that Eq. ( 129) is consistent with and generalizes expressions for the orbital susceptibility that appear in the literature [67,[123][124][125].It allows us to compute the magnetic susceptibility of insulators using our Kubo formula Eq. ( 73) in terms of the conserved current operator.

Magnetic Quadrupole Moments
In this section, we will derive an expression for the magnetic quadrupole moment in the ground state of a fi-  109)-( 112).The real part vanishes.Notice the two curves are similarly quadratic at leading order, and higher-order curvatures separate the curves with increasing q.We used model parameters of ϕ = π/2, M = (3 √ 3 − 3/4)t2, and t = 4t2, which puts the ground state in the topological phase at C = 1 as shown in Fig. 7.
nite system with zero external field using our manifestlyconserved operator j q defined in Eq. ( 31).Using Eq. ( 115), we see that evaluating the average magnetic moment will involve expanding our expression for the current operator in powers of q; we thus expect our conserved current operator will give different predictions for magnetic multipole moments as compared to the nonconserved currents in Eqs.(10) or (11).Since we showed in Sec.II that all definitions of the current operator agree to linear order in q, we begin our analysis with the magnetic quadrupole moment.In the interest of generality, we will here return to consider particles of charge e We start by following Refs.[126][127][128], and take the average of Eq. ( 115) to find in components, To find the higher-order moments (quadrupole, octupole, etc.), we expand both sides of the equation in orders of the wavevector, We point out that in the thermodynamic limit, ⟨j µ 0 ⟩ = 0 by Bloch's theorem [129,130].We now match up the left and right hand sides by orders of q.The magnetic dipole moment M γ0 has been extensively considered in the literature [73,127].We focus instead first on the quadrupolar term, ⟨j µνα 2 ⟩ = iϵ µνγ M α γ1 .Using the antisymmetry of the Levi-Civita symbol, we can rewrite this as Since we can use our definition of the conserved current in Eq. ( 31) to find This implies from Eq. ( 134) that the magnetic quadrupole moment takes the form Given that we know our conserved current in Eq. ( 31) differs from the non-conserved trapezoid [(32)] and midpoint [(33)] at order q 2 , we expect Eq. ( 135) to differ from conventional expressions for the magnetic quadrupole moment that have appeared in the literature.Indeed, if we insert Eq. ( 32) for the non-conserved midpoint current j operator into Eq.( 133), we find the alternative expression which also gives a quadrupole moment of Finally, if we use instead the midpoint current j mid,q of Eq. ( 33), we find the alternative expression By comparing Eqs. ( 135), (137), and ( 139), we see that the conserved current predicts a different value for the ground state magnetic quadrupole moment compared to the nonconserved currents conventionally used in the literature.This suggests that care must be exercised when computing the magnetic quadrupole moment in effective low-energy models, for which the non-conserved currents may not correspond to physically relevant observables.Furthermore, similar discrepancies will appear in the octupole magnetic moment, as well as all higher moments.

Gauge invariance and the Streda formula
Continuing with our study of response to timeindependent fields, we will now specialize to two dimensions, and show how the Streda formula [39,40] arises as a consequence of gauge invariance and Eq.(115).While this result is somewhat tangential to our main argument, it highlights the importance of defining the q-dependent conductivity using the conserved current operator.To begin, recall that although we derived the Kubo formula for the conductivity in the A 0 = 0 gauge, gauge invariance requires that the conductivity tensor also governs the response to gradients of A 0 .In particular, let us consider response to a longitudinal electric field E q (ω) = −iqA 0,q (ω).We have from Eq. ( 116) that We now take the ω → 0 limit at fixed nonzero q.Recall from our discussion preceding Eq. ( 115) that in this limit, the perturbing electric field is static, bounded, and spatially periodic.Thus, the system remains in a perturbed equilibrium state and the only currents that flow are magnetization currents.Hence, we have We can expand Eq. ( 141) to lowest order in q to obtain where we have introduced Finally, we note that a time and space independent scalar potential A 0 is indistinguishable from (the negative) of a uniform variation of the chemical potential −δµ.Taking derivatives on both sides then yields the generalized Streda relation Finally, for an insulator, we know that the longitudinal component σ µν (ω, q)q ν must be analytic and vanishing as q and ω both go to zero [45].This allows us to interchange the order of limits in Eq. ( 143), and identify the susceptibility σ µν (0, 0) with the DC Hall conductivity.We see then that for an insulator, the only nonvanishing component of the response to a longitudinal electric field at zero frequency is the Hall conductivity, and the components of the Hall conductivity are equal to the derivatives of the magnetization with respect to chemical potential.We note that this derivation makes clear the importance of considering an insulating system: it is only for an insulator that we can guarantee the regularity of σ µν (ω, q)q ν .

VI. SECOND-ORDER RESPONSE IN MOIR É MATERIALS
Impressive experimental advances in recent years have spurred a renewed interest in nonlinear electromagnetic responses in (topological) materials.Theoretical and experimental investigations into the nonlinear optical response of electrons in crystals has yielded new insights into band topology and geometry that can be probed using spatially uniform (q → 0) optical fields [13-16, 29-31, 34, 36, 38, 64, 65, 76, 131-133].So far, relatively little attention has been paid to the nonlinear response to finite q electromagnetic fields.As we have shown in Sec.V B, in moiré materials even the wavevector dependence of optical fields may play an important part in determining the electromagnetic response.
Recently, Ref. [43] examined the wavevector-dependent longitudinal nonlinear conductivity at low order in wavevector.Here, however, we will show how our diagrammatic perturbation theory of Sec.III can be used in conjunction with our definition of the conserved current in Eq. ( 53) to compute longitudinal and transverse nonlinear conductivities for arbitrary wavevector.We will focus primarily on the second-order response.In Sec.VI A we will write down and evaluate the Feynman diagrams for the second-order nonlinear conductivity σ µγν (ω 1 , ω 2 , q 1 , q 2 ) that determines the second-order response via This will require using all of the diagrammatic rules of Sec.III B, and will involve a careful treatment of the qdependent diamagnetic current vertices from Sec. II C. Next, in Sec.VI B we will apply our results to compute the second order response of our toy model of a moiré Chern insulator from Eqs. ( 109)- (112).Motivated by experimental considerations, we will focus on spatially rectified (q 1 = −q 2 ) second harmonic generation (ω 1 = ω 2 ).This component of the nonlinear conductivity can be probed using transient grating spectroscopy techniques [134][135][136][137][138]. In such measurements, a pair of coherent lasers are made to interfere at the sample to generate a sinusoidal electric field, whose wavevector can be tuned by changing the wavelength of the beams and the angle between the beams.The nonlinear response of the system to such a finite wavevector perturbation is then determined using a probe pulse.By using light in the visible to extreme ultraviolet range, grating wavevectors on the order 0.001-0.1 Å−1 can be achieved.For moiré systems with lattice constants a ∼ 100 Å, this means that the nonlinear response can be measured over a range of wavevectors spanning multiple moiré Brillouin zones; even for conventional materials with lattice constants on the order of 1 Å, the transient grating wavevector can span a large portion of the first Brillouin zone.In both cases, a formalism such as ours is necessary to compute the nonlinear response of the manifestly conserved current at large q.
A. Second-Order Response We may now turn to higher-order conductivities like the second-order response.We will proceed in a similar fashion as in the linear response case, by using the rules and conventions set out in Secs.II C and III B.
To compute the second-order conductivity, we begin by using the diagrammatic rules in Sec.III B to write down the four relevant Feynman diagrams, shown in Fig. 13.Note that in all cases, q 1 and q 2 flow "into" the fermion loop, while q 12 flows "out" of the fermion loop.We see that diagrams Figs.13(b) and (c) involve the two-photon diamagnetic current vertex, while diagram Fig. 13(a) involves the three photon vertex.Using rule #5 from the Feynman diagram rules in Sec.III B, we see that since Figs.13(b) and (c) are invariant under the exchange of (γ, q 1 ) with (ν, q 2 ), these diagrams enter into the expression for the conductivity with a multiplicity factor of 1/2.
To this end, we will examine a representative example of the input and output vertices in Fig. 13(c) in applying rules #5 and #6 before writing the full mathematical expression from the diagrams.First, the input vertex has a fermion line going into the solid dot vertex, so that k ′ = k + q 1 .The photon line associated with this line carries a value of −q 1 as well, so this input vertex goes as Similarly, we can apply Feynman diagram rules #5 and #6 to the output vertex of the same diagram.Since the fermion line coming into the vertex has momentum k, then k ′ = k.We also see there are two vector potential lines emanating from the open dot vertex: one is going "in" with momentum −q 2 , and the other is going "out" with q 12 .So the output velocity vertex goes as Similarly applying the Feynman rules to each diagram in Fig. 13 and summing the results, we can write the second-order conductivity in terms of Matsubara frequencies as ,n3n1 (k, q 12 ) + (γ, q 1 ) ←→ (ν, q 2 ) .
Each of these diagrams must be symmetrized under the exchange of (γ, q 1 ) with (ν, q 2 ), since as we see from Eq. ( 145), these are merely labels for components of the incident electromagnetic field.We indicate this explicit symmetrization in the last line of Eq. ( 148).Feynman rule #5 in Sec.III B ensures that this explicit symmetrization does not overcount diagrams, by including multiplicity factors like the 1 2 in the second line of Eq. ( 148) from Fig. 13(b).In this format, we can carry out the Matsubara frequency integrals and analytically continue back to real frequencies (using the prescription of Refs.[29,65,80] that iω 1 → ω + 1 , iω 2 → ω + 2 , and ω + 12 ≡ ω 1 + ω 2 + 2iη) in order to find This analytic expression in Eq. ( 149) also allows us to comment on the regularity of the conductivity at low frequencies.In the q 1 , q 2 → 0 limit, it was shown in Refs.[29,43,80] that the diamagnetic vertices in Figs. 13  (a-c) serve a crucial role in regulating the low frequency conductivity; in particular, proper definitions of the diamagnetic current is necessary to ensure that the secondorder conductivity for a band insulator satisfies Using our formalism, we can now examine the behavior of the static (ω 1 , ω 2 → 0) second-order conductivity as a function of q 1 and q 2 .Since we defined our diamagnetic current vertices to satisfy the generalized Ward identity in Eq. ( 64), we are guaranteed that for a band insulator the longitudinal components of the static second order conductivity go to zero as the wavevectors tend to zero.On the other hand, as we showed for the linear conductivity in Sec.V C 1, the transverse components of the second order conductivity can have singularities as we take ω → 0 at fixed q 1 and q 2 .These singularities correspond to second order magnetic and magnetoelectric responses, and are not unphysical.
B. Second-Order Response in a Moiré Chern Insulator: Harmonic Generation in Frequency and Self-Focusing in Wavevector Using Eq. ( 149), we can numerically compute the second-order conductivity as a function of frequency and wavevector for any noninteracting electron system.Typically, calculations of the second-order conductivity have focused on sum frequency generation (where the measured current oscillates as the sum of applied electric field frequencies) or difference frequency generation responses (where the measured current oscillates as the difference of applied electric field frequencies) [29,33,80,[139][140][141]. When we allow for spatially inhomogeneous electric fields, however, the response at second order becomes more complicated.In particular, when the applied electric field varies sinusoidally in space with wavevector q and in time with frequency ω, we can consider sum frequency generation (measured current oscillating at frequency 2ω) and "difference wavevector generation" in wavevector (measured current is spatially uniform with wavevector q − q = 0).Experimentally, this could be measured in a transient grating experiment by looking at the second harmonic generation signal.
To see an example of this in practice, we return to our model of a moiré Chern insulator from Eqs. ( 109)- (112).We can consider the "sum-in-frequency, differencein-wavevector" response.Since previous works [43] have presented a formalism for computing the longitudinal components of the response, we primarily focus on the transverse response σ yyy (ω, ω, qx, −qx) that quantifies the spatially uniform current j y (2ω) that flows in the y direction in response to a transverse y-polarized electric field.Such a response could be probed using inplane polarized THz radiation at oblique incidence.Referring to Fig. 8, we see that the wavevector of incident light is parallel to the vector ∆K = 4π 3 √ 3a separating the K and K ′ points in the moiré Brillouin zone.In Fig. 14 (a) and (b) we plot the real and imaginary parts of σ yyy (ω, ω, qx, −qx) respectively, as a function of frequency for five different values of q between q = 0 and q = |∆K|.We see that for all wavevectors the real part of the response initially grows as a function of frequency, peaks at ω ∼ 0.8t 2 = 0.8THz-such that 2ω is approximately equal to the topological gap-and then decreases at larger frequency.Intriguingly, as a function of wavevector we see that for small values of q the magnitude of the current at ω ∼ 0.8t 2 initially grows as q increases, before again decreasing: the peak value of σ yyy (ω, ω, qx, −qx) is similar for q = 0 and q = |∆K|.
To conclude, we can also study magnetoelectric response by examining the low frequency behavior of σ yyy (ω, ω, qx, −qx).
In line with our discussion in Sec.V C 1, we expect that as ω → 0 for fixed q, the there will be a 1/ω singularity in σ yyy (ω, ω, qx, −qx) whose weight disperses at least quadratically at small q; the weight of this pole quantifies the magnetization current that flows in response to the external magnetic and electric fields (with the magnetic field determined by Eq. ( 4)).At the same time, we expect the longitudinal response ωσ yyy (ω, ω, q ŷ, −q ŷ) to be regular as ω → 0 for any fixed q.We can see both properties in Fig. 14(c).We show the quantities lim ω→0 ωσ yyy (ω, ω, qx, −qx) (in blue) and lim ω→0 ωσ yyy (ω, ω, q ŷ, −q ŷ) (in red) computed for our moiré Chern insulator model.Thus, our formalism for computing transverse nonlinear electromagnetic responses correctly captures the low-energy behavior dictated by Maxwell's equations and gauge invariance discussed in Sec.V C 1 for the linear response.

VII. CONCLUSION
In this work, we have developed a formalism for computing spatially nonuniform (wavevector dependent) linear and nonlinear electromagnetic response functions for condensed matter systems.In Sec.II we introduced a definition of the current operator that can be defined using only the velocity and position operators, independent of detailed knowledge of the microscopic form of the Hamiltonian.Furthermore, unlike the approximations Eqs. ( 11) and (10), our current in Eq. ( 31) is manifestly conserved independent of the microscopic details of the Hamiltonian.It is important to emphasize that while Eqs.(11), (10) and (31) coincide for nonrelativistic Hamiltonians of the form of Eq. ( 8), only our conserved current in Eq. ( 31) remains conserved when Eq. ( 8) is approximated by truncating the Hilbert space.For nonrelativistic systems, we can thus view Eq. ( 31) as a conserving approximation to the total current applicable to effective models such as Wannier-based tight-binding models.This is crucial for calculations involving approximate models of non-superconducting systems with truncated Hilbert spaces, since physically relevant approximations must conserve charge if the total Hamiltonian also conserves charge.Additionally, unlike Eqs. ( 11) and (10), our conserved current (31) also remains conserved when relativistic corrections to the kinetic energy are taken into account.
Next, we showed how Ward identities could be iteratively applied to determine the diamagnetic current operator order-by-order in the electromagnetic vector potential.Using our conserved current as a starting point, we continued in Sec.III to develop a diagrammatic perturbation theory for computing spatially nonuniform linear and nonlinear conductivities, focusing on the case of noninteracting electrons for simplicity.Focusing first on the linear conductivity, we showed how our fully chargeconserving approach implies a generalized f-sum rule relating the density-density response function to the diamagnetic conductivity as a function of wavevector.
We also applied our formalism to compute the wavevector dependent Hall conductivity in a toy model of Weyl semimetal.To connect our formalism to experimentally relevant systems, In Sec.V B we introduced a model for a 2D Chern insulator in a moiré superlattice, such that the wavelength of THz radiation can be a nonnegligible fraction of the moiré lattice scale in certain geometries.We used this model to compute the Kerr angle and ellipticity as a function of frequency for oblique incidence, showing that the effects of spatial nonuniformity are potentially measurable in the next generation of experiments.We also showed how the low frequency transverse conductivity yields insights into the magnetic susceptibility, magnetic quadrupole moment, and Streda formula for insulating systems.Finally, we applied our formalism to study secondorder response to spatially nonuniform, time varying electric fields in two-dimensional systems.We calculated the experimentally relevant spatially uniform secondharmonic generation current that flows in response to a transverse, spatially varying AC electric field for a model of a moiré Chern insulator, which points towards future experimental work on transient grating nonlinear spectroscopy in moiré materials.
Our work opens up several avenues for future theoretical and experimental studies.First, while our explicit calculations in Secs.V and VI were carried out for noninteracting systems, our Feynman diagram formalism in Sec.III can naturally accommodate a treatment of interacting systems by including additional interaction vertices.Our manifestly charge conserving approach to transverse linear and nonlinear conductivity would thus allow a consistent quantitative treatment of magnetoelectric response in models for candidate axionic charge density wave materials such as (TaSe 4 ) 2 I [142,143].Along the same lines, our formalism can be systematically applied to systems with disorder.In the lowest order approximation, the effect of disorder on the nonlinear conductivities may be phenomenologically accounted for by replacing ω + → ω + i/τ , where τ is the quasiparticle lifetime [29].Our formalism, however, can treat disorder scattering more formally by including vertices for scattering of electrons by the disorder potential, and standard diagrammatic techniques [78] for averaging over disorder may be applied.In particular, we expect the (nonlinear) conductivities to depend on both the current vertices defined in Sec.II as well as (nonlinear generalizations of) the diffuson propagator obtained by averaging over disorder.Additionally, our formalism implies additional generalized sum rule relations between nonlinear density response functions and diamagnetic current vertices, generalizing our results of Sec.IV B and making contact with Ref. [3].
As we showed that commonly used approximations for the current operator do not conserve charge within the context of effective models, our work also prompts a reexamination of wavevector-dependent quantities computed using those approximations, such as the magnetic quadrupole moment in models of higher-order topologi-cal insulators [126,128].Furthermore, in calculations of (nonlinear) X-ray scattering processes where both core and valence electronic states must be treated on equal footing, our response formalism based on the conserved current of Eq. ( 31) can be used to ensure Ward identities are obeyed even when relativistic corrections to the kinetic energy cannot be ignored [144].Our work also motivates experimental studies of nonlinear optics in moiré systems, where the effect of spatial inhomogeneity of optical fields may be nonnegligible provided samples are large enough.
Finally, recent advances in superlattice and gate engineering open the door to experimentally studying (nonlinear) response to spatially inhomogeneous electromagnetic field outside of optics.In particular, gate tunable electronic superlattice potentials [145][146][147] could be modulated in time to create time-dependent electromagnetic fields.While recent theoretical [148] and experimental [149,150] progress along these lines has focused on tunable superlattices in untwisted bilayer graphene systems, applications of these techniques to twisted systems would allow the study of response to time-dependent electromagnetic fields with wavevectors comparable to the moiré lattice spacing.We thus emphasize that our results will be directly applicable to the next generation of transport experiments in gate-tunable superlattice devices.
Let us conclude by reemphasizing the importance of the conserved current from Eq. (31).Unlike the approximations of Eqs. ( 10) and (11) to the minimally coupled current in nonrelativistic systems commonly used in the literature, Eq. ( 31) is manifestly conserved for any model Hamiltonian.Since most Hamiltonians of interest in condensed matter systems arise as low energy approximations to (i.e.truncations of the Hilbert space of) complicated many-body semi-relativistic systems, electromagnetic response functions computed within these models will only be faithful approximations to what is experimentally measured if the current used in the calculation is conserved within the Hilbert space of the model.All of Eqs.(10), (11), and (31) accomplish this task to linear order in the wavevector, ensuring that they yield the same approximations in (paramagnetic) conductivities to quadratic order in wavevector.However, only Eq. ( 31) is generally conserved at quadratic order in wavevector and beyond, and only Eq. ( 31) allows for the determination of diamagnetic current vertices via Ward identities, which are essential for properly regularizing the low-frequency conductivity.Thus, we expect that Eq. ( 31) is a necessary starting point for obtaining consistent approximations to response functions beyond quadratic order in wavevector.
where P represents the path ordering (in λ) of the exponential.
Setting τ = 1, multiplying by e A , and using the definition of Eq. (A5) for C(τ ) we have which is equal to Eq. ( 29).

Appendix B: Example: Longitudinal and Transverse
Current Operator for Semi-Relativistic Free Electrons As an example and verification of our expression for the current operator derived from Eq. ( 31), we consider the semi-relativistic free electron Hamiltonian By minimally coupling H SR to a background vector potential, we can derive from Eq. ( 22) that the current operator is given by We can alternatively use our Eq.( 31) to obtain the current j q .The single-particle velocity operator for H SR is Inserting Eq. (B3) into Eq.( 31), we find that our formalism for the current operator yields Finally, using Eqs.( 32) and ( 33), we have that and Let us examine the longitudinal and transverse components of the current using the three possible definitions: Eqs.(B2), (B4) and (B5).For the longitudinal current, we find q•j q,min = q Eq. (B7) shows that the longitudinal component of the current j q , as defined in Eq. ( 31), agrees with the longitudinal component of the current j q,min defined through minimal coupling via Eq.( 21); both j q and j q,min satisfy the continuity equation, Eq. ( 16), as mentioned in the main text.The longitudinal components of jq and j mid,q from Eqs. (B8) and (B9) are distinct, implying that the (non-conserved) jq and j mid,q as defined in Eqs.(32) and (33) do not satisfy the continuity equation.Subtracting Eq. (B7) from Eqs. (B8) and (B9) we find This explicitly shows that jq and j mid,q , as defined in Eqs. ( 32) and (33), are not conserved for the semirelativistic Hamiltonian from Eq. (B1).
Let us now examine the transverse components of the currents Eqs.(B2), (B4), and (B5).Taking the cross product with q, we find that the minimally coupled current satisfies On the other hand, for our λ integral definition of the current from Eqs. ( 31) and (B4) we find This means that although the current j q is conserved, its transverse components differ from that of the minimally coupled current for the semi-relativistic Hamiltonian in Eq. (B1), implying a different definition for the magnetization current.
On the other hand, the transverse current for the midpoint definition of the current Eq. ( 33) is Curiously, if we examine the transverse component of the non-conserved trapezoid current, jq , from Eq. ( 32), we find that for the semi-relativistic system Thus, for the semi-relativistic Hamiltonian, the transverse component of the non-conserved current jq is equal to the transverse component of the minimally coupled current.This means that for the semi-relativistic Hamiltonian, Eq. (B1), we have In other words, for the semi-relativistic Hamiltonian in Eq. (B1), the minimally coupled current in Eq. ( 22) can be written entirely in terms of the velocity operator via Eqs.( 31) and (32).We could, if we desired a modeldependent formulation of the current, apply Eq. (B16) to the semirelativistic Hamiltonian which is applicable to heavy elements, where ⃗ σ is a vector of Pauli matrices acting on electron spin; ⃗ σ • (p i × ∇V (x i )) represents the spin-orbit coupling energy.
If we can assume that our Hamiltonians under study have the form of H SOC in Eq. (B17), then the current defined in Eq. (B16) entirely in terms of the velocity operator can be used to compute both longitudinal and transverse responses.If we cannot assume that Eq. (B17) holds for our system of interest, then we cannot determine the transverse component of the minimally coupled current entirely in terms of the velocity operator without detailed knowledge of the full microscopic Hamiltonian.In such a situation, the only guiding principle is that we define a conserved current, such as Eq. ( 31).In addition, we note that the nonrelativistic approximation to the Dirac equation that yields Eq. (B17) in the absence of an external electromagnetic field does not give a minimally coupled Hamiltonian in the presence of a nonzero vector potential [151].Furthermore, Eq. (B16) fails when higher-order relativistic corrections are taken into account.To see this, we can consider a toy model with the next highest power in momentum according to the semi-relativistic kinetic term: Minimally coupling H ′ to a vector potential via Eq.( 22) gives the minimally coupled current such that Eq. (B20) is not conserved.Importantly, using Eq.(B20) we also find that Thus, for a general Hamiltonian, neither the longitudinal nor the transverse components of the conventional current jq will correctly reproduce the minimally coupled current j q,min .and using Eq. ( 48), we observe that Eq. (C2) reduces to We can simplify Eq. (C4) by rewriting χ γk+(1−λ)q v χ ηk+(1−λ)q in terms of the Bloch Hamiltonian.First, note that where we have introduced In the tight-binding limit where we can simplify Eq. (C5) further as where is the tight-binding embedding matrix, and is the matrix Bloch Hamiltonian in the "non-periodic" gauge with h k+G = V † (G)h k V (G) for any reciprocal lattice vector G. Furthermore, in the tight-binding limit we also have (C13) Notice in going from the second to third line in Eq. (C13), we used an even stricter form of the tight-binding limit than in Eq. (C9): we require not just the first moment but all moments of the position operator to be diagonal in the tight-binding basis.Nevertheless, in this strict tight-binding limit we can combine Eqs.(C4), (C10) and (C13) to find j q → e kαβ 1 0 dλ V ((λ − 1)q)V (k + (1 − λ)q)∂ k h k+(1−λ)q V † (k + (1 − λ)q)V (−λq) αβ c † αk c βk+q Eq. (C14) allows us to compute matrix elements of the conserved current j q in the tight-binding basis entirely in terms of the tight-binding Bloch Hamiltonian h αβ (k).Note that Eq. (C14) could have been anticipated from Eqs. (C4) and (C5) by noting that the nonperiodic basis states in Eq. (C16) have vanishing Berry connection in the strict tight binding limit.
Following a completely analogous set of steps, we find for the nonconserved trapezoid current in Eq. ( 32  In this Appendix, we apply the Karplus-Schwinger relation from Eq. ( 60) to the specific case of the density operator ρ q in order to derive an iterative formula for the diamagnetic current vertices.
First, consider the commutator between a general operator A and the exponential e −iq•x that appears in the Fourier component ρ q .By introducing an auxiliary pa-rameter τ , we can write (D3) Eq. (D3) allows us to rewrite commutators with the density operator that appear in the generalized Ward Identity of Eq. ( 64) in terms of iterated integrals over auxiliary variables λ.
We can also use Eq.(D3) in conjunction with the tightbinding expression in App.C to express the generalized diamagnetic current operators in Sec.II C in the tightbinding basis.In keeping with notation from Sec. C, we will let an overline indicate the operator in the orbital basis.We have Furthermore, in analogy with Eq. ( 66) we can define the tight-binding velocity vertex Note that we can express any of the diagrams and subsequent conductivities in Sec.III B in either the orbital or tight-binding bases.

Appendix E: The Kerr Effect Derivation
In this Appendix, we present a derivation of the Kerr effect for systems with an anisotropic dielectric tensor.
In App.E 1 we start by analyzing electromagnetic scattering from a 3D interface.Then, in App.E 2 we apply these results to the experimentally-relevant situation of scattering from an encapsulated 2D sample considered in Sec.V B.

Kerr Effect Derivation in 3D Materials
Consider an interface between vacuum and a 3D semiinfinite slab of material, with the interface normal to the z-axis.We take the material dielectric tensor ϵ to have the form [152]  We will also make the approximation µ ∼ 1, which is tantamount to including all magnetic response in the transverse dielectric tensor [153][154][155].
We now consider an incident light beam with electric field of the form   e −iωt+iki(y sin(θi)+z cos(θi)) .(E2) Here, R(θ i ) denotes a rotation matrix about the x-axis by angle θ i (while this has no effect on x-polarized light such as in Eq. (E2), we retain it for ease of comparison with later steps in the derivation).Notice that this electric field is linearly polarized and propagating in the direction of θ i , measured from the normal of the plane of the material.The wavevector k i and frequency ω will be related later and E 0 is the real and positive amplitude.
We next consider the reflected field [155,156]:   e −iωt+ik R (y sin(θ R )+z cos(θ R )) .(E3) Before considering the form of the transmitted field, we consider the differential equations it must satisfy.Consider Maxwell's equations in a material [154,157,158]: Here, with the assumption µ ∼ 1, then H T = B T .Also, D T = ϵE T , and we take ρ f = J f = 0. We can now combine Eqs.(E5) and (E7) to obtain a wave equation, The fields that satisfy Eq. (E8) are of the form t ℓ E ℓ e −iω+ik ℓ •r where E ℓ is the ℓ-th eigenvector of ϵ and |k ℓ | is proportional to the corresponding eigenvalue.Given the form of the dielectric tensor in equation (E1), we can solve the eigenvalue equation ϵE = n 2 E [159], to obtain the eigenvalues with normalized (in units of E 0 ) eigenvectors, E p , E m , and E 3 respectively [160,161].The meaning of n 2 ℓ is the refractive index (squared) in the material.
Next, let θ T ℓ be the angle of propagation of each eigenmode of the transmitted wave, measured with respect to the normal.Satisfying the wave equation implies |k ℓ | 2 = n 2 ℓ ω 2 /c 2 .Note that when we assign an angle, this also changes our eigenvalue equation to R(θ T ℓ )ϵR −1 (θ T ℓ ) R(θ T ℓ )E ℓ = n 2 ℓ R(θ T ℓ )E ℓ but the eigenvalues remain the same even if the eigenvectors are now angle-dependent [155].Given the required form of |k ℓ | 2 to be a solution to Eq. (E8), we can now write down the generalized form of the solution of the transmitted electric field [155,158,160]: E T =t p E p e −iωt+ik T 1 (y sin(θ T 1 )+z cos(θ T 1 )) +t m E m e −iωt+ik T 2 (y sin(θ T 2 )+z cos(θ T 2 )) +t 3 E 3 e −iωt+ik T 3 (y sin(θ T 3 )+z cos(θ T 3 )) , (E10) where t ℓ and E ℓ are angle dependent.
To now solve this electromagnetic scattering problem, we can first derive Snell's law for each of the transmitted modes.Snell's law states that the argument in the exponential of all the electric fields must be equal at the boundary.Therefore, we will now have a set of four equations that results in: (E11) k i sin(θ i ) = k T 1 sin(θ T 1 ), (E12) k i sin(θ i ) = k T 2 sin(θ T 2 ), (E13) k i sin(θ i ) = k T 3 sin(θ T 3 ). (E14) It should be noted by the nature of reflection that k R = −k i and θ R = −θ i .
We may now consider the remaining boundary conditions derived by requiring that Eqs.(E4), (E5), (E6), and (E7) are satisfied by components of the incident, reflected, and transmitted waves at the boundary.The Snell's law equations ensure the position-dependent arguments in the exponential are the same, guaranteeing the reflection and transmission coefficients are position (and time) independent.Following Ref. [157], the boundary equations are [156] [ where n = [0, 0, 1], the normal to the 2D material plane.We also observe, as good check, that ∇ • E R = 0 and ∇ • E T = 0 after resolving the coefficients through the boundary conditions.Solving the transmitted coefficients in this basis results in t 3 = 0.
To solve for the Kerr angles, we first rotate our local coordinate system to the frame of the reflected wave, R(θ R )E R .Then we redefine the coefficients as r p = (r x + ir y )/2 and r m = (r y − ir y )/2, which describes the reflection coefficients for the right-hand and left-hand circular polarizations.Each of these coefficients is, in general, complex so we may break up each according to complex polar coordinates, that is r p = |r p |e iαp and r m = |r m |e iαm .Kerr rotation results in a phase shift between right-and left-hand circularly polarized components of the reflected wave.Furthermore the reflected wave can have an amplitude difference between right-and left-hand circularly polarized components.This allows us to define [22,152,154,161] Our derivation expresses the Kerr angle and ellipticity in terms of components of the dielectric tensor.This can be related to the frequency and wavevector-dependent conductivity using [152,154] where I is the 3 × 3 identity matrix.

Kerr Effect Derivation in 2D Materials with a Substrate
In this section we will extend our analysis of Kerr effect to 2D materials.Most of the derivation will remain the same, but there are a few subtle changes.In this setup there will be a vacuum, followed by a 2D thin film, and followed by a substrate which will have the permittivity ϵ R .We will take the ẑ axis to be normal to our 2D film, in the direction of the vacuum.

(E24)
The Snell's law found by comparing the exponentials is also similar to before Since we have a 2D material with nonzero conductivity tensor, a surface current is generated at the interface.As such, the boundary conditions get modified to include this surface current by way of the surface conductivity, j s = σ 2D •E T,2D , which enters into Maxwell's equations as a free current on the surface.Since ∇ • E R,2D = 0 and ∇ • E T,2D = 0 force t 3 = r z = 0, we arrive at four equations arising from two boundary conditions [153,162], Since σ 2D is a property of the 2D sample, it takes the form where the bulk conductivity is related to the surface conductivity by σ 2D = dσ 3D [152].The Kerr rotation and ellipticity can be extracted using Eqs.(E19) and (E20) respectively.
Here, v c (q) is the Fourier transformed Coulomb interaction (in 2D, v c (q) = 2πe 2 /κq and v c (q) = 4πe 2 /κq 2 in 3D) [163,164].The function Π(ω, q) is the polarizability or the electric susceptibility for noninteracting electrons.The electric susceptibility is related to the conductivity through the continuity equation.First consider a longitudinal electric field, E q = −iqA 0 (ω, q), defined in terms of a scalar potential A 0 .The current that flows in response to this field is j(ω, q) = −iσ(ω, q) • qA 0 (ω, q) Additionally, from the continuity equation, −e∂ t n(t, x) + ∇ • j(t, r) = 0. Expressing the density in terms of the response function Π, we find [165]: iωe 2 Π(ω, q) + q • σ(ω, q) • q = 0. (F2) Now Eq.F2 may be substituted into Eq.F1 to obtain the relation ϵ(ω, q) = 1 + i ωe 2 v c (q) (q • σ(ω, q) • q) .(F3) Note that since our conserved current in Eq. ( 31) satisfies the continuity equation, our formalism ensures that Eq. (F3) is obeyed.To solve for the plasmon dispersion, we set Eq. (F3) equal to zero and solve for ω(q), the plasma frequency as a function of q.Since our conserved current Eq. ( 31) obeys the continuity equation, we can use the conductivity Eq. ( 73) to determine the RPA plasmon dispersion in generic models.By contrast, since the midpoint and trapezoid definitions of the current, Eqs. ( 10) and (11), do not obey the continuity equation at large q, they will not give reliable estimates of the plasmon dispersion.

1 2 SpectrumFIG. 2 .
FIG. 2. Energy spectrum of the Weyl semimetal Hamiltonian Eq. (101) along the kx = ky = 0 line in the Brillouin zone, with parameter values m = 4tz, tx = ty = tz, γ = 0, γz = 0, and µ = 0.The dotted lines indicate two types of indirect transitions that can be driven by a time-and space-dependent electric field: the green dotted line shows a jump of of |q| = π/2 in momentum and ω = 2tz in energy, and the dotted yellow line show an jump of |q| = π in momentum space with zero energy transfer ω = 0.

FIG. 3 .
FIG. 3.The Hall conductivity versus the frequency of the perturbing electric field at different values of the wavevector (|q| = π/4, π/2, 3π/4) in the x direction.(a-c) show the real parts of the Hall conductivity for each wavevector while plots (d-f) show the imaginary part.The parameters m = 4tz, tx = ty = tz, γ = 0, γz = 0, and µ = 0 were used for these calculations.Black curves correspond to the conductivity computed with the conserved current jq using Eq.(77).For comparison, we also show the Hall conductivity computed with the non-conserved trapezoid current jq (blue) and the non-conserved midpoint current (pink).We see that at large wavevectors, the non-conserved currents give quantitatively different predictions for the Hall conductivity as compared to the conserved current.All plots are generated with natural units (i.e.all energy parameters in terms of tz and e = 1).

FIG. 4 .
FIG. 4. The real [(a)] and imaginary [(b)] parts of the Hall conductivity as a function of frequency for several values of the wavevector in the x direction.The conductivities in (a) and (b) were computed using the conserved current Eq.(53) and the Kubo formula Eq. (77).(c) shows the real part of the anomalous Hall conductivity [as defined in Eq. (100)] as a function of the wavevector changing in the x direction.The parameters m = 4tz, tx = ty = tz, γ = 0, γz = 0, and µ = 0 were used for these calculations.
, we focus on the anomalous Hall conductivity as a function of wavevector for several different values of the chemical potential.For |µ| < |2t z |, we see from Fig 2 that the Fermi surface will consist of two disjoint pockets centered on each Weyl node.At |µ| = |2t z | there is a Lifshitz transition, where the two Fermi pockets meet; for |µ| > |2t z | there is a single Fermi pocket.Since the density of states and topological charge of the Fermi surface changes drastically at the Lifshitz transition, we expect to see a drastic change in the Hall conductivity near µ = 2t z for all q.
FIG. 5.The real [(a)] and imaginary [(b)] parts of the Hall conductivity as a function of frequency for several values of the wavevector in the ẑ direction.(c) shows the real part of the anomalous Hall conductivity [as defined in Eq. (100)] as a function of the wavevector changing in the ẑ direction.The parameters m = 4tz, tx = ty = tz, γ = 0, γz = 0, and µ = 0 were used for these calculations.

FIG. 6 .
FIG. 6.Chemical potential dependence of the anomalous Hall conductivity [as defined in Eq. (100)] for the Weyl semimetal model.(a) shows the real part of the anomalous Hall response as a function of the wavevector in the ẑ direction, for several different values of µ up to the Lifshitz transition at µ = 2tz.(b) shows the real part of the anomalous Hall response as a function of chemical potential at qz = 0, and (c) shows the real part of the anomalous Hall response as a function of chemical potential at qz = π.The parameters m = 4tz, tx = ty = tz, γ = 0, and γz = 0 were used in these calculations.

2π 3 ( 1 √ 3 , 1 √ 3 , 1 )
1) and 2π 3 (− in units of the inverse moire lattice constant, respectively.The K and K ′ points are shown in the Brillouin zone in Fig. 8.By adjusting M , t 2 , and ϕ, the Dirac nodes can be gapped out.A gap closes at a single Dirac node at each phase boundary as illustrated in Fig. 7.At the K and K ′ points, the band gaps are 2|M − 3 √ 3t 2 sin(ϕ)| and 2|M + 3 √ 3t 2 sin(ϕ)| respectively; we refer to these as the topological gaps.We show the band structure of the model for a representative set of parameters t 2 = 1 THz, t = 4 THz, M = 2.2 THz, and ϕ = π/2.

FIG. 7 .
FIG. 7. The top is the topological phase diagram of the Haldane Chern insulator as described by Eqs.(109)-(112).We set t = 4t2 and t2 = 1 THz to adhere to the non-overlap condition of |t2/t| < 1/3.The Chern number, as illustrated in this figure, is calculated by integrating the Berry curvature across the Brillouin zone.The bottom is the spectrum for the Haldane Chern insulator, capturing both K and K ′ points when ky = 2π/3, which is at kx = ±2π/3 √ 3.At these points, the separation in energies follows as 2|M ± 3 √ 3t2 sin(ϕ)|, which take on separation values of 6 THz and 14.8 THz.

FIG. 8 .
FIG. 8. Brillouin zone for the Chern insulator model in Eqs.(109)-(112) measured in natural units of inverse lattice constants.The dotted green line denotes the momentum transfer between the K ′ and K points.An excitation between these two points can be excited by adjusting the wavevector of the incoming electromagnetic field.

θFIG. 10 .
FIG.10.Kerr effect in the moiré Chern insulator model.(a) shows the Kerr angle θK and ellipticity ϵK as a function of frequency computed using the conserved current, the non-conserved midpoint current, and the nonconserved trapezoid current.(b) shows the difference between the angle and ellipticity computed using conserved current compared to the uniform approximation σµν (q, ω) ≈ σµν (0, ω), effectively showing the importance of q-dependent modifications to the optical response.In (c) we show the differences between the conserved and midpoint current predictions for the Kerr angle and ellipticity.Similarly, in (d) we show the differences between the conserved and trapezoid current predictions for the Kerr angle and ellipticity.These figures used model parameters of M = (3 √ 3 − 3/4)t2, t = 4t2, t2 = 1 THz and ϕ = π/2.Vertical dashed lines denote the two topological gaps and the highest allowed transition frequency.

1 ((((
FIG.12.Plot showing the imaginary part of the diagonal components of the transverse conductivity scaled by the frequency, in the limit of zero frequency for the moiré Chern insulator model introduced in Eqs.(109)-(112).The real part vanishes.Notice the two curves are similarly quadratic at leading order, and higher-order curvatures separate the curves with increasing q.We used model parameters of ϕ = π/2, M = (3 √ 3 − 3/4)t2, and t = 4t2, which puts the ground state in the topological phase at C = 1 as shown in Fig.7.

) that jq = e 2 kαβ∂ k h αβ k+q + ∂ k h αβ
C14), (C17), and (C18) serve as our starting point for computing linear and nonlinear response coefficients from tight-binding models in the main text.
Appendix D: Generalized Integration Formulation of an Operator Using the Karplus-Schwinger Relation