From magnetoelectric response to optical activity

We apply a microscopic theory of polarization and magnetization to crystalline insulators at zero temperature and consider the orbital electronic contribution of the linear response to spatially varying, time-dependent electromagnetic fields. The charge and current density expectation values generally depend on both the microscopic polarization and magnetization fields, and on the microscopic free charge and current densities. But contributions from the latter vanish in linear response for the class of insulators we consider. Thus we need only consider the former, which can be decomposed into"site"polarization and magnetization fields, from which"site multipole moments"can be constructed. Macroscopic polarization and magnetization fields follow, and we identify the relevant contributions to them; for electromagnetic fields varying little over a lattice constant these are the electric and magnetic dipole moments per unit volume, and the electric quadrupole moment per unit volume. A description of optical activity and related magneto-optical phenomena follows from the response of these macroscopic quantities to the electromagnetic field and, while in this paper we work within the independent particle and frozen-ion approximations, both optical rotary dispersion and circular dichroism can be described with this strategy. Earlier expressions describing the magnetoelectric effect are recovered as the zero frequency limit of our more general equations. Since our site quantities are introduced with the use of Wannier functions, the site multipole moments and their macroscopic analogs are generally gauge dependent. However, the resulting macroscopic charge and current densities, together with the optical effects to which they lead, are gauge invariant, as would be physically expected.


I. INTRODUCTION
In a material that is optically active the plane of polarization of light rotates as the light propagates through the medium; the rotation is associated with a difference in the phase velocities of right-and left-handed circularly polarized light. The frequency dependence of the rotation is called optical rotary dispersion, and the associated difference in absorption of light of the different circular polarizations is called circular dichroism.
The study of optical activity has a long history. Pasteur was the first to associate it with structural dissymmetry [1], and as early as 1928 its first quantum mechanical description was given by Rosenfeld [2]. This phenomenon is most often observed in liquid solutions. The usual solvent, water, is not itself optically active, but the solution is optically active if the symmetry group characterizing the structure of the solute molecules contains no improper rotations. Early theoretical treatments involved models of solute molecules based on at least two coupled oscillators at different sites in each molecule [3], and it was natural to associate optical activity with the variation of the electromagnetic field across the molecule. However, an alternate approach [4] is to consider the electric and magnetic multipole moments of each molecule as a whole, and to describe their response to the electromagnetic field and its derivatives at a nominal center of the molecule. Optical activity is then typically associated with the response of the electric dipole moment to both such as σ il (q, ω) depends on q as well as ω. If timereversal symmetry holds before the medium is subjected to the electromagnetic field, then the optical activity has been called natural [7]. If time-reversal symmetry is broken, then there are generally additional contributions to σ ilj (ω), and as well a rotation of the plane of polarization of light can result from an asymmetric component of σ il (ω), as σ il (ω) = σ li (ω) in general. This latter phenomenon can be thought of as an "internal" Faraday effect.
Yet such a general treatment of linear optical properties of media based on σ il (q, ω) has its drawbacks. First, when using the minimal coupling Hamiltonian and directly calculating the expectation value of the electronic current density operator, "artificial divergences" can arise when the number of bands involved in the calculation is necessarily truncated; sum rules must be employed before such a truncation is performed to avoid these [9,10]. Second, although one can attribute different constituents of σ ilj (ω) to the purported response of different multipole moments [7], those multipole moments, and the physical insight they carry, do not directly arise in the calculation. And third, the bulk relation (1) and its expansion (2) give little direction on how to even approximately treat the subtleties that would arise if one considered a finite system and had to be concerned with effects at interfaces.
A strategy that is more physical is certainly available for crystalline systems in the "molecular crystal limit." In this limit we imagine molecules, here with no improper rotations in their symmetry group, positioned at lattice sites with a lattice constant sufficiently large that electrons can be considered essentially "bound" to one molecule or another, but still much less than the wavelength of light. Adopting the approach of molecular physics [11], multipole moments can be associated with each molecule and from these one can introduce macroscopic fields P i mol (x, t), M i mol (x, t), and Q ij mol (x, t), describing respectively the electric dipole, the magnetic dipole, and the electric quadrupole moment per unit volume of the "molecular crystal." The macroscopic charge and current densities are then given by mol (x, t) = −∇ · P mol (x, t), where the polarization and magnetization fields are given by with ". . ." indicating contributions from higher-order multipole moments. Neglecting local field corrections, from the response tensors associated with the multipole moments of the molecules themselves one can then identify bulk linear response tensorsχ il E (ω),γ ijl (ω),β il P (ω),β il M (ω), andχ ijl Q (ω) that relate the multipole moments to the macroscopic electric and magnetic fields, +γ ijl (ω)F jl (x, ω) +β il P (ω)B l (x, ω) + . . . , where the superscript (0) identifies the contribution to a net quantity from the unperturbed system, ". . ." here indicate contributions that are higher order in the macroscopic electric and magnetic fields and their derivatives, including the linear response of M mol to B, is the symmetrized (spatial) derivative of the macroscopic electric field evaluated at x, and the circle accents identify that these linear response tensors are valid in the molecular crystal limit.
Using (5,4) in (3), transforming to wave-vector space, and comparing with (2), we can construct σ ilj mol (ω) in terms ofγ ijl (ω),β il P (ω),β il M (ω), andχ ijl Q (ω). Such a calculation based on molecular response, done in terms of the multipole Hamiltonian familiar in molecular physics [12], does not suffer from the artificial divergences mentioned above; thus the resulting expression for σ ilj mol (ω) is well behaved. In addition, if time-reversal symmetry is broken before the molecules are subjected to the electromagnetic field, thenχ il E (ω) =χ li E (ω), which gives σ il mol (ω) = σ li mol (ω), leading to another source of the rotation of the plane of polarization of light as it propagates through the molecular crystal. Here the multipole moments of the molecules explicitly appear, and with the underlying macroscopic fields P i mol (x, t), M i mol (x, t), and Q ij mol (x, t) in hand one could begin to consider electrodynamics in the presence of interfaces.
But now what of more realistic models of crystalline materials, wherein the molecular crystal limit is not satisfied? Although there are no centers with which particular electrons are associated, in the "modern theory of polarization and magnetization" one can still define electric and magnetic dipole moments [13][14][15], albeit indirectly, through the response of the electronic charge and current densities to external electromagnetic fields. This approach is generally focused on the limit of static applied fields to insulators, the inclusion of higher-order moments in this framework is work in progress [16,17], and its generalization to optical fields is not obvious.
We recently introduced [18,19] a general approach to calculating both the static and the optical perturbative response of a medium based on the introduction of microscopic polarization and magnetization fields, p(x, t) and m(x, t). The usual macroscopic fields P (x, t) and M (x, t) are defined as the spatial averages of the corresponding microscopic fields. In general there are also microscopic free charge and current densities, the spatial averages of which are identified as the macroscopic free charge and current densities. However, at zero temperature, for the class of insulating crystals to which we restrict ourselves in this paper -which includes ordinary insulators [20] and Z 2 topological insulators -those free charge and current densities vanish in linear response, and the full microscopic response, to first order in the electromagnetic field, is captured by p(x, t) and m(x, t).
With the introduction of Wannier functions, the microscopic fields p(x, t) and m(x, t) can be decomposed into constituents associated with each lattice site, and these "site" contributions can be expanded in terms of a series of "site multipole moments." The spatial average of these microscopic fields then leads to an expansion of the macroscopic polarization and magnetization fields in the form (4), even if the molecular crystal limit does not hold. Further, the response of the multipole moments associated with each lattice site can be calculated in terms of the electromagnetic field and its derivatives evaluated at that site, and leads naturally to a description of the linear response that follows the form (5), again even though the molecular crystal limit does not hold. As well, the artificial divergences that can plague standard minimal coupling calculations are absent.
In this approach the site contributions to the electronic component of the microscopic polarization and magnetization fields, and thus to their multipole moments, depend on a modified form of the Wannier functions resulting from a generalized Peierls substitution [18]. There is also a well-known "gauge freedom" in choosing the original Wannier functions from which the modified functions are constructed, for they can be altered by adjusting the k-dependent unitary transformation relating them to the Bloch energy eigenstates. In general this leads to a "gauge dependence" of the site multipole moments, both initially and in their response to the electromagnetic field. And while exponentially localized Wannier functions (ELWFs) would of course be a natural choice for the original Wannier functions, we show that whatever choice is made the resulting electronic charge and current densities predicted are gauge invariant [21], as would be physically expected; the expressions we extract for σ il (ω) and σ ilj (ω) are thus gauge invariant. Even within the independent particle and frozen-ion approximations, which we adopt in this work, we believe this is the first derivation of σ ilj (ω) for an insulator that is valid not only at frequencies below the band gap, but also at frequencies above the band gap where absorption can occur. Thus, the σ ilj (ω) we present provides a description for both the optical rotary dispersion and the circular dichroism of crystalline insulators. At frequencies below the band gap we find agreement with an earlier calculation of σ ilj (ω) [7] that focused on that limit.
The special case of static and uniform electric and mag-netic fields is particularly interesting. In that limit the tensor describing the modification of the polarization due to the electric field becomes symmetric, even in the absence of time-reversal symmetry in the unperturbed crystal. But in the absence of both time-reversal and spatial inversion symmetry, a magnetic field can still induce a polarization and an electric field can still induce a magnetization. This phenomenon is called the magnetoelectric effect [22]; in an earlier work [19] we used our approach to derive the so-called orbital magnetoelectric polarizability (OMP) tensor that describes the magnetoelectric effect in the limit of fixed ion cores and with the neglect of spin contributions, and found agreement with earlier studies based on the "modern theory of polarization and magnetization" [23,24]. Optical activity can be understood as arising from the generalization of the magnetoelectric effect to finite frequencies, where the electromagnetic field is necessarily not uniform; time-reversal symmetry then need not be broken for the phenomenon to occur. And as our calculation is based on a microscopic identification of polarization and magnetization fields, we can identify a finite frequency generalization of the Chern-Simons contribution to the OMP tensor; this contribution is isotropic and thus does not lead to an induced electronic charge-current density in the bulk, which makes it inaccessible to approaches based on the bulk charge-current density response alone. Finally, since our calculation is based on the identification of site quantities, we can easily compare the general response of a crystal to that of a crystal in the molecular crystal limit mentioned above. In this paper we identify expressions for the response tensors χ il E (ω), γ ijl (ω), β il P (ω), β il M (ω), and χ ijl Q (ω) in both cases, indicating the response tensors that are valid in the molecular crystal limit by a circle accent as we have above. In particular, while the OMP tensor is identified with β il P (0) = β li M (0), the relationβ il P (ω) =β li M (−ω) continues to hold for finite frequencies in the molecular crystal limit, but it fails for a crystal more generally. Thus our approach is well positioned to explore the boundary between molecular physics and condensed matter physics in their descriptions of optical activity, and indeed of other optical phenomena.
The structure of this paper is as follows. In Section II we present the basic expressions for the microscopic polarization and magnetization fields, identify the site multipole moments, and present their relation to the macroscopic response functions; some of the details are relegated to Appendices A and B. The linear response of a crystalline insulator, within the independent particle approximation, is presented in Section III. Here for simplicity we neglect the spin degree of freedom and treat the ion cores as fixed. The response of the site multipole moments is detailed in Section IV, where we also consider some of the symmetries of the response tensors. In Section V we construct the linear response of the macroscopic charge and current densities, and identify σ il (ω) and σ ilj (ω); their constituent tensors are listed in Ap-pendix C, and in Appendices D and E we confirm that the response is gauge invariant and thus that σ il (ω) and σ ilj (ω) are as well. We also consider the special case of frequencies below the band gap of the insulator and confirm, using a result presented in Appendix F, that we have agreement with earlier work for σ ilj (ω) [7]. In Section VI we consider the molecular crystal limit and show that in this limit our general crystalline expressions reduce to what would be expected. We discuss and conclude in Section VII.

