Magnetoelectric polarizability: A microscopic perspective

We extend a field theoretic approach for the investigation of the electronic charge-current density response of crystalline systems to arbitrary applied electromagnetic fields. The approach leads to the introduction of microscopic polarization and magnetization fields, as well as free charge and current densities, the dynamics of which are described by a lattice gauge theory. The spatial averages of such quantities constitute the fields of macroscopic electrodynamics. We implement this formalism to study the orbital electronic response of a class of insulators to applied uniform dc electric and magnetic fields at zero temperature. To first-order in the applied fields, the free charge and current densities vanish; thus the response of the system is characterized by the first-order modifications to the microscopic polarization and magnetization fields. Associated with the dipole moment of the microscopic polarization (magnetization) field is a macroscopic polarization (magnetization), for which we extract various response tensors. We focus on the orbital magnetoelectric polarizability (OMP) tensor, and find the accepted expression as derived from the"modern theory of polarization and magnetization."Since our results are based on the spatial averages of microscopic fields, we can identify the distinct contributions to the OMP tensor from the perspective of this microscopic theory, and we establish the general framework in which extensions to finite frequency can be made.


I. INTRODUCTION
Interest in describing the response of insulators to external electromagnetic fields dates back to the earliest studies of electricity and magnetism.In pioneering work near the start of the twentieth century, Lorentz [1] based his definition of the macroscopic polarization and magnetization fields on a physical picture of molecules with electric and magnetic moments [2], and from that perspective addressed the response of the macroscopic quantities to applied fields.
Near the end of the twentieth century a new approach, called the "modern theory of polarization and magnetization," was introduced [3][4][5][6].Largely focused on the static response of crystalline materials to uniform fields, the microscopic underpinning was now electronic Bloch eigenfunctions [7], or alternately the spatially localized Wannier functions that could be constructed from them.But a macroscopic perspective was taken to define the polarization and magnetization.For example, if one imagined a slow variation in material parameters leading to a macroscopic current density J , the polarization (or at least its change) could be defined through J = dP /dt [8].Thus, instead of a microscopic picture of the underlying position and motion of charges leading to the definition of macroscopic quantities, as it did for Lorentz, links to a microscopic picture follow from the definitions.In the "modern theory," the macroscopic polarization was found to be related to the dipole moment of a Wannier function and its associated nucleus [3].The ambiguity of which nucleus to associate with a given Wannier function -the "closest," or one some number of lattice spacings away?-leads to a "quantum of ambiguity" in the macroscopic polarization itself.Such ambiguities are inherent to the "modern theory," and can generally be related to the behavior and description of charges and currents at the surface of a finite sample [9].
We have recently argued [10] that it is useful to expand upon the approach of Lorentz by introducing microscopic polarization and magnetization fields in bulk crystals, and defining the corresponding macroscopic fields as their spatial averages; in general, microscopic "free" charges and currents, and their spatial averages, are also introduced, and the resulting description takes the form of a generalized lattice gauge theory.The strategy employed is an extension of that used to introduce microscopic polarization and magnetization fields for atoms and molecules [11,12], which itself is an extension of Lorentz' characterization of molecules by a series of multipole moments.Such microscopic polarization and magnetization fields allow for the visualization of electronic dynamics in the sense that perturbative modifications to these microscopic fields arising due to the applied electromagnetic field can be found and exhibited if one has the Wannier functions in hand; existing schemes, which are primarily ab initio based [13,14], can be used to construct such Wannier functions.In the usual "long-wavelength limit" of optics, where the applied electric field is varying in time but its variation in space is neglected, we recover the standard results for crystalline solids.In other instances where comparison with the "modern theory" is possible, such as the response of the polarization to a static or uniform applied electric field, or the response of systems expected to exhibit the quantum anomalous Hall effect, we also find agreement [10].
The approach we implement has the advantage that it can be employed to describe the effects of generally spatially varying, time-dependent applied electromagnetic fields [15].Exploring its characterization of this general-ized optical response, in the linear and nonlinear regimes, is our main program.There is, however, an interesting overlap of our program with that of the "modern theory," and that is in calculating the static response of the polarization to an applied, uniform magnetic field, and the static response of the magnetization to an applied, uniform electric field.This phenomenon is termed the magnetoelectric effect, and for a class of insulators [16] in the "frozen-ion" approximation, where spin contributions are also neglected, both of these responses are characterized by the orbital magnetoelectric polarizability (OMP) tensor [9,[17][18][19][20][21], to first-order in the applied fields.Unlike its generalization to finite frequency, the OMP tensor is non-vanishing only when both spatial-inversion and time-reversal symmetry are broken in the unperturbed system [22], and it is composed of two distinct terms: the Chern-Simons contribution and the cross-gap contribution.The former is isotropic and entirely a property of the subspace of originally occupied states, while the latter involves both occupied and excited states of the unperturbed system.The Chern-Simons contribution has generated particular interest in the literature because of its topological features; there is a discrete ambiguity in its value, which can be used to identify Z 2 topological insulators [9,23,24].As well, while the analytic structure of the cross-gap contribution is of the form one would expect to find from the usual treatment of linear response using a Kubo formalism, the Chern-Simons contribution takes a rather unexpected form.
The expression for the OMP tensor found via the "modern theory" is well established [19,20].In this paper we present a calculation of the OMP tensor within our framework of identifying microscopic polarization and magnetization fields.It is a special case of our general approach, in which by using a set of orthogonal functions that are well-localized spatially one can associate a portion of a total quantity with the point about which each of these functions is localized; a total quantity can be decomposed into "site" contributions.With this, our goal is to formulate the relation of these "site" quantities to the applied electric and magnetic fields evaluated at that site.The OMP tensor is extracted upon taking these electric and magnetic fields to be uniform in space and independent of time.Our results are in complete agreement with those of the "modern theory," as would be expected, but in the process we achieve some insight into the microscopic origin of the distinct contributions to the response.In particular, we can compare our results for the OMP tensor with what would be expected in a "molecular crystal limit," a model in which on each lattice site there is a molecule with orbitals that share no common support with the orbitals of molecules on other lattice sites.And with the development of the formalism presented here we position ourselves to extend this approach to describe material response at finite frequency.
After some preliminary discussion to begin Sec.II, we extend the formalism [10] where necessary in order to calculate the response of the site quantities to arbitrary applied electromagnetic fields.This is a very general development, and only in the later sections do we restrict ourselves to the limit of uniform and static applied electric and magnetic fields.The calculations are made in Sec.III; in Sec.IV we show that the accepted expression for the OMP tensor is reproduced.We also calculate the OMP tensor in the molecular crystal limit by two approaches.The first is a direct molecular physics calculation, and the second is by taking the appropriate limit of our general expressions; they agree, as they should.
We also discuss the nature of both the Chern-Simons and cross-gap contributions from the perspective of this microscopic theory.In Sec.V we conclude.

II. PERTURBATIVE MODIFICATIONS TO THE SINGLE-PARTICLE DENSITY MATRIX
The electronic response of a crystalline insulator is a consequence of the evolution of the fermionic electron field operator, ψ(x, t).We assume that, in the Heisenberg picture, the dynamics of this object is governed by where e = −|e| is the electron charge, and H 0 x, p(x) is the differential operator that governs the dynamics of the electron field in the unperturbed infinite crystal.The presence of an applied classical electromagnetic field, characterized by its vector and scalar potentials, A(x, t) and φ(x, t), has been included through the usual minimal coupling prescription.In writing (2) the independent particle approximation is made, neglecting any interactions other than between the electron field and the applied electric and magnetic fields, which are taken to be the macroscopic Maxwell fields.These external fields are assumed to be applied at times greater than an initial time at which the system is taken to be in its unperturbed ground state, and the expectation values of pairs of (Heisenberg) field operators ψ(x, t) and their adjoints in the unperturbed ground state are used to construct the minimal coupling Green functions [10].
Implementing the frozen-ion approximation, we take where V (x) is the spatially periodic lattice potential that characterizes the crystal structure and satisfies V (x) = V (x + R) for all Bravais lattice vectors R, and In (4) we have allowed for the presence of an "internal," static, cell-periodic magnetic field described by the vector potential A static (x), where A static (x) = A static (x + R).
The inclusion of such an "internal" field respects the discrete translational symmetry of the crystal, but generically leads to a Hamiltonian (3) with broken time-reversal symmetry, which will be important in what follows.In future publications we plan to include both the spin-orbit and Coulomb interactions that we neglect here.
In periodic systems, a set of exponentially localized Wannier functions (ELWFs), {W αR (x) ≡ x|αR }, can generally be constructed [25][26][27][28][29] via where Ω uc is the unit cell volume, and are eigenfunctions of (3) that are normalized over the infinite crystal such that ψ mk |ψ nk = δ nm δ(k − k ).A periodic gauge choice is made such that the Bloch eigenvectors |ψ nk and the unitary matrix U (k) are periodic over the first Brillouin zone [30].Associated with each Bloch eigenfunction ψ nk (x) is an energy E nk and a cellperiodic function u nk (x) ≡ x|nk satisfying the orthogonality relation (mk|nk) = δ nm ; we adopt the notation for functions g(x) ≡ x|g and h(x) ≡ x|h that are periodic over a unit cell, where the integration is over any unit cell.Also, we restrict our study to three-dimensional systems [31].Here and below n is a band index, α is a type index, and k denotes a crystal-momentum within the first Brillouin zone.
In this paper we initiate our considerations with the zero temperature ground state of an insulator, and consider the class of insulators for which the set of occupied (unoccupied) energy eigenfunctions map to a set of occupied (unoccupied) ELWFs.This class includes ordinary insulators [32] and Z 2 topological insulators [27], but excludes, for example, Chern insulators [29].Therefore, U nα (k) = 0 only if both n and α label occupied or unoccupied band and orbital types.We associate a filling factor f n with band n that is either 0 or 1, and a filling factor f α with orbital type α that is also either 0 or 1; the latter can be inferred directly from the occupancy of the bands used in the construction of a particular ELWF.
The ELWFs (5) can generally be expressed as where another set of cell-periodic functions, { x|αk }, have been introduced and are formed via Here the { x|αk } are in general not the cell-periodic parts of energy eigenfunctions.The sets of vectors {|αk } and {|nk } are related by a (in general) multiband gauge transformation characterized by U (k) [9].Although we use Roman and Greek subscripts in the notation U nα (k), as we will primarily be considering a transformation from the cell-periodic part of energy eigenfunctions to cellperiodic functions that are not associated with energy eigenfunctions, this need not necessarily hold.Indeed, the simplest type of gauge transformation, although of course it would not generally lead to |αk associated with ELWFs, is one that involves each band individually, and is of the form U nα (k) = δ nα e −iλn (k) , where, for a given n and k, λ n (k) ∈ R. Rather than a special limit of the general multiband transformation, this could be considered as simply a new choice of Bloch eigenstates; under Bloch's theorem, energy eigenstates are uniquely defined only within a k-dependent phase [26], even at k points where there is no degeneracy.Here, however, it is considered as one type of U nα (k), associated with a gauge transformation of the U(1) type.That is, we consider the vectors {|nk } fixed at the start, and use the term "gauge dependent" for quantities that depend generally on the U nα (k) and their derivatives, including U nα (k) of the U(1) type.
The ELWFs are an important element of our approach, because we use them to introduce "site" quantities, and define the macroscopic polarization and magnetization in terms of their moments, as we discuss in detail below.Of course, not all gauge transformations of an initial set {|nk } will lead to ELWFs via (5), as indicated by the example given above.Nonetheless, the response tensors, such as α il , must be such that the resulting charge and current densities in the bulk are not only invariant with respect to the choices of U nα (k) that lead to ELWFs, but in fact to all choices of U nα (k); that is, the charge and current densities in the bulk must be gauge invariant.This is plausible because it would be possiblealthough we would argue much less convenient and less interesting physically -to calculate those charge and current densities directly from the minimal coupling Hamiltonian without ever introducing Wannier functions.And we shall see that this gauge invariance does indeed hold.
A useful identity [26] is where is the non-Abelian Berry connection associated with the set {|αk }; here and below, superscript indices indicate Cartesian components [33], repeated Cartesian components are summed over, and we adopt the shorthand ∂ a ≡ ∂/∂k a .The object ( 9) is related to the non-Abelian Berry connection associated with the set {|nk }, Here we have defined the Hermitian matrix W a [9], populated by elements which, for the class of insulators we consider, is nonzero only if f m = f n .Under the aforementioned periodic gauge choice, all objects appearing in (11) are periodic over the first Brillouin zone.In what follows, the k dependence of the preceding objects is usually kept implicit.
A consequence of the Hamiltonian (3) and the resulting dynamics (2) of the electron field operator is that the differential operators associated with the spatial components of the conserved current take the form in the usual fashion [34], where J a x, p(x) are the analogous differential operators arising for the unperturbed system.Another useful identity is where the matrix elements are This can be shown by breaking the integral in ( 14) into the sum of integrals over unit cells; the sum over lattice vectors yields the Dirac delta function in (14), and the expression for (15) follows from the form of the integral over the unit cell and the use of (10).Indeed, a more general form of (14,15) can be derived involving the matrix elements of p a (x) in the basis of the cell-periodic functions x|αk using the same strategy [35].
In previous work [10] we considered the calculation of the expectation values of the electronic charge and current density operators for a crystalline insulator perturbed by an electromagnetic field.Noting that the lesser, equal time single-particle Green function can be used to find such quantities, we employed a set of spatially localized, "adjusted Wannier functions" as a basis in which to expand the electron field operator in an effort to associate portions of the full electronic Green function with individual lattice sites.Upon identifying such "site" Green functions, and thereby identifying "site" charge and current densities, we defined microscopic polarization and magnetization fields associated with each lattice site using the same functions that are used in atomic and molecular physics to relate the microscopic polarization and magnetization fields of atoms and molecules to the microscopic charge and current densities [36]; we call these functions "relators."The full polarization and magnetization fields are given by summing the respective site contributions.We then insisted that these microscopic polarization and magnetization fields, together with the electronic charge and current density expectation values, satisfy the expressions arising in classical macroscopic electrodynamics relating such quantities.This led to the identification of microscopic "free" electronic charge and current densities, which take predictable forms.At zero temperature the first-order modifications to both the free charge and current densities due to an applied electromagnetic fields vanish for the class of insulators considered here, even for applied fields in the x-ray regime [10].As a consequence, the first-order perturbative modifications to the expectation values of the electronic charge and current density operators resulting from an applied electromagnetic field can be found directly from the firstorder modifications to the microscopic polarization and magnetization fields.
A quantity central to the calculation of both the microscopic polarization and magnetization fields is the singleparticle density matrix, η αR;βR (t).Thus, a starting point in describing the effect of an applied electromagnetic field to a crystalline insulator is identifying how the applied fields affect this object.The single-particle density matrix evolves according to [10] i ∂η αR;βR (t) ∂t = µνR1R2 F µR1;νR2 αR;βR (t)η µR1;νR2 (t), (16) where The definitions of ∆(R 1 , R 2 , . . ., R N ; t), Ω 0 R (R, t), and HαR;µR1 (t) are as given earlier [10], and are provided in the Appendices.The first of these quantities is related to the magnetic flux through the surface generated by connecting the points (R 1 , R 2 , . . ., R N ) with straight lines, when the usual choice of straight-line paths for the relators is adopted.The second is related to the electric field along the path connecting points (R , R).The third quantity can be understood as a generalized "hopping" matrix element.Each of terms appearing in (17) is gauge invariant in the electromagnetic sense, and consequently so is η αR;βR (t).
In Appendix A we show that e i∆(R,R ,R ;t) HαR;µR (t) = e i∆(R,Ra,R ,R ;t) HαR;µR (R a , t) and for any lattice sites R a and R b .We have defined HµR1;νR2 (R a , t) ≡ χ * µR1 (x, t)e i∆(R1,x,Ra;t) H Ra (x, t) + 2 ∂∆(R a , x, R 2 ; t) ∂t e i∆(Ra,x,R2;t) χ νR2 (x, t)dx and where as before [10].Here Ω Ra (x, t) is related to the magnetic field along the path connecting points (R a , x), and is defined in Appendix C. The functions in the set {χ αR (x, t)} are generally not orthogonal, but they are related to the mutually orthogonal "adjusted Wannier functions" introduced earlier [10].In the limit of a weak applied magnetic field, a perturbative expansion can be constructed for each χ αR (x, t) [10], and the lowest order terms are Choosing R a = R b , one can then re-express (17) as To prepare for later perturbative analysis, we expand all quantities in powers of the applied electromagnetic field such that , where the superscript (0) denotes the contribution to the total quantity that is independent of the applied fields, the superscript (1) denotes the contribution that is linear in the electric and magnetic fields, and so on.Using (16), and equating terms appearing with the same powers of the applied fields, the zeroth-order term η (0) αR;βR is found to satisfy the same equation of motion as the unperturbed single-particle density matrix, and so while from (20) it is found that Next, implementing the usual Fourier series analysis via the first-order modification to the single-particle density matrix due to an applied electromagnetic field can be identified, via (16), as recall that for the class of insulators of interest here there are well defined filling factors associated with the orbital type indices, but in general this is not so.In identifying (28), we have also introduced Often in optics one is interested in applied fields that vary little over electron correlation lengths which, for the class of insulators we consider, are on the order of the lattice constant.In Appendix C we show that if this approximation is made, and if |y − x| and |z − x| are on the order of lattice constants, then where E(y, t) is the applied electric field, B(y, t) is the applied magnetic field, and Furthermore, under this approximation, and implementing ( 3) and ( 27), the first-order modification to ( 21) is found to be In the above we have used straight-line paths in the relators, and in what follows, this choice is always made.
Eq. ( 28)-( 34) make clear the motivation for introducing the arbitrary lattice site R a ; it serves as a reference site for the applied electromagnetic field.This will prove useful when considering the response of a quantity associated with site R to a spatially varying electromagnetic field; choosing R a = R, the modification to that site quantity is related to the applied field strength at that site.
However, in this paper we restrict ourselves to uniform applied electric and magnetic fields in the dc limit.Thus F ab (x, t) vanishes and η (1) αR;βR (ω) = 0 only if ω = 0. We keep only the non-vanishing, first-order perturbative modifications arising from the applied electric and magnetic fields, which we denote by η (E) αR;βR and η Implementing (28), the first-order perturbative modification to the single-particle density matrix due to a uniform applied dc electric field is found to be (Appendix B), where E ≡ E(R a , ω = 0), for any R a , is now the uniform dc electric field.The result ( 36) is consistent with previous work [10], where this expression was derived via a different method.Notably ( 36) is written as a single Brillouin zone integral, unlike (28).This feature is expected upon comparison to the usual perturbative treatment, in which, when spatial variation of the electric field is neglected, a k-conserving interaction term arises that results in perturbative modifications being given by single Brillouin zone integrals [37].While here we assume the applied electric field is uniform, more generally this reduction to a single k integral holds as a consequence of the approximation that the applied electric field is to vary little over electron correlation lengths, allowing a Taylor series expansion about each lattice site.In the limit of uniform applied fields of interest here, the R a dependence of (36) vanishes.
There are two distinct contributions to the corresponding first-order perturbative modification arising from a uniform applied dc magnetic field, where we have defined and where lab is the Levi-Civita symbol; B ≡ B(R a , ω = 0), for any R a , is the uniform dc magnetic field.We mention that the two terms appearing in (37) are not simply the individual contributions of the terms of (28).Moreover, in the first-order perturbative modifications of the quantities considered below, the first term of (37) gives rise to gauge invariant contributions, while the final term gives rise to gauge dependent contributions [38].

III. ELECTRONIC RESPONSE TO APPLIED UNIFORM DC E AND B FIELDS
We now use the perturbative modifications to the single-particle density matrix (35) to calculate the firstorder response of the electric and magnetic dipole moments -found from the microscopic polarization and magnetization fields, respectively -to applied uniform dc electric and magnetic fields.In this section we restrict earlier [10], more general expressions to this limit.Thus a function previously written in terms of frequency components g(ω) (27), will simply be given by the single non-vanishing component g ≡ g(ω = 0).Furthermore, as previously mentioned, to first-order in the applied fields, the free charge and current densities vanish at zero temperature for the class of insulators we consider [10], as would be expected physically, and so those quantities do not appear here.

A. Summary of formalism
The microscopic polarization and magnetization fields can be decomposed into site contributions [10], such that where the sums range over all Bravais lattice vectors R.
In the limit of applied uniform dc electric and magnetic fields treated perturbatively, the discrete translational symmetry that was present in the unperturbed Hamiltonian ( 3) is now lost (cf.( 2)), but the polarization and magnetization fields associated with each lattice site are still physically equivalent; by this we mean for any R and R .This is a result of the fact that it is the electric and magnetic fields, not the vector and scalar potentials, that enter in the expressions that follow.Thus, in its perturbative response, the system retains its periodic nature in the limit of applied uniform dc fields.Each site polarization field, p R (x), is related to the electronic charge density associated with that site, ρ R (x), via the "relator" s i (x; y, R) (see Appendix C), as shown previously [10].It is given by where above and below the sums range over all lattice vectors and orbital types, and Meanwhile, there are two contributions to each site magnetization field The first of these, mR (x), corresponds to the "local" or "atomic-like" contribution to each site magnetization field, and is related to the electronic current density associated with that site, j R (x), via the "relator" α ib (x; y, R) (see Appendix C) [39].It is given by where It is not obvious that there is another contribution to each site magnetization field, as each atomic-like contribution is found from a site electronic current density, and these collectively compose the total current.However, in extended systems where the charge-current density associated with each site need not be conserved, there is in general an additional term, mR (x) [10].Adopting the terminology of the "modern theory," it corresponds to the "itinerant" contribution to each site magnetization field, and is given by Here with where and Typically it is the response of the electric and magnetic dipole moments to external perturbations that one is interested in studying; these correspond to the spatial integrals of the microscopic polarization and magnetization fields [10] respectively, From ( 40) we find and ν R = νR + νR , where, from (43), is the atomic-like contribution to the magnetic dipole moment and, from (45), is the itinerant contribution to the magnetic dipole moment.The macroscopic polarization and magnetization can be found from their respective dipole moments introduced above, and are taken to be In the limit of applied uniform dc fields, both the electric and magnetic dipole moments are independent of R, and as a consequence, the macroscopic polarization and magnetization are uniform.

B. Unperturbed expressions
We begin by confirming that our microscopic treatment yields the standard expressions for the unperturbed ground state macroscopic polarization and magnetization, more usually constructed from macroscopic arguments [40].While (51,52,53) have been defined to include only valence and conduction electron contributions, the contributions from ion cores can be identified as well (see Mahon et al. [10]).However, we focus only on the former contributions here.
Expanding (51) in powers of the applied electromagnetic field, we find the zeroth-order term to be where we have used ( 8)- (11).This corresponds to the unperturbed ground state electric dipole moment, and upon implementing (54), the usual expression [3] for P (0) is reproduced.It is well known that the unperturbed macroscopic polarization is unique modulo a "quantum of ambiguity."This ambiguity originates from the gauge dependence of (55), and it has been shown that the gauge dependent term of (55) contributes only to this "quantum" [3].Importantly, it is only the diagonal elements of the W i matrix that appear in (55), and as a result, even a U(1) gauge transformation can give rise to this "quantum of ambiguity."This is discussed further in Sec.IV.Turning to the magnetic dipole moment, we expand (52) and (53) in powers of the applied electromagnetic field.The zeroth-order terms are found to be and which, together, form the unperturbed ground state magnetic dipole moment, ν and (57) are "multiband gauge dependent."By this we mean that only if there are intersections amongst occupied bands is it gauge dependent; in the isolated valence band limit both (56) and (57) become gauge invariant.Nonetheless, even if the occupied bands are not isolated the sum of (56,57) is generally a gauge invariant quantity and thus there is no ambiguity in the value of the unperturbed macroscopic magnetization.Implementing (54), the usual expression [4,5] for M (0) is reproduced.

C. First-order perturbative modifications
We now turn to the first-order modifications to µ R and ν R and thus, through (54), to P and M due to an applied electromagnetic field.Generally, the "site quantities" we consider are of the form where η αR ;βR is the single-particle density matrix, and where Λ i βR ;αR (R), which we call a "site quantity matrix element," is given by see (41,51) for µ R , and (44,46,52,53) for ν R .

Dynamical and compositional modifications
The first-order modification to (58) due to an applied electromagnetic field thus has two types of contributions, the first arising from the combination of the unperturbed expression for Λ i βR ;αR (R) and the first-order modification to η αR ;βR due to the applied fields, and the second from the combination of the first-order modification to Λ i βR ;αR (R) due to the applied fields and the unperturbed expression for η αR ;βR , (25), We refer to the first contribution (61) as "dynamical" because it arises from modifications to the single-particle density matrix, which characterizes the electron hopping amplitude between various lattice sites and orbital types, due to the applied fields.Notably in the response of the site quantity associated with R this lattice site always appears as at least one lattice site in the single-particle density matrix components.This is expected physically; the site quantity associated with R is affected by electrons moving between different orbital types at that lattice site, and by electrons moving from the region nearest R to regions nearest other R .Contributions to (61) arising from η βR ;αR and η αR;βR , with R = R, are a consequence of the extended nature of the system, in which electrons are not confined to regions of space; such contributions vanish in the limit that the crystalline solid is considered simply as a periodic array of "isolated molecules," which we call the "molecular crystal limit" (see Sec. IV B).Conversely, contributions to (61) arising from η αR;βR take the form of single-site modifications.
The second contribution (62) to the first-order modification of a site quantity associated with R depends on the first-order modification to the site quantity matrix element associated only with lattice site R, Λ i(1) αR;αR , and with orbital types α that are originally occupied.It is not associated with any change in the single-particle density matrix, but rather with the dependence of the associated site quantity matrix elements on the applied fields themselves; thus we call it a "compositional" modification.It leads to a dependence of the moments µ R and ν R on those fields, even though in this contribution (62) the populations remain as they were before the fields were applied.There is a familiar analog to this in the response of an atom to an applied magnetic field.Considering a single electron, the initial operator for the magnetic dipole moment ν atom = (e/2mc)X × P, where here X and P are the position and momentum operators of the electron, becomes ν atom → (e/2mc)X × (P − eA(X)/c) when a magnetic field is applied.For a uniform applied magnetic field we can take A(X) = (B × X)/2, giving The second term gives a contribution when the expectation value is taken, even in the ground state (say a 1s orbital), and gives the diamagnetic response of the atom.The contributions from the compositional modification to (60) are of this form.Nonetheless, in the extended systems we consider it is important to note that during the perturbative analysis many lattice sites and all orbital types may be involved; for instance, observe that ( 46)-(49) would be used in constructing (53).
We also distinguish between the first-order modifications arising from the applied electric field and the applied magnetic field, such that where each of these modifications is composed of a dynamical and a compositional term, 2. Induced polarization a. Response to electric field.We begin by considering modifications due to an applied electric field.From (41), and the fact that χ αR (x) only depends on the applied magnetic field and ∆(R, x, R) = 0, it is clear that and so there is no compositional modification to µ R .The first-order modification is entirely dynamical, as described above, and given by Implementing (54), we find the usual result from perturbation theory [41].This modification is gauge invariant, in that the final line of (66) is independent of the unitary transformation U nα (k) [42].b.Response to magnetic field.The first-order modification due to an applied magnetic field has nonvanishing dynamical and compositional modifications.Considering first the compositional modification, we find We later simplify this term; notably it is gauge dependent.Recalling the magnetic field modifications to the single-particle density matrix (37), we find which is also gauge dependent.The two sets of terms appearing in the final equality of (68) originate individually from the two terms of (37).In going from the first to the final equality we have implemented (8,11), and used W i mn = 0 only if f m = f n , which holds for the class of insulators considered here.Similar arguments are used in the following subsection when finding (70,72).In Sec.IV we explicitly combine (67) and (68), and show that the usual OMP tensor is reproduced [20].

Induced magnetization
In this work, we only consider modifications to the site magnetic dipole moment arising from an applied electric field.We defer to a later study the response of the magnetic dipole moment to a magnetic field, as considered in this framework.
a. Response of atomic-like contribution to electric field.This involves the first-order modification to (52) due to an applied electric field.There is no compositional modification, as from (44) it is clear that following the argument leading to (65), and so this modification is entirely dynamical.Using (36) we find b. Response of itinerant contribution to electric field.This involves the first-order modification to (53) due to an applied electric field, and has non-vanishing modifications of both compositional and dynamical origin.The compositional modification is and the dynamical modification is Both ( 70) and (71) are gauge dependent in general, while (72) is multiband gauge dependent.Very generally, there is a simplification that occurs when (70,71,72) are summed to form the total magnetic dipole moment: the term appearing in the final line of (72) cancels with the term appearing in the final line of (70), and as a result the gauge dependent terms appearing in the total ν (E) R do not explicitly depend on the energies E nk .

A. Constructing the OMP tensor
The OMP tensor, which characterizes the first-order response of the macroscopic polarization to an external uniform dc magnetic field, is defined through and, from (54), (67), and (68), is found to be after some manipulation of the band indices.If we now define an analogous tensor characterizing the first-order response of the macroscopic magnetization to an applied uniform dc electric field, we find that, upon combining (70), (71), and (72), this response is described by the same α tensor introduced above, but with the order of the indices switched, such that This is the usual result in the ω E gap limit, which is effectively the condition we have initially assumed.
In the previous section we made a concerted effort to identify dynamical and compositional contributions to the various first-order modifications; distinguishing between these proves useful here.Because the α il tensor characterizes the linear response of both P i to B l , and M l to E i , we can focus on the terminology associated with only the magnetization.In what follows, we illustrate how the three types of non-vanishing terms -the dynamical atomic-like, the compositional itinerant, and the dynamical itinerant modifications -combine to give an OMP tensor having the usual form where α il G contains only cross-gap contributions to the response and α il CS , the Chern-Simons contribution, is a property of the originally occupied subspace.We find that the cross-gap response, α il G , originates from a combination of the dynamical atomic-like (70) and dynamical itinerant (72) terms, whereas α il CS has contributions from the dynamical atomic-like term as well, but also from the compositional itinerant term (71).
In order to proceed, we first re-express (73) in terms the cell-periodic functions x|nk , and re-write the sums to be over occupied states (the set {|vk }) and unoccupied states (the set {|ck }).Taking the term in (73) involving the ratio of energy differences to start, which can be traced back to (70), and adopting the shorthand |n ≡ |nk , we find The first set of ... 's are identified as cross-gap contributions and will be included in α il G .The second set of ... 's, together with the penultimate and final terms of (73), form α il CS ; the penultimate term of (73) originates from the gauge dependent term arising in the dynamical atomic-like modification (70) that does not explicitly depend on energy, and the final term of (73) originates from the compositional itinerant modification (71).We find (Appendix D), which is the usual Chern-Simons contribution to the OMP tensor [17][18][19][20], and is multiband gauge dependent in the sense introduced after Eq. ( 57).Due to this gauge dependence, this contribution to the OMP tensor is multivalued, but, like P (0) , has been shown to be unique modulo a quantum of indeterminacy [9].Furthermore, this contribution is isotropic, and the corresponding quantity vanishes in systems of spatial dimension less than three.
The remaining terms compose α il G , in accordance with (74).These terms originate from the dynamical atomic-like modification (70), and the dynamical itinerant modification (72).We find which is in agreement with the usual expression for the cross-gap contribution [19,20].Importantly this contribution is found to be gauge invariant; this is a result of the fact that the multiband gauge dependent terms appearing in (70) and (72) cancelled one another.It has been pointed out that the net gauge dependence of the OMP tensor, through the Chern-Simons contribution, has no effect on the charge-current density in the bulk [18,20].We present a slightly different formulation of that argument here, but starting with the same assumption used in earlier arguments: While the derivation we have presented here, as well as that in the approach of the "modern theory," holds strictly only for uniform dc applied fields, we can expect that for electric and magnetic fields varying very slowly in both space and time the same response tensors α il CS and α il G can be used to good approximation [43].For such slowly varying fields, the first-order response to the applied fields in the bulk of a medium can be expressed as where we have taken α il CS = δ il α CS , and at the level of linear response we have , where χ il E and χ il B are the usual linear electric and magnetic susceptibilities.Inside the bulk material we then immediately find that the macroscopic charge and current densities that result, ρ (1) (x, t) = −∇ • P (1) (x, t), can also be written as That is, the contributions from the Chern-Simons response coefficient, α CS , completely cancel each other.In deriving (80) from (79) we have used Faraday's law and Gauss' law for magnetism, which of course must be assumed to hold for E(x, t) and B(x, t), no matter how slowly they are varying in space and time.Interestingly, the relation between {P (1) (x, t), M (1) (x, t)} and {P (1) (x, t), M (1) (x, t)} can be understood as another kind of "gauge dependence."Very generally, such sets of fields lead to the same charge-current densities when for a general vector field a(x, t) and a general scalar field b(x, t); here the sets of fields {P (1) (x, t), M (1) (x, t)} and {P (1) (x, t), M (1) (x, t)} of (78) are related by b(x, t) = 0, which can be easily confirmed using (81), under the condition that at t = −∞ the system is unperturbed.
In the past, the origin of the ambiguity arising in the OMP tensor has, like that in the unperturbed macroscopic polarization, been understood in the context of finite sized systems [9,44,45].At a surface where α il cannot be treated as uniform the argument in the above paragraph breaks down, and the surface current that will arise shows an ambiguity that reflects the well-known gauge dependence of α CS [17,18].However, it seems there should be an equivalent bulk interpretation.That certainly holds for the unperturbed macroscopic polarization; in a calculation where the energy eigenstates are chosen and fixed at the start, a "quantum of ambiguity" arises from the gauge dependent term of (55).This term has been shown [3] to depend only on the phase of the determinant of U (k), and as such, is not qualitatively different whether the bands are isolated or not.Indeed, this ambiguity can be understood at the level of a U(1) gauge transformation, which in the simplest of cases takes the form U nα (k) = δ nα e −ik•R ; this corresponds to changing the site with which each Wannier function is associated and, in turn, this changes P (0) by a discrete amount, proportional to R. However, such an interpretation cannot be used to understand the ambiguity associated with the OMP tensor; the gauge dependent term of (76) vanishes in the limit of isolated bands, and thus must interpreted on the more general grounds of a multiband gauge transformation.Yet, at some level these ambiguities appear to be related, since the terms giving rise to them are constructed from the same object, the W a matrix; (55) contains only diagonal matrix elements, while (76) contains only off-diagonal matrix elements.Perhaps it is from this perspective that a bulk interpretation of the discrete ambiguity associated (76) can be formulated.

B. Limiting cases
In this section we explore the magnetoelectric response in the limit of isolated molecules.First we construct the response tensor for a single molecule, and then use that to construct the response tensor of a crystal in the "molecular crystal limit," which we take to be a model where there is a molecule at each lattice site with orbitals that have no common support with orbitals of molecules associated with other lattice sites.Finally, we show that the molecular crystal limit so obtained is in agreement with the appropriate limit of our general expressions (76,77).

A single molecule
As pointed out earlier (Appendix D of Ref. [10]), the response of a molecule to an applied electromagnetic field can be treated via the same approach we have used here for a crystal, by constructing microscopic polarization and magnetization fields from the electronic Green function, with the expectation values of the electric and magnetic dipole moments following from the single-particle density matrix (51,52,53).However, for a molecule (or atom) it is also possible to follow a more common strategy in molecular physics [11,12], where microscopic polarization and magnetization operators are introduced, leading to operators associated with the electric and magnetic dipole moments.We present that approach here (Appendix C of Mahon et al. [10]) to better make the connection between this calculation and molecular physics.
We take the initially unperturbed system to be described by (3,4), where we include an A static (x) -which of course need not be periodic, since we are considering a localized system -to guarantee the breaking of timereversal symmetry, and consider a V (x) that vanishes as |x| → ∞ and that does not satisfy inversion symmetry about any point.The latter condition could arise in a molecule because of a non-centrosymmetric configuration of the nuclei, or even in an atomic system because of an imposed dc electric field that is considered part of the unperturbed Hamiltonian.We consider a "special point" R = 0, which for a molecule could be taken to be, say, the center of mass of the ion cores and for an atom could be taken as the position of the ion core.In the frozen-ion approximation the contribution of the ions to the multipole moments of a molecule will not affect the perturbative response calculation we make, so we neglect them.We introduce a "special point" electron field operator (see Appendix C of Mahon et al. [10]) ψ sp (x, t), ψ sp (x, t) = e −iΦ(x,0;t) ψ(x, t), where Φ(x, 0; t) ≡ e c s a (w; x, 0)A a (w, t)dw, and A(x, t) is again the vector potential describing the applied electromagnetic field.Then the Hamiltonian operator governing the evolution of ψ sp (x, t) is where p(x, 0; t) is given by (22) and Ω 0 0 (x, t) by (C3); we have We begin the evolution at a time t 0 before the electromagnetic field is applied, taking as the (Heisenberg) ket the ground state |G ; at this time we have ψ sp (x, t 0 ) = ψ(x, t 0 ) ≡ ψ(x).The dynamics can equivalently be described in a Schrödinger picture where the field operator is fixed at ψ(x) and the ket evolves from |G at t 0 according to a Hamiltonian operator H(t) having the form of (84), but with ψ sp (x, t) replaced by ψ(x).Using the approximate expressions (30,31) for Ω 0 (x, t) and Ω 0 0 (x, t), and neglecting the variation of the electric and magnetic fields over the atom or molecule, we can write H(t) as where E(t) ≡ E(0, t) and B(t) ≡ B(0, t), and the operator for the electric dipole moment µ, and operators for the paramagnetic ( ν P ) and diamagnetic ( ν D (t)) contributions to the magnetic dipole moment operator ν(t) = ν P + ν D (t) are given by µ = e ψ † (x)x ψ(x)dx, Note that the last expression is the field theoretic analogue of (63).Taking H 0 |G = E G |G , we then calculate the response, to first-order, to applied static fields E and B. The diamagnetic term in (85) will make no contribution to the response to first order, and if we denote the excited states by |H (with energies E H ), then by standard perturbation theory the response of the expectation value of the electric dipole moment operator to the magnetic field, µ atom , and the response of the expectation value of the magnetic dipole moment operator to the electric field, ν (E) atom , are given by where We expand the field operator ψ(x) in terms of orbitals W v (x) that are occupied in the ground state and orbitals W c (x) that are not, where c v and c c are fermionic annihilation operators generating single-particle eigenfunctions of H 0 (x, p(x)) with energies E v and E c respectively, with where and Here we have inserted a complete set of states {|n }, with single-particle energies {E n }, into the first expression, where the label n ranges over all v and c, and in going to the third line we have used the commutation relation of x and H 0 (x, p(x)) to write the matrix element of p b (x) in terms of x b nv in the usual way.From (89) we then have As in a solid, this vanishes unless both time-reversal symmetry and inversion symmetry are broken: For if there is time-reversal symmetry the orbitals W n (x) can be chosen to be real and the quantity inside the brackets is purely imaginary, while if there is inversion symmetry the matrix elements x i vc themselves vanish.Splitting the sum over n in (91) into a sum over filled states v and a sum over empty states c we can write where and In the second of these expressions we sequentially put where both n and n range over all states; noting that the sums over n and n give no net contribution, we have where the quantity in brackets is real.In three dimensions at least two of i, l, a, b must be identical; if i = l the expression is found to vanish, and in general we can write In the presence of uniform applied fields E and B, and in the frozen-ion approximation, each magnetoelectric response tensor component αil of a molecule can then be written (92) as the sum of a Chern-Simons-like term (94) and a cross-gap-like term (93).As in a solid, the Chern-Simons-like contribution depends only on the occupied orbitals, and describes an isotropic response, regardless of how complicated might be the structure of the molecule.Note also that in a situation where all relevant initially occupied (unoccupied) orbitals are degenerate, , the cross-gap-like term vanishes, in line with the vanishing of the cross-gap term in a solid when all relevant initially occupied (unoccupied) bands are degenerate [20].

The molecular crystal limit
We can now construct the "molecular crystal limit," in which we consider a periodic array of molecules where the orbitals associated with a molecule at a given lattice site share no common support with those of molecules associated with other lattice sites; again, we take the external electric and magnetic fields to which the molecules respond to be the macroscopic Maxwell fields.Since the first-order modifications to the moments associated with each molecule are given by (86), using the expressions (54) for P and M , together with the defining equation (1) for the orbital magnetoelectric polarizability tensor, from (92,93,94) we have where simply with αil G and αil CS given by (93,94) respectively.Alternately, rather than building up the molecular crystal limit by assembling a collection of molecules, we can imagine starting with a full bandstructure calculation and taking the limit where the Wannier functions associated with each lattice site have no common support with the Wannier functions associated with a different lattice site.We also take the ELWFs to be eigenfunctions of the unperturbed Hamiltonian, (3), and so E nk → E n .These conditions lead to simplifications in the general expressions (76,77), and when they are employed the result should reproduce the molecular crystal limit (95,96).A first simplification is that, since the bands are flat, ∂ a (E ik + E jk ) → 0. Further, taking the orbitals introduced in (88) to be the ELWFs W v0 (x) and W c0 (x), the flat bands can be identified by taking U nα (k) = δ nα in (5) relating |αR with |ψ nk .Hence the Hermitian matrices W a (12) vanish, and we need not distinguish between the connections (9,10); the lack of common support for orbitals associated with different lattice sites, together with (8), also implies that ξ a cv , etc., are independent of k, and from ( 8) and (9,10) we take where we have used (90) and Applying these molecular crystal conditions to the general expressions (76,77), we indeed recover (96), as expected.Note that in this limit there is no gauge dependence in α il CS , since the bands are isolated, and as well both α il CS and α il G arise solely from dynamical contributions.For while a calculation of the magnetic susceptibility would involve a compositional contribution due to the diamagnetic response, even in the molecular crystal limit (see the discussion around (63)), the contributions to the magnetoelectric response in that limit are purely dynamical; were the calculation for a single molecule here done in terms of the Green function strategy used for a crystal, the response would result from changes in η α0;β0 .Of course, in the molecular crystal limit all responses are what we have called atomic-like rather than itinerant.

C. Microscopic origin of α il
G and α il

CS
While the qualitative features of the two contributions to the OMP tensor have been discussed earlier [20], the microscopic nature of the approach implemented here can provide further insight into the character of the cross-gap and Chern-Simons contributions.
Since both of these contributions are non-vanishing in the molecular crystal limit, neither can simply be understood as entirely a consequence of the delocalized nature of Bloch electrons.For the Chern-Simons contribution, this is in agreement with earlier work [20] where a particular model for a molecule at a lattice site was constructed that exhibits a Chern-Simons-like response; our expression (94) for the Chern-Simons-like response of an arbitrary molecule generalizes that result.However, neither contribution can be understood as a purely "localized molecule-like contribution" either, because the full expressions for α il CS (76) and α il G (77) contain terms that vanish in the molecular crystal limit; this is again in agreement with earlier arguments [20].
When moving from the molecular crystal limit of the Chern-Simons and cross-gap tensors, where the only contributions are atomic-like, to the full crystal expressions, both acquire itinerant contributions.In addition, while the cross-gap tensor is purely dynamical in nature, both in the molecular crystal limit and more generally, the Chern-Simons tensor acquires a compositional modification when moving from the molecular crystal limit to the general expression for a crystal.This suggests that perhaps it is through this compositional modification that a bulk interpretation of the discrete ambiguity associated with the Chern-Simons tensor can be constructed.We plan to explore this conjecture in a future publication.

V. CONCLUSION
We have implemented a previously developed [10] microscopic theory of polarization and magnetization to study the zero temperature, orbital electronic response of a class of insulators to applied uniform dc electric and magnetic fields.To first-order in the applied fields, the free charge and current densities vanish [10]; thus the perturbative modifications to both the electronic charge and current density expectation values due to the applied fields can be found directly from the modifications to the microscopic polarization and magnetization fields.Associated with the dipole moment of the microscopic polarization (magnetization) field is a macroscopic polarization (magnetization), for which various response tensors can be extracted.
A quantity central in any calculation implementing this microscopic theory is the single-particle density matrix.We began by re-expressing the equation governing its dynamics to include an arbitrary lattice site, R a , which is to be used as a reference site for the applied electromagnetic field when calculating modifications to site quantities.This strategy will be useful in future work, where we plan to take into account the spatial variations of the electric and magnetic fields.However, in the calculation reported here we have restricted ourselves to the limit of uniform applied dc electric and magnetic fields.We began by reproducing the usual electric susceptibility [41], and then found the orbital magnetoelectric polarizability tensor.Generally, the OMP tensor is written as a combination of the isotropic, Chern-Simons contribution and the cross-gap contribution; this microscopic theory reproduces the usual result [17][18][19][20].
In the course of the perturbative analysis it became evident that there are generally two distinct types of modifications that contribute to response tensors.The first type arose from modifications to the single-particle density matrix due to the applied fields, and were termed "dynamical."The other type arose from modifications to the diagonal elements of site quantity specific matrices, and were termed "compositional."The electric susceptibility was found to arise from only a dynamical modification.In our analysis of the magnetoelectric response, we found three terms with distinct microscopic origins combine to form the OMP tensor: a dynamical modification to the atomic-like contribution, and dynamical and compositional modifications to the itinerant contribution.We found the Chern-Simons contribution to arise from a combination of parts of the atomic-like dynamical modification, and the itinerant compositional modification.The cross-gap contribution was found to arise from the remainder of the atomic-like dynamical modification, and the itinerant dynamical modification.We have also compared our expressions with those that arise in the molecular crystal limit, a model in which a periodic array of isolated molecules is considered.
As is well known, the first-order macroscopic bulk charge and current densities that result from applied uniform dc electric and magnetic fields are gauge invariant, but at a surface where the bulk arguments underlying this result are not valid, gauge-dependent currents emerge.The ambiguity associated with such surface currents can be studied in detail via the microscopic theory that underpins the approach taken here; we plan to consider this in future work.We also intend to extend the calculations presented in this paper to take into account the frequency dependence of the response tensors in general.This will lead to a description of frequency-dependent magnetoelectric effects, such as optical activity, for which the breaking of both spatial-inversion symmetry and timereversal symmetry in the underlying material system is not required.A microscopic understanding of the mechanisms giving rise to such effects is now accessible via this formalism.
space space where E ≡ E(R a , ω = 0) for any R a , in this limit.Then, recalling (5), we find where we have used the identity as well as ( 8) and (11), in going to the final expression.Notice that the approximation of an applied electric field that varies little on the scale of the lattice constant gives rise to a simplified form of (29), which in turn allows η (E) αR;βR to be written as a single Brillouin zone integral.
To derive the expression for η (B) αR;βR a similar procedure is followed; however, one must also use (15).In this case, the approximation that the magnetic field varies little on the scale of the lattice constant allows η (B) αR;βR to be written as a single Brillouin zone integral.so we have Ω 0 y (x, t) (x i − y i )E i (y, t) + 1 2 (x i − y i )(x k − y k ) ∂E i (y, t) ∂y k = (x i − y i )E i (y, t) ∂E k (y, t) ∂y i = (x − y) • E(y, t) + 1 2 (x i − y i )(x k − y k )F ik (y, t), where we have used (33).Finally, we look at ∆(x, z, y; t).This can be done "by hand" when the magnetic field is uniform, but in what follows we work it out formally.c e ∆(x, z, y; t) A i (y, t) s i (w; x, z) + s i (w; y, x) + s i (w; z, y) dw + ∂A i (y, t) ∂y j (w j − y j ) s i (w; x, z) + s i (w; y, x) + s i (w; z, y) dw + . . .
So, in all (w j − R j ) s i (w; x, z) + s i (w; y, x) + s i (w; z, y) dw = 1 2 (x i − y i )(z j − y j ) − 1 2 (z i − y i )(x j − y j ), and c e ∆(x, z, y; t) ∂A i (y, t) ∂y j (x i − y i )(z j − y j ) − (z i − y i )(x j − y j ) = 1 2 ∂A i (y, t) ∂y j − ∂A j (y, t) ∂y i (x i − y i )(z j − y j ), and since ∂A j (y, t) ∂y i − ∂A i (y, t) ∂y j = kij B k (y, t), We collect all of these approximate expressions in (31,30,32).The expansions of Ω j y (x, t) and Ω 0 y (x, t) derived here can also be derived using a formal expansion of the relators (C1) about u = 0.Here we outline the steps in going from the first to the second equality of Eq. (76).The final term of the first line of (76) can be re-expressed, using (11) where in going from the second to third line we have used the fact that, in three-dimensions, at least two of i, l, a, b must be identical; if i = l the expression is found to vanish.This sort of argument is often used in what follows.The contribution to (76) that is linear in W (notice the penultimate term of (76) also contributes here) is proportional to the Brillouin zone integral of where we have used an integration by parts on the initially linear in W term.This "miraculous cancellation" is also presented in Appendix C of Vanderbilt [9].Thus, the cubic in W contribution is the only gauge dependent term that has not vanished, or been cancelled.Its contribution to (76) is proportional to the Brillouin zone integral of The contribution of this term to (76) is proportional to the Brillouin zone integral of the well-known term arising from the gauge-transformation of the Chern-Simons 3-form (see, e.g., Eq.C.19 of Vanderbilt [9]).The contribution independent of W is proportional (notice the first term of (76) contributes here) to the Brillouin zone integral of which is equivalent to Eqs. (A11a)+(A11b) of Essin et al. [20].Then, in all, (76) results from the combination of (D1)+(D2).
Appendix D: Constructing the Chern-Simons contribution to the OMP tensor

f
ji Re iξ i ji W b in ξ a nj = lab Re nms if n W i nm ξ b ms ξ a sn − ξ i sm ξ a mn W b ns + lab ijn f i Re iξ i ji W b in ξ a nj = δ il lab Re nms if n ξ b ms ξ a sn W l nm + ξ l sm ξ b mn W a ns + ξ a sm ξ l mn W b ns = δ il cab Re nms if n ξ b nm ξ a ms W c sn .Now the combined contribution of the linear and quadratic in W terms is proportional to

lab Re nms if n ξ i nm ξ b ms ξ a sn + lab 2 cvv
Re (∂ i v|c) (c|∂ a v ) (v |∂ b v) + cv Re (∂ i v|c) (∂ a c|∂ b v) = lab Re vm (∂ i v|m) (∂ b m|∂ a v) − (∂ i v|c) (∂ b c|∂ a v) + 2 lab cvv Re (∂ i v|c) (c|∂ a v ) (v |∂ b v) = lab Re vv (∂ i v|v ) (∂ b v |∂ a v) + 2 lab cvv Re (∂ i v|c) (c|∂ a v ) (v |∂ b v) ,(D2) , as proportional to the Brillouin zone integral of (ξ i ps+ W i ps )U sγ ∂ b U † γm (ξ a mn + W a mn )U nα U †where we have used the identities iab ∂ b ξ a mn = i iab s ξ b ms ξ a sn and iab ∂ b W a mn = −i iab s W b ms W a sn ; in the following we will also often use W a mn = 0 only if f m = f n .We now consider the contributions to (76) at each order in W. We first consider the terms quadratic in W; their contribution to (76) is proportional to the Brillouin zone integral of ns = δ il cab Re nms if n ξ a nm W c ms W b sn ,