II. MULTIPOLE MOMENTS
In earlier work [18,19] we showed how the (total) microscopic charge and current densities can be written as where, in this work, with ρ ion (x) the charge density associated with fixed ion cores, and ρ(x, t) and ĵ (x, t) the expectation values of the microscopic electronic charge and current density operators, respectively [25]. These operators are obtained from the minimal coupling Hamiltonian via Noether's theorem and involve the electron field operators and their adjoint, which we take to be the dynamical degrees of freedom of the crystalline system; they evolve under the minimal coupling Hamiltonian, which results in the (assumed classical) electromagnetic field entering (6), and thus in both p(x, t) and m(x, t) generally having a nontrivial dependence on time [18,19]. These microscopic fields can generally be decomposed as a sum of constituent fields [18], one associated with each Bravais lattice vector R characterizing the structure of the unperturbed crystalline system, Each "site" polarization p R (x, t) is related to a portion ρ R (x, t) of the (total) charge density that is associated with the lattice site R, and each "site" magnetization m R (x, t) is related to a portion j R (x, t) +j R (x, t) of the electronic current density that is associated with the lattice site R, where the "relators" s i (x; y, R) and α ib (x; y, R) have been introduced and discussed previously [18]; they are presented in Appendix A. In general the microscopic "free" charge and current densities, ρ F (x, t) and j F (x, t), are also relevant. However, in this paper we assume the crystal to be in its zero temperature ground state before the electromagnetic field is applied and so, for the class of insulators considered here and specified below, both the unperturbed free charge and current densities, and their linear response to the electric and magnetic fields vanish [18]. This is as would be expected physically, and we can henceforth neglect those fields. The macroscopic polarization and magnetization fields, P (x, t) and M (x, t), can be identified as spatial averages of the microscopic fields (7), as discussed in Appendix B. Anticipating the integration over each site contribution (8) associated with such spatial averaging, we perform a formal expansion of each site contribution in terms of Dirac δ functions and their derivatives about that site, as we detail in Appendix A. The expansions are characterized by their dependence on a parameter u, and explicitly retaining the terms that are at most linear in that parameter we find where is the electric dipole moment, is the electric quadrupole moment, and is the magnetic dipole moment, each associated with lattice site R; here iab is the Levi-Civita symbol. Terms that are higher order in u, indicated by ". . ." in the expansions (9), involve the electric octupole moment, the magnetic quadrupole moment, and higher-order moments. For the sort of systems considered here, within the independent particle approximation one can physically expect the response of the moments µ i R (t), q ij R (t), and ν i R (t) to the microscopic electric and magnetic fields to depend on those fields in the neighborhood of R. The approximation of neglecting "local field corrections", which we adopt here, involves taking those fields to simply be the macroscopic fields E(x, t) and B(x, t) that are the spatial averages of the microscopic electric and magnetic fields; we call these macroscopic fields the "Maxwell fields" (see Appendix B). With this approximation, we show in Section IV that the linear response of each site moment (10,11,12) can be related to the Maxwell fields evaluated at that site, E(R, t) and B(R, t), and their spatial derivatives there. Then, implementing the usual Fourier series analysis, we find that the relevant terms are We have chosen to introduce a unit cell volume Ω uc here because, with the neglect of local field corrections, the response tensors χ il E (ω), γ ijl (ω), β il P (ω), β il M (ω), and χ ijl Q (ω) appearing here reduce to those of (5), in the molecular crystal limit. We show in Appendix B that macroscopic multipole moments, analogous to those appearing in (5), can be constructed from the corresponding site multipole moments (14) (see (B5) and (B8)), such that where the unperturbed contributions simply acquire a factor as they are in fact independent of R. Further, the macroscopic charge and current densities are given by with macroscopic polarization and magnetization fields even far from the molecular crystal limit. Notably in the systems we consider here, the unperturbed contributions to (15) vanish when implemented in (16). Hence, the lowest-order charge and current densities arise due to the linear response tensors χ il E (ω), γ ijl (ω), β il P (ω), β il M (ω), and χ ijl Q (ω). In the next two sections we turn to the calculation of these response tensors.

III. LINEAR RESPONSE
The charge and current densities associated with each lattice site that were mentioned above can be written as where ρ ion R (x) is the static contribution to the charge density associated with lattice site R due to the appropriate ion core(s) and where the ρ βR ;αR (x, R; t), j βR ;αR (x, R; t), andj βR ;αR (x, R; t) are generalized (electronic) "site quantity matrix elements" that have been presented earlier [18]. These quantities can be reasonably expected to vanish unless x is "close" to R, guaranteeing that ρ R (x, t), j R (x, t), andj R (x, t) have that property as well. To avoid possible confusion we note that the total microscopic charge and current densities are given by ρ( depends on the net charge density that is associated with R, while m R (x, t) is sensitive to only a portion of the current density that is associated with R.
It is clear from (18) that the single-particle density matrix, η αR ;βR (t), is central in the identification of electronic "site" quantities and in describing their dynamics [18]. This object captures the electronic transition amplitude from a particular Wannier orbital of type α associated with lattice site R to a Wannier orbital of type β associated with R , at time t (see Eq. (33,36) of Ref. [18]).

A. Dynamical and compositional contributions to the multipole moments
The site quantities of primary interest are the Cartesian components of the lowest-order multipole moments (10,11,12) that are associated with lattice site R. Indicating such a site quantity generally by Λ R (t), it is clear that upon inserting the relevant term(s) (18) in the desired site multipole moment expression (10,11,12), is generally of the form where Λ βR ;αR (R; t) is a general (electronic) site quantity matrix element and Λ ion R involves ρ ion R (x). In addition to the dependence of the single-particle density matrix on time, which would be expected in the presence of a time-dependent electromagnetic field, the site quantity matrix elements appearing in (18) also have a time dependence -and thus so do the Λ βR ;αR (R; t) associated with the various site multipole moments -because they themselves depend on the electromagnetic field. This sort of dependence is not unexpected in the response of systems to the full electromagnetic field. The diamagnetic response of an atom, for example, is not due to a change in its wave function when a magnetic field is applied, which would be captured by the single-particle density matrix, but rather arises because the expression of the charge velocity in terms of the canonical momentum is modified.
We begin by expanding all objects in powers of the electromagnetic field, such that etc. Again, the superscript (0) denotes the contribution to the quantity that is independent of the Maxwell fields; this is the value the object would take in the unperturbed system. The superscript (1) denotes the linear response of the quantity to the Maxwell fields [26]. Here ". . ." represent terms that are higher than first order in the Maxwell fields and will later be neglected. Also, for n = 0, ρ ion(n) R (x) = 0 and consequently Λ ion(n) R = 0 as the ion cores are assumed fixed; thus, in describing the electronic response, the net response of the system is captured. From (19) it is clear that there are two (electronic) contributions to the linear response of a general site quantity to the Maxwell fields, We have called [19] the first term on the right-hand side, a "dynamical" contribution to the linear response, because it arises from modifications to the unperturbed single-particle density matrix due to the Maxwell fields, and the other term, a "compositional" contribution, because it arises due to the way in which the site quantity matrix elements themselves depend on the Maxwell fields. As we will show, (22) only describes first-order modifications of single-site properties as a result of the electromagnetic field. Moreover, we will generally decompose the linear response of a site quantity (20) as a sum of the contributions from the Maxwell electric field, its symmeterized derivative, the Maxwell magnetic field, and higher-order derivatives of these fields, such that In general each of the constituents on the right-hand side of (23) is composed of a dynamical contribution and a compositional contribution; for instance, However, for each site multipole moment that is considered only a limited number of the constituents in (23) are retained; this is detailed in Section IV.
In the remainder of this section we determine the evolution of η (1) αR ;βR (t) from the initial η (0) αR ;βR , and in the following section we combine those results with the Λ (0) βR ;αR (R, t) and Λ (1) βR ;αR (R, t) appropriately, to find the linear response of the site multipole moments.
B. Evolution of the single-particle density matrix In the independent particle approximation, the equations of motion governing the evolution of the (electronic) single-particle density matrix elements take the form [18] i where The quantitiesH µR1;νR2 (R a , t) can be understood as generalized "hopping" matrix elements and are as previously defined [19]. With the neglect of local field corrections they involve the Maxwell fields in the neighborhood of the lattice sites appearing, including lattice site R a . This lattice site can be arbitrarily chosen [19], and we discuss its choice below. The Maxwell field B(x, t) also enters in the quantities ∆(R , . . . , R ; t), which are proportional to the magnetic flux through the surface generated by connecting the points (R , . . . , R ) with straight lines, when the usual choice of straight-line paths for the relators is adopted (see Appendix A). In this work, this choice is always made.
An expansion of the hopping matrix elements in powers of the electromagnetic field [19] gives where W αR (x) ≡ x|αR is the ELWF identified by its type α and the lattice site R with which it is associated, and H 0 x, p(x) is the differential operator that governs the dynamics of the electron field operators in the unperturbed infinite crystal; we take where we allow for a static and periodic magnetic field described by a vector potential satisfying A static (x) = A static (x + R) for any lattice vector R [27]. The eigen- x|nk being a cell-periodic function, and are identified by a band index n and an index k identifying the associated crystal momentum k; we denote the corresponding eigenvalues by E nk . These energy eigenfunctions can be used to construct ELWFs [28][29][30][31][32] via where the vectors |αk are related to the vectors |nk by a (unitary) "multiband gauge transformation", Generally for an insulating crystal in its zero temperature ground state there is a filling factor f n associated with each |nk that is either 0 or 1. And in this paper we restrict ourselves to the class of insulators characterized by the property that the sets of occupied and unoccupied cell-periodic functions x|nk can be used separately to construct sets of ELWFs; this class contains both ordinary insulators and Z 2 topological insulators [29,32]. Thus we can associate an analogous filling factor f α with each |αk that is also either 0 or 1 depending on the occupancy of the |nk used in the construction of that particular |αk , and so U nα (k) = 0 only if f n = f α . That is, (27) is a unitary transformation between elements of the (un)occupied subspace of the electronic Hilbert space alone. Associated with the set of vectors {|nk } is a non-Abelian Berry connection, and with the set of vectors {|αk } is another, where we have defined the Hermitian matrix W a (k), populated by elements and in general we adopt the shorthand ∂ a ≡ ∂/∂k a . For the class of insulators we consider, W a mn (k) is nonzero only if f m = f n . In what follows, the k dependence of the above introduced objects is usually kept implicit.
At zeroth order in the Maxwell fields, the unperturbed expression (25) can be implemented into (24); with this, the elements of the unperturbed single-particle density matrix for the zero temperature ground state are identified to be Following the same procedure, but now retaining only the terms that are first order in the Maxwell fields, the linear response of the single-particle density matrix is identified as [19] η (1) (recall (13) Here H Ra (x, ω) involves the electromagnetic field via the scalar quantity Ω 0 Ra (x, ω) and the vector quantity Ω Ra (x, ω). Very generally, Ω 0 y (x, ω) involves a line in-tegral involving E(z, ω) from y to x, while Ω y (x, ω) involves a more complicated line integral involving B(z, ω) from y to x [18], which also appears in (31). For |x − y| on the order of a lattice constant, an expansion [19] gives and similarly for |x − y| and |x − z| on the order of a lattice constant we have (see Appendix A). With the approximations (33,34) we find we can write [19] H (1) Thus R a acts as a natural point about which to expand the electromagnetic field, and a natural choice of R a for use in (31) would be a site "close" to R or R . Still leaving that choice open, we implement (36) in (31) to identify the contributions to η (1) αR ;βR (ω) from the electric field, the symmetrized derivative of the electric field, and the magnetic field. We write this decomposition as which will allow for the identification of the dynamical contributions to the constituents of (23). Notably we will neglect contributions related to the spatial variation of the magnetic field, since we identify any such terms as higher-order modifications (see Appendix A).
C. Linear response of the single-particle density matrix By implementing the first and second terms of (36) in (31) via (32), and noting (35) is independent of E l (x, ω) and F jl (x, ω), the linear response of the single-particle density matrix to the Maxwell electric field and its symmetrized derivative are found to be and where we have introduced The linear response to the Maxwell magnetic field involves the third term of (36), and using (35) it is found to be where we have introduced In the limit of uniform dc Maxwell fields, both (39) and the final term of (41) vanish trivially, and B ab mn (k, ω = 0) reduces to the previously defined B ab mn (k); the above expressions are thus consistent with past results [19]. Also, (38,39,41) are written as single Brillouin zone integrals. In past work [19] we showed explicitly how this reduction to a single k-integral arises when implementing (31) to find the ω = 0 component of (38). However, this reduction emerges more generally as a consequence of expressing the variation of the electromagnetic field over the unit cell through the expansion, following from (33,34,35), in powers of the length of the unit cell divided by the wavelength of light. Upon implementing the resulting expressions in (31), via (32,36), and using previously introduced identities [19], the reduction to a single k-integral occurs.

IV. THE RESPONSE TENSORS
The linear response of the single-particle density matrix (37) allows for the identification of the dynamical contributions (21) to the linear response of the site multipole moments (10,11,12) to the Maxwell fields and their derivatives. For the compositional contributions (22), we implement (30) and (13) to write which can also be decomposed into contributions due to the Maxwell fields and their derivatives. The decomposition of the net linear response is given by (23). However, for a given site multipole moment of interest, µ i R (ω), , we consider only those constituents of the associated Λ (1) αR ;αR (R; ω) and of η (1) αR ;βR (ω) that lead to the explicitly included first-order terms in (14); these are µ We have justified the retention of only these terms in Appendix A. From (10,11,12) follow the relevant site quantity matrix elements associated with the site multipole moments of interest in terms of the site quantity matrix elements appearing in (18). These have been presented earlier [18], and we now use them to determine the desired response tensors.
We are finally in a position to set R a . When considering the dynamical contribution to the linear response of a particular multipole moment associated with lattice site R to a particular Maxwell field or its derivative, we always choose R a = R in the constituent of η (1) αR ;βR (ω) being implemented. For instance, when considering µ αR ;βR (ω). In the expressions that follow, one of the matrix element indices R or R always equals R, making the use of the expansions (33,34,35) in deriving (38,39,41) sensible. As well, we are always able to manipulate the expressions for the Λ (1) αR ;αR (R; ω) that appear in such a way that the Maxwell fields are evaluated at R. Collectively, this results in the net linear response of the moments associated with R being related to the electric field, the magnetic field, and the symmetrized derivative of the electric field evaluated at R, and facilitates the passage to a relation between the macroscopic polarization and magnetization and the Maxwell fields (see (14,15) and Appendix B).
A. Linear response of the electric moments

Dipole response to the electric field
We begin with the linear response of a site electric dipole moment to the Maxwell electric field, µ (E) R (ω). While ρ (1) αR ;αR (x, R; t) depends on the magnetic field, it does not depend on the electric field; the compositional contribution µ i(E;II) R (ω) vanishes. The linear response of this quantity to the electric field is thus entirely dynamical -it is solely due to η (E) αR ;βR (ω) -and is given by From (14) we identify which is gauge invariant in that it is independent of U nα (k). In general (44) is not symmetric under exchange of Cartesian components i and l. However, if the unperturbed crystal is time-reversal symmetric, then one can show χ il E (ω) is equal to χ li E (ω); in the ω → 0 limit the exchange of these indices is always symmetric, and if absorption is neglected, then χ il E (ω) is equal to χ li E (−ω) [33]. To obtain (43) we have implemented previously introduced [19] identities, and in the remainder of this section we often do so; Eq. (8,14,15) of that work are particularly relevant.
We now take into account the spatial variation of the Maxwell electric field. The compositional contribution vanishes, as in the response calculated above; the linear response is entirely dynamical, and it is given by The two distinct contributions appearing in the braces of (45) originate individually from the first and second lines of (39), while the contribution from final line of (39) vanishes. Via (14) we again identify the relevant response tensor. We explicitly symmeterize the indices labeling Cartesian components that are contracted with the symmeterized derivative of the Maxwell electric field, and we find Notably γ ijl (ω) = γ ilj (ω); the underlying presence of this symmetry -even if the indices j and l were not contracted with an object symmetric in those indices, F jl (x, ω), in (45) -can be recognized by identifying that the objects carrying these indices in (46) originate from the second term of (36), which was used in (31) to obtain (39). There j and l are clearly symmetric as the components of x and R a commute. Unlike χ il E (ω), γ ijl (ω) is gauge dependent.

Dipole response to the magnetic field
As ρ αR ;αR (x, R; t) does depend on the magnetic field, there are nonvanishing compositional and dynamical contributions to µ βR ;αR (x, R; ω) be the part of ρ (1) βR ;αR (x, R; ω) that is proportional to the magnetic field, the compositional contribution is given by Note that, in going from the first to the final equality, we ensure (using Eq. (29) of Ref. [18]) the ∆(R 1 , y, R; ω) that enters via ρ (B) αR ;αR (y, R; ω) (see Eq. (28,45) of Ref. [18]) is in a form that, upon implementing (35), the magnetic field is evaluated at R. Writing β il(II) P as the compositional contribution to β il P (ω) (see (15)), we have which is again gauge dependent. Interestingly, β il(II) P is independent of frequency and is identical to the compositional contribution to the tensor describing the linear response of the electric dipole moment to a uniform dc magnetic field (see Ref. [19]).
The form of the dynamical contribution is similar to (43) but with η αR ;βR (ω); denoting its contribution to β il P (ω) by β il(I) which is also gauge dependent. The two distinct terms appearing in the braces of (49) originate individually from the first two lines of (41), and the contribution from final line of (41) vanishes. We separate out the frequencyindependent terms that appear in (49) and combine them with (48). Together these terms give rise to the previously found [19,23,24] OMP tensor, α il = α il G + δ il α CS , where α CS is termed the Chern-Simons contribution and α il G the cross-gap contribution; the expressions for these are given in Appendix C. The remaining terms are used in the construction of an explicitly frequency-dependent response tensor, α il P (ω), which vanishes in the ω → 0 limit. In all, then, we find where we have defined Notably α il P (ω) arises due to the linear response of the single-particle density matrix η αR ;βR (ω) alone, making it entirely the result of a dynamical contribution. Here α CS and α il P (ω) are gauge dependent, but α il G is not.

Quadrupole response to the electric field
The compositional contribution to q ij(E) R (ω) vanishes as ρ (1) αR ;αR (x, R; t) does not depend on the electric field. The dynamical contribution involves η (E) αR ;βR (ω) and again takes the form of (43), except that it will be the second moment (see (11)) of ρ (0) βR ;αR (y, R) that will appear rather than the first. Using the expression for q ij(E) R (ω) that results, from the second of (14) we identify another gauge-dependent response tensor, with χ ijl Q (ω) = χ jil Q (ω). This symmetry of the response tensor is a consequence of the symmetry in the definition (11) of q ij R (t). Notably, both χ ijl Q (ω) and γ ijl (ω) arise from dynamical contributions alone, and are of similar form apart from an energy derivative term that appears in γ ijl (ω).

B. Linear response of the magnetic dipole moment to the electric field
The expression (12) for a site magnetic dipole moment shows that there are two contributions; an "atomic-like" contribution arising due to j R (y, t), and an "itinerant" contribution arising due toj R (y, t) [18]. We denote the contribution of the first of these to the linear response of the site magnetic dipole moment to the Maxwell electric field ν    (14) We now identify these contributions.

Response of the atomic-like contribution
As j αR ;αR (y, R; t) does not depend on the electric field,ν R (ω), we compare to (14) and extract

Response of the itinerant contribution
In contrast, sincej αR ;αR (y, R; t) does depend on the electric field, there will be a nonvanishing compositional We find the compositional contribution to bẽ which, like (48), does not depend on frequency. To ensure that the electric field is evaluated at R, in reaching (56) we have used the form of F µR1;νR2 αR ;βR (t) presented above in the expression forj αR ;αR (y, R; t) (Eq. (60,61,62) of Ref. [18]), and set R a = R. The dynamical contribution is β il(I) While (55) and (56) are generally gauge dependent, (57) is only gauge dependent if there are degeneracies present in the unperturbed system. Very generally, there is a simplification that occurs when (55,56,57) are summed to form the total response tensor (54); the gauge-dependent terms appearing in (57) cancel with terms appearing in (55), and as a result the gauge-dependent terms appearing in the total β il M (ω) do not explicitly depend on the energies E nk . In all we have where we have separated out the dc-like terms, α li = α li G + δ il α CS , as in (50), and defined Like α il P (ω), α li M (ω) is entirely a consequence of a dynamical contribution. The form of (59) is similar to that (51) found for α il P (ω), apart from a term related to an energy derivative. Also, like α il G , α il P (ω) and α li M (ω) are "cross-gap" contributions; that is, they depend on both initially occupied and unoccupied Bloch energy eigenstates, and their corresponding energies. Unlike α il G , however, both α il P (ω) and α li M (ω) are gauge dependent.
A qualitative feature shared by the response tensors γ ijl (ω), α il P (ω), α li M (ω), and χ ijl Q (ω) is that they are all gauge dependent. Moreover, the explicitly gaugedependent terms within these tensors are of a similar form; the terms that involve the objects W a nm are all linear in W a nm , and also involve the energies E nk and the non-Abelian Berry connection ξ b nm . This is in contrast to what is found at the level of uniform and static Maxwell fields, where the only gauge dependence of such a tensor enters via the Chern-Simons contribution (C3) to the OMP tensor [19,23,24]. There the explicitly gaugedependent term of α CS involves the W a alone and gives rise to a discrete ambiguity associated with the OMP tensor.

V. MACROSCOPIC CHARGE AND CURRENT DENSITIES
We now construct expressions for the linear response of the macroscopic charge and current densities to the Maxwell fields, and as well identify the effective conductivity tensor σ il (q, ω) to first order in q.

A. The macroscopic current density
Retaining only the contributions to the multipole moments that are linearly induced by the Maxwell fields and that are explicitly included in (15), implementing them into the expressions (16,17) to obtain the linear response of the current density and, following (13), writing this as where α il = α il G +δ il α CS . Of the response tensors appearing here, only χ il E (ω) and α il G are gauge invariant. The rest, which are α CS , α il P (ω), α il M (ω), γ ijl (ω), and χ ijl Q (ω), are all gauge dependent. Yet the linear response of the current density J (1) (x, ω) is in fact gauge invariant. To see this, first note that α CS appears in (60) in the form vanishing via Faraday's law. So in considering the bulk response (60) we can discard α CS , replacing α il by α il G . For the other gauge-dependent terms, we re-express each response tensor as a sum of a gauge-invariant contribution, denoted by a breve accent, and a gauge-dependent contribution. We then find (see Appendix D); that is, the sum of the gaugedependent contributions vanishes. Thus the linear response of the current density is gauge invariant, as expected.

B. The macroscopic charge density
A similar analysis holds for the linear response of the charge density to the Maxwell fields, where from (16) we have (1) Again, retaining only the contributions to P (x, t) that are explicitly included in (17), those involving the electric dipole and quadrupole moments, and retaining only the contributions to the electric dipole and quadrupole moments that are linearly induced by the Maxwell fields and that are explicitly included in (15), for the frequency components we have (1) Again the Chern-Simons coefficient α CS makes no contribution, since it appears in the form vanishing since the Maxwell magnetic field necessarily satisfies ∇·B(x, t) = 0. This is analogous to the scenario for J (1) (x, ω). As was the situation there, we expect (62) to be gauge invariant as a whole. Separating out the explicitly gauge-dependent terms as before, we find (see Appendix E), which is gauge invariant, as expected. We note that the expressions (61,63) satisfy continuity as also expected.

C. The effective conductivity tensor
Finally, we can identify the linear dependence on q of the effective conductivity tensor σ il (q, ω). Fourier transforming (2) to position space, we have Comparing with (61) we can identify which agrees with the usual optical conductivity tensor found via the Kubo formula in the long wavelength limit [18]. Then, defining the dc limit of σ ilj (ω) as and implementing Faraday's law, we can identify All of (65,66,67) are gauge invariant, as expected.
In the absence of time-reversal symmetry, the σ il (ω) of (65) is nonsymmetric and can lead to the rotation of the plane of polarization of light as it propagates through the medium; this can be thought of as an "internal" Faraday effect, as illustrated by the discussion of the molecular crystal limit in the next section. The σ ilj (ω) of (67) is generally nonvanishing and nonsymmetric with respect to the exchange of any of its indices, even in the presence of time-reversal symmetry. But if that symmetry is present, then σ ilj DC will vanish and the resulting σ ilj (ω) describes what has been called natural optical activity [7]. In general the tensor σ ilj (ω) can be evaluated at frequencies above the band gap, and thus can be used to describe both optical rotary dispersion and circular dichroism. Earlier work [7] considered σ ilj (ω) at frequencies below the band gap, where E ck − E vk = ω for all c, v, and k; here (c) v are the band indices labeling Bloch energy eigenstates of the unperturbed Hamiltonian that are initially (un)occupied. To compare our results with theirs, in our expression (67) for σ ilj (ω) we can take the 0 + limit immediately without introducing any divergences, and we follow them [7] in adopting the notation " . =" to identify equalities that only formally hold in this limit. Introducing the shorthand E cvk ≡ E ck − E vk , and putting σ ilj (ω) = Re[σ ilj (ω)] + iIm[σ ilj (ω)], we find and where we have adopted the previously introduced [7] B ab nm ≡ − (see Appendix F). This is in agreement with the orbital electronic contribution to σ ilj (ω) found by Malashevich and Souza [7], as expected. Notably the only nonvanishing contribution to σ ilj (ω) in the ω → 0 limit is due to σ ilj DC , which is purely imaginary, as α il G is real. Thus, in this limit, (68) is expected to vanish, which it does.

VI. THE MOLECULAR CRYSTAL LIMIT
We now consider our response tensors χ il E (ω), γ ijl (ω), , and α li M (ω) in the molecular crystal limit. That is, 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, neglecting any local field corrections. We denote the response tensors in this limit by a circle accent.
We discussed the approach to this limit from the full crystalline expressions earlier [19]; in essence, this limit can be reached by taking the ELWFs (26) to be eigenfunctions of H 0 x, p(x) , in addition to the condition on the common support of these functions mentioned above [34]. The former condition can be achieved by taking E nk → E n and U nα (k) → δ nα , and, consequently, Again restricting ourselves to frequencies below the band gap, as in the second part of Section V C, and implementing these substitutions, (44) becomes where E n n ≡ E n − E n . In the presence of time-reversal symmetry this tensor is symmetric under the exchange of Cartesian components i and l, but in general it is not, and we have onlyχ il E (ω) .
=χ li E (−ω). These results follow the pattern of the corresponding tensor for the more general crystalline system (see text surrounding Eq. (44)). Note that even were it the only response tensor present, an asymmetric χ il E (ω) would be sufficient to lead to the rotation of the polarization of light as it propagates through a medium, as can be easily confirmed. In this molecular crystal limit it is easy to give an example of how this might arise. Suppose, for example, that the breaking of time-reversal necessary for the asymmetricχ il E (ω) occurs because each molecule -or, simpler, atom -is subject to a dc magnetic field that is incorporated in the unperturbed atomic Hamiltonian. Then, if light is propagating in the direction of the dc magnetic field, then the rotation of its plane of polarization that results is just the Faraday effect, which is well known in atomic systems and indeed has a variety of applications [35]. Next, (46) simplifies to and (53) to Recall from previous work [19] α il = e 2 2mcΩ uc lab vcn where, in this limit, Further, (51) becomes Then, combining this with the dc-like contribution, the full response of the polarization to the magnetic field is given bẙ Finally, (59) simplifies to Combining this with the dc-like contribution, the full response of the magnetization to the electric field is given bẙ Physically one expects that an equivalent way to derive these expressions would be to solve for the linearly induced moments of the individual molecules; since local field corrections are being neglected, the fields to which they respond are the Maxwell fields, and the limiting response tensors above should be equal to the appropriate molecular response tensors multiplied by the number of molecules per unit volume, here equal to Ω −1 uc . The molecular calculations can be made with the usual multipole moment Hamiltonian [11,12], which including the moments relevant here can be written aŝ whereĤ 0 mol is the Hamiltonian in the absence of any Maxwell fields; E i (t), F ij (t), and B i (t) are the Cartesian components of the electric field, its symmeterized derivative, and the magnetic field evaluated at the position of the molecule; andμ i ,q ij ,ν i P , andν i D are the indicated components of the operators for the electric dipole and quadrupole moments, and the paramagnetic and diamagnetic dipole moments of the molecule. The diagmagnetic dipole moment is neglected here since it is not involved in optical activity, but the matrix elements of the other moments can be written in terms of the "position" and "momentum" matrix elements (71,75) involving the {W v0 (x)} and the {W c0 (x)}, now identified with the filled and empty orbitals of a molecule fixed at the origin.
The result is that (72,73,74,76,77) are indeed the appropriate molecular response tensors divided by Ω uc . The molecular calculation also clarifies certain symmetries in the expressions in the molecular crystal limit. For example, in this case one can immediately identifẙ . =γ lij (−ω), and the equivalent expressions . =α il M (−ω). The first of these holds because the response calculations leading to both quantities involve the different-time commutator of the electric dipole and the electric quadrupole moment operators, while the second holds because the response calculations leading to both involve the differenttime commutator of the electric dipole and the paramagnetic dipole moment operators. These symmetries no longer hold in the full crystal calculation, where the site multipole moments are not the result of the expectation values of site operators, but rather are evaluated in terms of the single-particle Green function.

VII. CONCLUSION
We have presented a theory for the effective conductivity tensor σ il (q, ω) of a class of insulating crystalline solids at zero temperature. In retaining terms that are at most linear in q, we extract tensors σ il (ω) and σ ilj (ω) that describe phenomena involving the rotation of the plane of polarization of light as it propagates through a medium; the former contributes through its antisymmetric part only when time-reversal symmetry is broken in the unperturbed system, and can be considered as describing an "internal" Faraday effect, while the latter contributes more generally and describes optical activity. Although we have restricted ourselves to the independent particle approximation, and have neglected spin effects and the motion of ion cores, within these approximations our expression for σ ilj (ω) describes both optical rotary dispersion and circular dichroism.
Our approach is based on introducing microscopic polarization and magnetization fields, from which the charge and current density expectation values can be found. The corresponding macroscopic fields of elementary electrodynamics can then be defined as the spatial averages of those microscopic fields; the "free" macroscopic charge and current densities that can generally arise vanish in the linear response of the class of insulators we consider. With the use of a set of Wannier functions, we associate portions of these microscopic fields with each lattice site, thereby introducing site polarization and magnetization fields from which site multipole moments are extracted.
We then construct macroscopic multipole moments from these site multipole moments, and from their linear response to the electromagnetic field we identify the tensors describing the response of the electric dipole moment per unit volume P i (x, t) to the electric field E l (x, t), to the symmetrized derivative of the electric field F jl (x, t), and to the magnetic field B l (x, t); the response of the electric quadrupole moment per unit volume Q ij (x, t) to E l (x, t); and the response of the magnetic dipole moment per unit volume M i (x, t) to E l (x, t). From these tensors we construct σ il (ω) and σ ilj (ω). Due to its focus on identifying site quantities, our strategy allows for an easy comparison with results in the "molecular crystal limit", where the electrons associated with a molecule at one lattice site cannot move to another site. But it certainly does not require that idealization.
In the limit of uniform and static electric and magnetic fields we recover the magnetoelectric effect described earlier by others [23,24] and us [19], the latter calculation using the approach implemented here. There the first-order modifications of both P i due to B l and of M l due to E i are described by the orbital magnetoelectric polarizability (OMP) tensor α il , which is nonvanishing only if both time-reversal and spatial inversion symmetry are broken in the unperturbed system. At finite frequencies the previously identified contributions to the OMP tensor, α il G and δ il α CS , remain as contributions to the response of P i (x, ω) to B l (x, ω) and of M l (x, ω) to E i (x, ω). However, additional explicitly frequency-dependent contributions, α il P (ω) and α il M (ω), to the total response tensors emerge, generally resulting in the tensors describing the linear response of P i (x, ω) to B l (x, ω) and of M l (x, ω) to E i (x, ω) to differ. These additional contributions are classified as "cross-gap" contributions, like α il G , but are gauge dependent. Thus, the net cross-gap contributions would be given by α il G + α il P (ω) and α il G + α il M (ω), respectively. The terms α il P (ω) and α il M (ω) that arise and differentiate the responses result from contributions that we identify as "dynamical." Furthermore, as the finite frequency generalization of the "compositional" contributions to the response tensors is trivial, and because α il P (ω) and α il M (ω) are manifestly "cross-gap" contributions, the Chern-Simons contribution to the finite frequency response tensors remains unchanged; that is, the finite frequency generalization of the Chern-Simons contribution to these response tensors is identical to that in the limit of uniform and static Maxwell fields.
In the molecular crystal limit the Chern-Simons contribution, which does not contribute to the bulk macroscopic charge and current densities that are linearly induced by the Maxwell fields, becomes gauge invariant [36]. As well, in that limit the response tensor characterizing the finite frequency linear response of P i (x, ω) to B l (x, ω), and the response tensor characterizing that of M l (x, ω) to E i (x, ω), are related; this relation does not hold in general beyond the molecular crystal limit. Similarly, the relations between the tensors describing the linear response of P i (x, ω) to F jl (x, ω) and of Q ij (x, ω) to E l (x, ω) that hold in the molecular crystal limit do not hold generally. This is because in the molecular crystal limit the site multipole moments can be associated with expectation values of associated operators familiar from molecular physics, whereas for a crystal in which charges can move more freely a Green function approach was used to define them.
Generally, these macroscopic multipole moments were introduced with the use of Wannier functions associated with each lattice site, and thus P i (x, t), Q ij (x, t), and M i (x, t) are "gauge dependent" in the sense that they depend on the choice of these Wannier functions. A natural choice, of course, would be a set of ELWFs. However, we showed that whatever choice is made the expressions for the linear response of the macroscopic charge and current densities to the Maxwell fields are gauge invariant. Thus our expression for σ ilj (ω), as well as that for σ il (ω), can be evaluated without any calculation -or any thought -of the Wannier functions than underpin our approach. At frequencies below the band gap we found agreement with earlier work restricted to that frequency range [7].
Yet, while they do not appear explicitly in the final expression for σ il (ω) or σ ilj (ω), it is the site multipole moments that can be introduced with the aid of these Wannier functions, and the microscopic polarization and magnetization fields on which the approach is based, that make possible the natural connection and comparison with the molecular crystal limit. This should lead to an understanding of which features of the optical activity of any material of interest can be associated with physics beyond that limit. As well, the use of such site quan-tities in our approach offers the possibility of considering the optical response of a finite system, where simply identifying the bulk tensors σ il (ω) and σ ilj (ω) is not sufficient, and will lead to the description of other linear and nonlinear optical response features that depend on the variation of the electromagnetic field throughout a finite crystal. We plan to turn to these generalizations in future publications. The "relators" allow us to obtain the microscopic polarization and magnetization fields from the charge and current density expectation values. They also arise in the way the Maxwell fields enter the dynamical equations governing such quantities. Thus, an expansion of the relators is relevant for the identification of the electric and magnetic moments and in expanding the equations of motion of quantities associated with the electron field in terms of powers of the Maxwell fields and their derivatives. As a consequence, the expansion parameter u appearing in the relator expansions can be used to identify which perturbative modifications to the various site multipole moments due to a particular Maxwell field, or derivative of that field, appear at the same "order." We now show this.

VIII. ACKNOWLEDGMENTS
The expansions of Ω j y (x, t) and Ω 0 y (x, t), (33,34), derived previously can be more easily derived using a formal expansion of the "relators", s i (w; x, y) and α ij (w; x, y), about u = 0. We begin with the definition, under the choice of a straight-line path; see Ref. [18], where we find Recall we have previously defined and found for nearly uniform Maxwell fields Ω a y (x, ω) which we have implemented in this work. We now find these approximate expressions in a different way. We write the first of (A1) as Used in (A3), following a partial integration with respect to w, immediately gives (A5). Notice that the first term of (A5) originates from the O(u 0 ) term of the s i -relator expansion (A6), and the second term from the O(u) term. Similarly, we expand the second of (A1) to the same order, O(u), and find α ij (w; x, y)  Using this in (A2) we immediately arrive at (A4). Then the explicitly retained term of (A4) originates from an O(u) term of the α ij -relator expansion. Thus, (A4) and the second term of (A5) appear at the same order of the expansion parameter u. This is consistent with Faraday's law as the spatial derivatives of E(x, ω) are related to frequency factors times B(x, ω). Thus such terms appear at the same "order" with respect to the Maxwell fields and their derivatives kept in an expansion. It appears that the expansion parameter u captures this information. Now the site electric and magnetic multipole moments can be found from the "site" polarization and magnetization fields, respectively, using these same relator expansions. The site electric dipole moment (10) originates from the O(u 0 ) term of the s i -relator expansion (A6), while the site electric quadrupole moment (11) originates from the O(u) term of the s i -relator expansion. The site magnetic dipole moment (12) originates from the O(u) term of the α ij -relator expansion (A7). However, it is not only via the expansion of those relators that relate the microscopic charge and current densities to the microscopic polarization and magnetization fields that the expansion parameter u enters. When finding the linear response of the single-particle density matrix, (38,39,41), the quantities (A4,A5) are used. Thus, modifications to the site electric dipole moment appear at least at order O(u 0 ), but not all modifications to this quantity appear at this order; for instance, (43) appears at O(u 0 ), while (47) and (49) appear at O(u). Furthermore, modifications to the site electric quadrupole and the site magnetic dipole moments appear at least at order O(u); for example, (53) and (55)-(57) appear at O(u). In this work, we only consider the contributions to the linear response of a site quantity appearing at most at O(u); we neglect higher-order modifications, such as those leading to the magnetic susceptibility, which would appear at O(u 2 ), those related to spatial derivatives of the magnetic field, or those related to the linear response of higher-order multipole moments.

Appendix B: Microscopic and macroscopic fields
In this Appendix we describe approaches to constructing macroscopic fields from the microscopic fields appearing in (6).
One approach (I) often adopted to treat infinite crystals is to start from the Fourier transform to wavevector space of all the quantities of interest. For the current density, for example, we would have etc. If j(q, t) is nonzero only for a single, small q, then the variation in the current density is sinusoidal. If one wants to consider less trivial variations, then one needs to treat a range of qs. To do this one can introduce a macroscopic field associated with each microscopic field -e.g., J (q, t) with j(q, t) -by setting J (q, t) = j(q, t) for some restricted region of reciprocal space near the origin -say, |q| ≤ ∆ −1 , where ∆ is a length satisfying ∆ a, with a on the order of a lattice constant -and J (q, t) = 0 for other q in reciprocal space. Also choosing ∆ λ, where λ characterizes a typical range of variation of the fields that one is trying to capture, the macroscopic fields can describe excitations in the crystal characterized by typical length scales much larger than the lattice constant.
Another approach (II) with the same goal starts in position space rather than reciprocal space, and introduces a smooth weighting function w(x) to extract a macroscopic field L(x) from the associated microscopic field l(x) [37], identifying the macroscopic field at a point x with the average of the associated microscopic field in the neighborhood of x. We take w(x) to be a smooth positive function, peaking at x = 0 and spherically symmetric about that point, dropping off continuously as |x| → ∞ with a characteristic length scale ∆ satisfying the conditions given above, and with an integral over all space equal to unity. A typical example would be a Gaussian function, w(x) = w II (x), where w II (x) = e −|x| 2 /∆ 2 ∆ 3 π 3/2 .
The two approaches can be formally related, of course, because from (B1) we have L(q) = w(q)l(q), and by formally setting w I (q) = θ(∆ −1 − |q|), where θ(q) is the Heavyside step function, we recover the first approach. It has the advantage that constructing a macroscopic field from its associated microscopic field is a projection in wavevector space; thus, choosing w(x) = w I (x), if the operation (B1) is repeated there is no additional change. On the other hand, the w I (x) that results extends far beyond |x| = ∆, and as well takes on negative values. Indeed, any w(q) which, like w I (q), has a vanishing second derivative in some directionq about q = 0 will lead to a w(x) which must take on negative values, since it has a vanishing second moment. Thus the second approach, where one begins with a smooth and wellbehaved averaging function in position space (and note w II (q) = exp − |q| 2 /(4∆ 2 ) ), seems a better choice if one wants to understand the averaging physically, and with it one can envision a treatment of finite media and interfaces. In this paper we only concern ourselves with nominally infinite crystals, so the two approaches lead to essentially the same results; we indicate the small differences below, but most of what we say would apply to either.
We adopt the semiclassical approximation, where the electromagnetic field is treated classically, and in the Maxwell equations for the microscopic electric and magnetic fields, e(x, t) and b(x, t), we take ρ(x, t) and j(x, t) (6) as the microscopic charge-current density of the crystal. Using the averaging procedure (B1) to identify the macroscopic fields from their microscopic counterparts, we immediately find that those macroscopic fields satisfy the macroscopic Maxwell equations in the form ∇ · D(x, t) = 4π F (x, t), c∇ × H(x, t) = 4πJ F (x, t) + ∂ ∂t D(x, t), where D(x, t) = E(x, t)+4πP (x, t), H(x, t) = B(x, t)− 4πM (x, t), and etc. As mentioned in the text, we refer to the macroscopic fields E(x, t) and B(x, t) as the "Maxwell fields." Using the expansions (9) in the expression (7) for the total p(x, t) and m(x, t), and then spatial averaging using (B4), we find (4), where the macroscopic electric dipole moment per unit volume, electric quadrupole moment per unit volume, and magnetic dipole moment per unit volume are given by respectively. Since F (x, t) and J F (x, t) vanish in the problem at hand, upon implementing (4) in the macroscopic Maxwell equations, (B3), P i (x, t), Q ij (x, t), and M i (x, t) serve as the only source terms at this level of analysis. The remaining task is to establish the constitutive relations (5). We can do this by inserting (14) in (B5). The terms that will appear involve R w(x − R)L(R, ω), where here L(R, ω) is one of the macroscopic fields E l (R, ω), B l (R, ω), or F jl (R, ω). To investigate this kind of sum we note that Ω uc e iq·x w(q) + 1 Ω uc G =0 w(q + G)e i(q+G)·x , (B7) where the G are reciprocal lattice vectors, and we have used If we choose w(q) = w I (q), then the q that will contribute to L(R, ω) are such that the second term in the final equality of (B7) rigorously vanishes; from the first term in that expression we see that, since w I (q) acts as a projector, we will have R w(x − R)L(R, ω) = 1 Ω uc L(x, ω), exactly. If we choose w(q) = w II (q), then there will be corrections to this, since w II (q) does not act as a projector. However, the corrections will be small given that the inequalities (B2) are assumed to be satisfied, and we can redefine our local field corrections to include them. We then find that (B5,14,B8) lead to (5), the form of our constitutive relations.

Appendix C: List of response tensors
We here list all the response tensors that were found in this work. The derivation of these response tensors, including the acknowledgment of the assumptions that have been made, and the identification of the quantities they relate, is presented in Section II-IV. The response tensor χ il E (ω) is gauge invariant. For all the other tensors, the portion indicated with a breve accent is gauge invariant. .
ns ξ a sm ξ l mn .