Multipole theory of optical spatial dispersion in crystals

Natural optical activity is the paradigmatic example of an effect originating in the weak spatial inhomogeneity of the electromagnetic field on the atomic scale. In molecules, such effects are well described by the multipole theory of electromagnetism, where the coupling to light is treated semiclassically beyond the electric-dipole approximation. That theory has two shortcomings: it is limited to bounded systems, and its building blocks - the multipole transition moments - are origin dependent. In this work, we recast the multipole theory in a translationally-invariant form that remains valid for crystals. Working in the independent-particle approximation, we introduce"intrinsic"multipole transition moments that are origin independent and transform covariantly under gauge transformations of the Bloch eigenstates. Electric-dipole transitions are given by the interband Berry connection, while magnetic-dipole and electric-quadrupole transitions are described by matrix generalizations of the intrinsic magnetic moment and quantum metric. In addition to multipole-like terms, the response of crystals at first order in the wavevector of light contains band-dispersion terms that have no counterpart in molecular theories. The full response is broken down into magnetoelectric and quadrupolar parts, which can be isolated in the static limit where electric and magnetic fields become decoupled. The rotatory-strength sum rule for crystals is found to be equivalent to the topological constraint for a vanishing chiral magnetic effect in equilibrium, and the formalism is validated by numerical tight-binding calculations.


Introduction
As the wavelength of optical radiation is large compared to atomic dimensions, the interaction of light with matter is generally well described by taking the long-wavelength limit (electric-dipole approximation). In that approximation, the response of the medium to an electromagnetic perturbation is treated as local in space. When nonlocality is taken into account, the response acquires a dependence on the wavevector q of light, and this is known as spatial dispersion [1,2].
Although the effects of spatial dispersion can often be treated as small corrections, they are significant in that they lead to qualitatively new phenomena. One example is natural optical activity [1], whose most familiar manifestation is the rotation of the plane of polarization of linearly-polarized light travelling through chiral molecules in solution. Lesser-known manifestations of spatial dispersion include gyrotropic birefringence and nonreciprocal directional dichroism [3,4]; these are magneto-optical effects that occur in acentric magnetic materials without requiring a net magnetization.
Because of the fundamental and industrial importance of chiral molecules, molecular quantum theories of natural optical activity -and, by extension, of other spatiallydispersive optical effects -have been developed over many decades, on the basis of the multipole theory of electromagnetism [5][6][7][8]. This has led to the development, starting in the mid 1990s, of several ab initio methods for calculating optical rotatory dispersion and natural circular dichroism spectra of molecules. Some of those methods rely on sum-overstates formulas [9]; in others, the explicit summation over a truncated set of excited states is avoided using either static-limit [10] or finite-frequency [11] linear-response schemes, or real-time propagation approaches [11,12] (see Refs. [13][14][15] for reviews).
In this work, we develop a microscopic theory of optical spatial dispersion in crystals that is firmly rooted in the molecular multipole theory. We work in the independentparticle approximation neglecting local-field effects [21], and focus on the electronic response with frozen ions. We proceed by evaluating the optical conductivity at first order in q, including both interband and intraband contributions, and arrive at a sum-over-states expression in terms of well-defined multipole-like transition moments.
To set the stage, let us introduce the multipole transition moments as defined in the standard molecular theory [7,8]. The electric dipole (E1) appears at leading order in the multipole expansion, followed at the next order by the magnetic dipole (M1) and electric quadrupole (E2). These are the needed ingredients to describe natural optical activity, gyrotropic birefringence, and nonreciprocal directional dichroism. In the independentparticle approximation, they take the form q ab nl = −e⟨ϕ n |r a r b |ϕ l ⟩ , where ϕ n (r) and ϕ l (r) are occupied and empty energy eigenstates of the molecule, respectively, and −e is the electron charge. (The M1 transition moment also has a spin part; we omit it for now, but it will be included later.) In trying to extend the multipole theory to periodic crystals, one is faced with the problem of how to define the transition moments when the molecular orbitals ϕ n (r) are replaced by Bloch eigenstates ψ nk (r) = e ik·r u nk (r), given that the matrix elements in Eq. (1) involve the nonperiodic position operator r. For E1 transitions between nondegenerate bands n and l, there is a well-known prescription, namely d nl (k) = −e⟨u nk |i∂ k u lk ⟩ [38,39].
The situation is less clear when it comes to defining M1 and E2 transition moments in the Bloch representation. Already for molecules, their definitions in Eq. (1) are somewhat problematic, as they give values that change under a rigid shift of the coordinate system. Molecular properties should be origin independent, and for spatially-dispersive optical coefficients that is generally ensured by a cancellation between the origin dependences of different terms of the same order in the multipole expansion [6][7][8][9]. This is not entirely satisfactory from a formal standpoint, and moreover it leads to slightly origin-dependent numerical results, because the cancellation is not exact for incomplete basis sets [8,9]. In our independent-particle formulation, the optical conductivity at first order in q is written in terms of "intrinsic" multipole transition momentsĒ1,M1, andĒ2 that are origin independent and well defined for both molecules and crystals. For molecules, these modified transition moments take the form q ab nl = −e⟨ϕ n | r a − r a n +r a l /2 r b − r b n +r b l /2 |ϕ l ⟩ .
As they are defined relative to an intrinsic origin located halfway between the centersr n = ⟨ϕ n |r|ϕ n ⟩ andr l = ⟨ϕ l |r|ϕ l ⟩ of the two orbitals, the matricesd,m, andq are manifestly origin independent. For crystals, we find that they acquire "quantum-geometric" forms when expressed in terms of the cell-periodic Bloch eigenstates:d is given by the interband Berry connection, whilem andq are described by generalizations -with both intraband and interband parts -of the intrinsic orbital moment [40] and of the quantum metric [41], respectively. The intraband orbital moment and quantum metric were already known to contribute to the spatially-dispersive optical response in metals [23][24][25]; we clarify that the quantum metric appears quite generally and not just in two-band models, and identify additional Fermi-surface terms. As for the interband counterparts of the intrinsic orbital moment and quantum metric, they had not been clearly identified in previous theoretical studies of spatial dispersion in band insulators [19][20][21][22], where the optical matrix elements were written in velocity form, and without isolating the magnetic and quadrupolar parts.
The paper is organized as follows. In Sec. 2, we introduce some basic definitions and relations. Section 3 contains the derivation and analysis of our main result: an expression for the bulk optical conductivity at first order in q. The derivation is split into several steps. We start in Sec. 3.1 from the Kubo formula for the frequency-and wavevectordependent optical conductivity in the velocity gauge. That formula suffers from apparent divergences at zero frequency, and we recast it in a form that is manifestly divergence free. We then expand the Kubo formula to linear order in q, treating the q dependence coming from the band dispersion and from the matrix elements in Secs. 3.2 and 3.3, respectively. In Sec. 3.4, we convert the optical matrix elements from velocity to length form, and express them in terms of intrinsic multipole transition moments between Bloch eigenstates. In Sec. 3.5 we collect terms, and report the final expression for the optical conductivity at first order in q. Finally, in Sec. 3.6 we briefly discuss its decomposition into magnetoelectric and quadrupolar parts, and the associated physical effects in the static limit. In Sec. 4 we consider the molecular limit of our formalism, and its relation to the standard multipole theory. The formal part of the paper ends in Sec. 5 with an analysis of optical sum rules, and in Sec. 6 we present numerical results for a tight-binding model. We conclude in Sec. 7 with a summary and discussion.

Basic definitions and relations
Consider the response of a medium to a monochromatic electromagnetic field. We work in the "temporal gauge" [2], where the electromagnetic field is fully described by the vector potential A(t, r) = Re A(ω, q)e i(q·r−ωt) . To linear order in the field amplitude, the induced current density reads so that one may define an effective conductivity as σ ab (ω, q) = (1/iω)Π ab (ω, q) [1][2][3]. If the spatial dispersion is weak, the effective conductivity can be expanded as where a summation over the repeated Cartesian index c is implied. The zeroth-order term is the optical conductivity in the long-wavelength limit, that is, in the electric-dipole approximation. The next term in the expansion captures the effects of spatial dispersion to the order of magnetic dipoles and electric quadrupoles, and will be the focus of our study. Let us split it into symmetric (S) vs antisymmetric (A), and into Hermitian (H) vs anti-Hermitian (AH) parts with respect to its first two indices, The H and AH parts are absorptive and reactive, respectively [1,2]. Regarding the S and A parts, they transform differently when time reversal T is applied to the material, that is, under reversal of its magnetic order parameter. According to the Onsager reciprocity relation [1,2], the A part is T even, and the S part is T odd; the former describes natural optical activity, and the latter describes spatially-dispersive magneto-optical effects. The entire σ ab,c tensor is odd under spatial inversion P, and hence it vanishes in centrosymmetric systems. If both P and T are broken but the combined PT symmetry is present, σ A ab,c vanishes but σ S ab,c can be nonzero.
3 The bulk formula for σ ab,c (ω) 3.1 Kubo formula for the effective conductivity σ ab (ω, q) We now specialize to a three-dimensional crystal described by a single-particle Pauli Hamiltonian H [38] with a local external potential, and introduce the electromagnetic perturbation via an interaction Hamiltonian H I expressed in the velocity gauge [42]. To linear order in the vector potential A(t, r), the interaction Hamiltonian can be written as [23] H (η) where m e is the electron mass, v = (1/iℏ)[H, r] is the unperturbed velocity operator, and S is the spin operator. In addition, we have definedÃ (η) (t, r) = e ηt A(t, r), where the parameter η is formally a positive infinitesimal that controls the adiabatic turning on of the coupling between the electromagnetic field and the crystal [8,[43][44][45][46]. A standard perturbative calculation yields the following Kubo formula for the optical conductivity in the spectral representation [23,44], Here k = BZ dk/(2π) 3 , ω lnk (q) = ω l (k + q/2) − ω n (k − q/2), where ℏω n (k) = ε n (k) is the band energy, and f lnk (q) = f l (k + q/2) − f n (k − q/2), where f n (k) = f [ω n (k)] is the Fermi-Dirac occupation factor. In the first term, N is the total number of electrons per unit volume; in the second, the matrix element is defined as where I lnk (q) is a sum of orbital and spin contributions [23], In Eq. (9a), v(k) = (1/ℏ)∂ k H(k) with H(k) = e −ik·r He ik·r , and in Eq. (9b), g s ≈ 2 is the spin g-factor of the electron. Henceforth, the index k will be omitted for brevity. The 1/ω prefactors in Eq. (7), inherited from Eq. (3), make it singular at ω = 0. That singularity is only apparent [46], and it can be removed as follows. First, split Eq. (7) into reactive and absorptive parts using lim η→0 + (x − iη) −1 = 1/x + iπδ(x). Next, notice that the reactive part can be rewritten by invoking the Kramers-Krönig relation while in the absorptive part the factor 1/ω can be replaced with 1/ω ln thanks to the delta function. Finally, recombine the two parts to obtain where we have definedω = ω + iη and At zeroth order in q Eq. (11) reduces to Eq. (25) of Ref. [46], and at first order it becomes, in the notation of Eq. (4), In Sec. 5, we discuss how the equivalence between the Kubo formulas (7) and (11) at zeroth and first orders in q is related to the oscillator-strength and rotatory-strength sum rules, respectively. In the context of tight-binding calculations the diamagnetic term in Eq. (7) changes form [47], while Eq. (11) remains unchanged. In the following subsections, we evaluate the expansion coefficients of the E and M matrices appearing in Eq. (13). We start in Sec. 3.2 with the expansion of E, and then devote Secs. 3.3 and 3.4 to the expansion and subsequent manipulations of M, which is where our treatment differs more substantially from that of previous works.
The terms in the resulting expression for σ ab,c (ω) can be classified as either "band dispersive" or "molecular," depending on whether or not they vanish for a crystal composed of nonoverlapping units. The first term in Eq. (13) is clearly band dispersive, because the quantity E ,c ln (ω) involves the band velocities v n = ∂ k ω n (see Eqs. (14c) and (14d) below). While less obvious, the second term in Eq. (13) is not purely molecular; as we will see, it has a band-dispersion component that went unnoticed in previous works [19,22].

Expansion in q of the band-energy terms
When expanding Eq. (12) in powers of q, the intraband (l = n) and interband (l ̸ = n) parts must be treated separately. To first order in q, one finds where f ′ n = ∂f n /∂ω n , f ln = f l − f n , ω ln = ω l − ω n , and Z ln (ω) = 1/(ω 2 ln −ω 2 ). For the intraband identities, we used which follows from ω nn (q) = (∂ a ω n )q a + O(q 3 ) and

Nondegenerate bands
Energy eigenstates are only defined up to overall phase factors, and observable quantities cannot depend on this phase arbitrariness. In the case of nondegenerate Bloch bands, physical observables must remain invariant under single-band quantum-mechanical "gauge transformations" of the form |u n ⟩ → e −iβn |u n ⟩, where β n is a real function of k.
As the M matrix defined by Eqs. (8) and (9) is clearly gauge invariant, the same must be true for its expansion coefficients entering Eq. (13) for σ ab,c (ω). When evaluating those coefficients, we would like to insist that each individual contribution -and not just their sum -is gauge invariant. Doing so will lead to a physically transparent and numerically robust expression for σ ab,c (ω) in terms of origin-independent quantities.
The coefficient M ab (0) appearing in the first term of Eq. (13) is trivially gauge invariant, as it involves a single term, which is a product of gauge-covariant velocity matrix elements (clearly, those matrix elements are also origin independent). Instead, the coefficient M ab,c appearing in the second term comprises several terms, not all of which are individually gauge invariant. The problematic terms are those that contain matrix elements such as ⟨u n |v a |∂ c u l ⟩, because Blochstate derivatives transform noncovariantly as This can be fixed by writing |∂ k u n ⟩ as |D k u n ⟩ − iA n |u n ⟩, where |D k u n ⟩ is the covariant derivative [39] and A n = ⟨u n |i∂ k u n ⟩ is the intraband Berry connection. The terms containing Berry connections cancel out, leaving where every term is a gauge-invariant product of gauge-covariant matrix elements, just like in Eq. (16).

Degenerate bands
Gyrotropic birefringence and nonreciprocal directional dichroism occur in antiferromagnetic crystals such as Cr 2 O 3 [3,4], where the energy bands are doubly degenerate at every k as a result of the combined PT symmetry [48]. To treat such cases, we introduce degeneracy indices λ and ν for the Bloch states in bands l and n, respectively. The Kubo formula (11) remains unchanged, but the matrix element therein becomes a trace over the degeneracy indices, M ab nl = λ,ν (I a lλ,nν ) * I b lλ,nν . The reasoning leading up to Eq. (17) follows through, provided that the covariant derivative is generalized . The object |D k u nν ⟩ transforms covariantly under multiband gauge transformations of the form where U n is a k-dependent unitary matrix in the degeneracy indices. To alleviate the notation, from now on we will assume nondegenerate bands.

Conversion to length (multipole) form
As we started out from the Kubo formula in the velocity gauge, the optical matrix elements (16) and (17) entering Eq. (13) for σ ab,c (ω) are written in terms of the velocity operator. Now, we would like to recast those matrix elements in a "length form" that brings out their multipole character. In the case of molecules [9,49], this is achieved by means of identities such as ⟨ϕ l |v|ϕ n ⟩ = iω ln ⟨ϕ l |r|ϕ n ⟩.
In periodic crystals, where the velocity operator is given by the gradient of the Hamiltonian, the conversion from velocity to length form follows from the identity which can be obtained by differentiating H|u n ⟩ = ε n |u n ⟩, and then writing |∂ k u n ⟩ as |D k u n ⟩ − iA n |u n ⟩. Contracting with ⟨u l | gives v ln = δ ln v n + iω ln A ln , where v n is the band velocity, and is the interband Berry connection. With Eqs. (18) and (19), one can split Eq. (17) for M ab,c nl as where we have defined The first term in Eq. (20) is molecular for l ̸ = n and band dispersive for l = n, whereas the second term is purely band dispersive and vanishes for l = n (the distinction between molecular and band-dispersive contributions was introduced in Sec. 3.1). From the gauge-covariant and Hermitian matrices A a and B bc , we can now define for crystals the intrinsic multipole transition moments that were introduced in Eq. (2) for molecules. The intrinsic electric-dipole matrix isd a = −eA a , while the intrinsic magneticdipole and electric-quadrupole matrices are related to the antisymmetric and symmetric parts of B bc as follows,m a ln = e 2 ϵ abc B bc ln , (22a) wherem ln comprises orbital and spin contributions. By expanding the covariant derivatives and then setting l = n, one finds thatd nn = 0 [see Eq. (19)], and thatm nn and −q bc nn /e are respectively the intrinsic magnetic moment m n [40] and the quantum metric g bc n [41] of a Bloch eigenstate.
Using Eq. (19) for A ln together with the completeness relation, Eq. (23) can be recast in a more convenient form for numerical work, As they are written in terms of matrix elements of the velocity operator, these expressions are manifestly origin independent. The correspondence with the molecular expressions in Eq.
In the case of degenerate bands, Eq. (23) gets modified in the manner described in Sec. 3.3.2. The modified Eq. (24) remains nonsingular, as its energy denominator only contains energy differences between nondegenerate bands.
There is at present considerable interest in quantum-geometric quantities associated with interband optical responses [50]. In this regard, we note that the quantity −q bc ln /e is distinct from the "band-resolved quantum metric" g bc nl /2 that has been introduced in connection with nonlinear optical responses [51][52][53]. The quantity g bc ln is gauge invariant for every l and n, and when summed over l it gives the quantum metric of band n, g bc Instead, −q bc ln /e is gauge covariant for l ̸ = n, and for l = n it reduces to g bc n .

Final expression
We have now gathered all the needed ingredients to evaluate Eq. (13) for σ ab,c (ω), namely the expansion coefficients of E in Eq. (14), and those of M in Eqs. (16) and (20). In Eqs. (26) and (31) below, we break down the resulting expression into antisymmetric (Teven) and symmetric (T -odd) parts. The real and imaginary parts of those two equations are either absorptive or reactive, as per Eq. (5).
To arrive at Eqs. (26) and (31), several terms containing double band summations were eliminated by exchanging the l and n indices (note also that the l = n terms therein vanish, because f nn = ω nn = A a nn = 0). Those equations are written in terms of the A a and B bc matrices, which in turn are related to the intrinsic multipole transition moments by According to Eq. (24), the A a and B ab matrices depend exclusively on band energies and interband velocity matrix elements. From the restrictions imposed on these quantities by the presence of P, T , or PT symmetry [48], the restrictions on A a and B ab can be deduced. In this way, it may be verified that the expressions given below satisfy the symmetry constraints discussed at the end of Sec. 2.

Antisymmetric (time-even) part
The antisymmetric part of σ ab,c (ω) takes the form The five terms in this expression can be classified as follows. The first is molecular, while the others are band dispersive; the four inside curly brackets are interband, while the fifth is intraband; and the first three are Fermi-sea-like, while the last two are Fermi-surface-like.
In insulators and cold semiconductors, only the Fermi-sea terms survive, and one can compare with previous treatments of optical activity in nonconducting crystals. In Refs. [19,22], the sole band-dispersion contribution to σ A ab,c (ω) came from differentiating the E matrix, that is, from the first term of Eq. (13) for σ ab,c (ω). It went unnoticed in those works that the other term in that equation -where one differentiates the M matrix instead -is not purely molecular, as shown in Eq. (20). This is why we have not two but three Fermi-sea terms in Eq. (26), one molecular and two band dispersive.
In conductors, the Fermi-surface terms contribute as well. Using Eq. (25), the last term in Eq. (26) becomes (ϵ acd K bd − ϵ bcd K ad ) /(eω), where This intraband contribution to optical activity involving the intrinsic magnetic moment of conduction electrons was identified in Refs. [23,24], and was evaluated for p-doped tellurium in Refs. [28,54]. The fourth term in Eq. (26) gives an additional interband contribution to the optical activity of conductors that was overlooked in previous works. The low-frequency behavior of the optical rotatory dispersion is different in insulators and in conductors. For simplicity, let us consider the propagation of light along the optical axis z of a uniaxial crystal. The rotatory power is given by [1] where ϵ 0 is the vacuum permittivity and c is the speed of light. To deal with absorption, the positive infinitesimal η inω = ω + iη has been reinterpreted heuristically as a phenomenological scattering rate τ −1 [6,46,55]. For frequencies and scattering rates well below the threshold for interband transitions, ω, τ −1 ≪ ω gap , Eq. (26) yields ρ(ω, τ ) = (ωτ ) 2 1 + (ωτ ) 2 a + bω 2 .
The coefficient b comes from the interband terms which haveω prefactors, and a = −(e/c 2 ϵ 0 ℏ)K xx (with K xx = K yy ) comes from the intraband term with a 1/ω prefactor. In insulators the coefficient a vanishes, and hence the rotatory power displays the familiar ω 2 dependence at low frequencies [6]; in conductors that coefficient is nonzero, and one can distinguish two different regimes as follows, In Sec. 6, we will illustrate these low-frequency profiles for a concrete tight-binding model.

Symmetric (time-odd) part
The symmetric part of σ ab,c (ω) reads The first three terms are Fermi-sea-like, and can be compared with the expressions obtained for insulators in Ref. [22]. The third is band dispersive, and it corresponds to Eq. (31) in that work, while the first two add up to Eq. (30) therein, revealing its mixed molecular/dispersive character. The remaining two terms in Eq. (31) are Fermi-surface-like, and they can be compared with the expressions obtained for metals in Ref. [25]. The last term was identified in that work. Writing ω 2 ln Z ln (ω) as 1 +ω 2 Z ln (ω) and noting that Re l A a nl A b ln is the quantum metric g ab n = −q ab nn /e [this can be seen from Eqs. (19) and (24c)], the fourth term in Eq. (31) splits into intraband and interband parts as follows, The intraband piece is similar to the quantity K ab in Eq. (27), but with the intrinsic magnetic moment replaced by the quantum metric (intrinsic quadrupole moment). An equivalent result was obtained in Ref. [25] for two-band models, but without invoking the identities leading up to Eq. (32), which is what allowed us to isolate a quantum-metric contribution in the general multiband case.

Magnetoelectric and quadrupolar responses
As recognized in Ref. [3], the physical basis for σ ab,c (ω) is provided by quadrupolar and magnetoelectric couplings. The quadrupolar couplings are described by a T -odd totallysymmetric tensor γ abc (ω), and the magnetoelectric ones by a T -odd traceless tensorα ab (ω) together with a T -even tensorα ab (ω). Those tensors are defined as [3,22] γ abc (ω) = 1 3i σ S ab,c (ω) + σ S bc,a (ω) + σ S ca,b (ω) , (33a) so that Inserting Eq. (34) in Eq. (3) for the induced current density and then using the Maxwell-Faraday equation yields the constitutive relation where X ω,q denotes X(ω, q)e i(q·r−iωt) , and M ω,q a = [α ba (ω) +α ba (ω)] E ω,q b . In the quasi-static limit the electric and magnetic fields become decoupled, and the three terms in Eq. (35) describe separate physical responses. The first corresponds to a magnetization current induced by a uniform electric field; the second to a current induced by a time-varying magnetic field; and the third to a current induced by a spatially-varying electric field. The first two are direct and inverse magnetoelectric effects [1], and the third is an electric quadrupolar effect.
The T -even magnetoelectric tensorα ab describes a "kinetic" magnetoelectric effect in gyrotropic conductors [56], while the T -odd tensorα ab describes magnetoelectric effects in both insulators and conductors. In the case of insulators, the inverse magnetoelectric response can be expressed as a polarization current, ∂ t P ω,q a =α ab (ω)∂ t B ω,q b ; integrating the adiabatic current in the quasi-static limit, one obtains the familiar form P a =α ab B b of the inverse magnetoelectric effect [1]. The full T -odd magnetoelectric tensor includes an additional trace ("axion") piece (θe 2 /2πh)δ ab ; the axion angle θ is boundary sensitive [39], and this is why it is not captured by the present formalism, which is based on the electromagnetic response of a bulk medium.
Quantum-mechanical expressions for the static magnetoelectric and quadrupolar susceptibilities can be obtained by inserting in Eq. (33) the formulas given in Sec. 3.5 for σ S ab,c (ω) and σ A ab,c (ω), evaluated at ω = 0. This has been done in previous works for selected terms only: in Ref. [22], the Fermi-sea orbital contribution toα ab (0) was shown to give the traceless part of the orbital magnetoelectric susceptibility tensor of insulators [57]; and in Ref. [23], the intraband contribution toα ab (0) was shown to reproduce the kinetic magnetoelectric susceptibility obtained from Boltzmann transport in the relaxation time approximation combined with the modern theory of orbital magnetization [58].
Compared to those previous works, the present formalism provides a more complete description. In particular, it captures T -odd magnetoelectric and electric-quadrupolar effects in conductors that so far have only been treated semiclassically [25,[59][60][61], and which were found to involve the quantum metric. A detailed account of magnetoelectric and electric-quadrupolar responses on the basis of the present formalism will be given in a separate work.

Molecular limit
In this section, we analyze the molecular limit of our formalism. First, we show how in that limit the bulk expressions for the intrinsic transition momentsd,m, andq reduce to those in Eq. (2). We then show how the formulas for σ A ab,c (ω) and σ S ab,c (ω) reduce to the standard molecular expressions in terms of the ordinary transition moments d, m, and q in Eq. (1).
Consider an idealized molecular crystal composed of nonoverlapping units. For such a crystal, the cell-periodic Bloch states assume the form [62,63] u nk (r) where ξ(r) = r − R(r) is the intracell coordinate, R(r) is the lattice vector that folds the absolute coordinate r into the home unit cell, ϕ n (r) is vanishingly small outside that cell, and .
= denotes an equality that only holds in the molecular limit. The intraband Berry connection can now be easily evaluated by integrating over the home cell, and the covariant derivative of a Bloch state reduces to for r in the home cell. Using this identity in Eq. (23) ford,m, andq, we recover after some manipulations the expressions in Eq. (2). In the case ofm, it is necessary to invoke the operator identity [v a , With these relations in hand, we can address the molecular limit of Eqs. (26) and (31) for σ A ab,c (ω) and σ S ab,c (ω). Since the energy bands become dispersionless in that limit, all band-dispersion terms in those equations vanish, leaving only the first term in each of them; and since the transition moments also become independent of k, we can set k . = 1/V c in those terms (V c is the cell volume) to find where we dropped theω dependence for brevity, and defined the (extensive) molecular tensorsḠ (40d) Writing f nl as f n (1 − f l ) − f l (1 − f n ), the f nl factors in the expressions above can be replaced with 2f n (1 − f l ). In that form,Ḡ,Ḡ ′ andā become single-particle versions of the multipolar susceptibility tensors G, G ′ and a defined in Eqs. (2.83), (2.85) and (2.86) of Ref. [7], with one difference: the ordinary transition moments d, m, and q have been replaced withd,m, andq. It is not immediately clear that the same is true forā ′ , since Eq. (40d) contains a factor of ω 2 ln /ω in place of theω factor appearing in Eq. (2.84) of Ref. [7] for a ′ . However, those factors are interchangeable in the expression for a ′ , as can be seen in the manner described around around Eqs. (2.75-2.78) of Ref. [7]. Thus, (Ḡ,Ḡ ′ ,ā,ā ′ ) are origin-independent versions of the molecular tensors (G, G ′ , a, a ′ ) entering the standard multipole theory.
Next, let us consider the propagation of light inside our idealized molecular crystal. For a given propagation directionn, we define (intensive) optical-activity and gyrotropicbirefringence tensors as β A ab (n) = −σ A ab,cn c and β S ab (n) = iσ S ab,cn c , respectively. Using Eq. (39), we obtain Eqs. (5.8) and (5.9) of Ref. [7] for those tensors, but with (Ḡ,Ḡ ′ ,ā,ā ′ ) in place of (G, G ′ , a, a ′ ). Inserting Eq. (2) in Eq. (40), the terms containing the orbital centers drop out from the combinations of molecular tensors appearing in Eq. (39). Thus, (d,m,q) can be safely replaced by (d, m, q) for the purpose of evaluating the optical properties of an idealized molecular crystal. This completes the proof that our formalism correctly reduces to the standard single-particle multipole theory in the molecular limit.
In summary, we have in Eqs. (39) and (40) a reformulation of the molecular multipole theory of optical spatial dispersion at linear order in q in terms of translationally-invariant property tensors. This is in contrast to the standard formulation, where translational invariance is achieved by a cancellation between the origin dependences of the magneticdipole and electric-quadrupole terms [6][7][8][9].

Sum rules
In Sec. 3.1, we wrote two alternative Kubo formulas for σ ab (ω, q), namely Eqs. (7) and (11). The former displays apparent 1/ω divergences at ω = 0, whereas the latter is explicitly divergence free. In this section, we scrutinize the mathematical identities that underlie the equivalence between those two formulas at zeroth and first order in q, and relate those identities to the oscillator-and rotatory-strength sum rules.

Equivalence between the two forms of the Kubo formula
Let us denote as (ie 2 /ω)∆ ab (q) the difference between the reactive parts of the Kubo formulas (7) and (11).
For the two Kubo formulas to be equivalent, ∆ ab (q) must vanish identically, and according to the derivation in Sec. 3.1 this is guaranteed by the Kramers-Krönig relations. To analyze the behavior of ∆ ab (q) at zeroth and first order in q, we expand it as Let us start with the zeroth-order term in the expansion. Writing the electron density N in the first term of Eq. (41) as n k f n , and using the identity in Eq. (15), followed by an integration by parts, to deal with the l = n contribution to the second term, we obtain The quantity in square brackets is formally real and symmetric in accordance with Eq. (43a), and it vanishes identically by virtue of the effective-mass theorem. We note that the same theorem was invoked in Ref. [64] to remove the apparent divergence at ω = 0 in the dielectric function of insulators and semiconductors.
In preparation for analyzing the first-order term in Eq. (44), let us compare Eq. (41) for ∆ ab (q) with the ω → 0 limit of Π ab (ω, q) = iωσ ab (ω, q) [see Eq. (3)], evaluated using the Kubo formula (7). This gives −e 2 ∆ ab (q) = Π ab (0, q), and so −e 2 ∆ ab,c = ∂ qc Π ab (0, q)| q=0 . From the analysis of this quantity in Ref. [23] (see Sec. III.A.1 of its Supplemental Material), we conclude that where Ω n is the Berry curvature. In agreement with Eq. (43b), the expression above is purely imaginary and antisymmetric in a and b. It vanishes identically for topological reasons, and that amounts to a no-go theorem for the chiral magnetic effect in equilibrium [23].
In a recent work [36], an expression was derived for σ ab,c (ω) that contains a term diverging as 1/ω. The authors found that the prefactor of that term was neglibible for a specific tight-binding model, but they were unable to confirm analytically that it vanishes in general. The removal of that apparent divergence can be achieved by means of Eq. (46).
The vanishing of ∆ ab (q) dictates the high-frequency behavior of the optical conductivity as follows. Suppose there is a frequency ω max above which the system does not absorb [6]; setting ω ≫ ω max in Eqs. (11) and (12) and comparing with Eq. (41), we can deduce that Thus, at high frequencies the optical conductivity reduces to the diamagnetic term; and since that term is independent of q, we conclude that σ ab (ω, 0) decays as 1/ω. In Sec. 5.2.4 of Ref. [6], the high-frequency behavior of the optical activity of molecules was inferred from the rotatory-strength sum rule. This is consistent with the present analysis, because that sum rule is a direct consequence of the vanishing of ∆ ab,c , as we will now show.

Optical sum rules
Consider the sum rules obtained by integrating over positive frequencies the absorptive part of the optical conductivity, taking into account both interband and intraband absorption. Writing ∞ 0 f (ω)dω as ⟨f (ω)⟩ and using Eq. (5) yields σ H ab (ω, q) = Re σ S ab (ω, q) + i Im σ A ab (ω, q) .
To evaluate this quantity, we begin by taking the Hermitian part of Eq. (11) for η → 0 + , Making the substitution and noting that at zero temperature only the second term contributes to Eq. (49) when ω > 0 and q ≈ 0, we obtain where we have defined Let us split R ab (q) and I ab (q) in Eq. (51) into even and odd parts in q. Using the identities in Eq. (42), one finds that the even part of R ab (q) plus the odd part of iI ab (q) is proportional to the second term in Eq. (41) for ∆ ab (q). Therefore, where we keep track of the vanishing quantity ∆ ab (q). The expansion of Eq. (53) in powers of q generates a series of sum rules. Since, according to Eq. (43), the terms Re ∆ ab and Im ∆ ab only contribute (formally) at even and odd orders in q, respectively, and since the reverse is true for the second terms in Eqs. (53a) and (53b), we obtain to linear order in q. Below, we consider each of these identities in turn.
Equation (54b) is the rotatory-strength sum rule for magnetic circular dichroism. At q = 0, the intraband part of Eq. (52b) vanishes because M ab nn (0) is real, and from the interband part we recover the bulk expression given in Ref. [65] for that sum rule. If a single band is occupied, the integrated magnetic circular dichroism spectrum is proportional to the intrinsic orbital magnetic moment of the Bloch states in that band [65,66]. Equation (54c) is a sum rule for nonreciprocal directional dichroism. An explicit expression can be obtained by expanding Eq. (52a) to first order in q.
Finally, by setting ∆ ab,c = 0 in Eq. (54d) in accordance with Eq. (46), we arrive at the rotatory-strength sum rule for natural circular dichroism, This sum rule is well known for molecules in solution [5], as well as for oriented molecules [6].
Here, we have relied on a topological argument [23] to show that it remains valid for crystals, both insulating and conducting. Alternative discussions restricted to insulators are given in Refs. [19,20]. The above derivation highlights the connection between the oscillator-strength and natural rotatory-strength sum rules for crystals, and the equivalence between the Kubo formulas (7) and (11) -that is, the vanishing of ∆ ab (q) -at zero and first order in q, respectively. More generally, the expansion of Eq. (53) in powers of q yields at each order two optical sum rules, one of which relies on the vanishing of ∆ ab (q) at that order.

A tight-binding example
In this section, we use numerical tight-binding calculations to validate our formalism, and to illustrate the distinctive low-frequency profiles of the rotatory power in insulators and conductors.
As a simple model of a bulk crystal with nonzero σ A ab,c , we take the tight-binding model of Ref. [58], which consists of honeycomb layers coupled by a chiral pattern of interlayer hoppings. To break time-reversal symmetry, so that σ S ab,c becomes nonzero as well, we add complex intralayer hoppings. The resulting model is depicted in Fig. 1(a), and its Hamiltonian reads The first term is a staggered on-site potential, with ξ i = ±1 for the two sublattices in each layer. The second and third terms describe intralayer hoppings between nearest-neighbor sites i and j: the second is the complex hopping responsible for breaking time reversal, and the third is a spin-orbit coupling term; therein, σ is the vector of Pauli matrices and δ ij is the vector taking from site j to site i. The last term is the helical pattern of interlayer hoppings that renders the model chiral, with d ij the vector taking from site j to site i in adjacent layers. We choose the distance a between nearest-neighbor sites on the same layer as the unit of length, and the nearest-neighbor hopping amplitude t 1 as the unit of energy. For our tests, we set c = 1, ∆ = 0.5, λ 1 = −0.06, and λ 2 = 0.05. Exploiting the translational symmetry of the crystal, we replace the site indices {i} with {Ri}, where the lattice vector R labels the cell, and i is now an intracell site index. The Hamiltonian matrix elements are denoted by H ij (R) = ⟨ϕ 0i |H|ϕ Rj ⟩, where ϕ Rj (r) = φ j (r − R − τ j ) is a basis orbital centered at R − τ j [39]. The tight-binding Hamiltonian in k space is constructed as leading to the eigenvalue equation H k · C nk = ε nk C nk .
The energy bands of the model are displayed in Fig. 1(b). There are two composite groups with two bands each, separated by a gap. We treat the lowest group as occupied, and calculate σ A ab,c and σ S ab,c at zero temperature using Eqs. (26) and (31), respectively. The A a and B ab matrices are evaluated from Eqs. (24) and (25), with effective velocity matrix elements given by [47] v nl (k) = 1 As the system is insulating, the only nonzero contributions to σ ab,c come from the Fermisea terms in Eqs. (26) and (31); we will restrict our calculations to frequencies well below the threshold for interband absorption, where one can safely set η = 0 in those equations. The magnetic point group of the model is 32. Of the four independent tensor components that are allowed by symmetry [67,68], σ A yz,x , σ A xy,z , σ S xx,y and σ S xz,y , only the first three are actually nonzero when the Fermi level ε F lies in the gap. Converged results, obtained by sampling the Brillouin zone on a uniform mesh of 50 × 50 × 50 k points, are shown as solid lines in Fig. 2.
For comparison, we show as filled circles in Fig. 2 the results obtained from calculations on finite crystallites. We treat them as "molecules," and evaluate σ A ab,c and σ S ab,c from where f 0 is the extrapolated value to be compared with the result of the bulk calculation, and f 1 /L, f 2 /L 2 , and f 3 /L 3 account for face, edge, and corner corrections, respectively [69]. The excellent agreement seen in Fig. 2 between the two types of calculations confirms the validity of our formalism for band insulators. In Fig. 3, we take the σ A xy,z curves from Fig. 2 and split them into origin-independent contributions. For the extrapolated crystallites (left panel), there are two types of contributions in Eq. (39): those containingm are denoted asM1, and those containingq are denoted asĒ2. For the bulk crystal (right panel), there are, in addition toM1 and E2 contributions from the first term in Eq. (26), band-dispersion contributions from the second and third terms, which are denoted as v.
A comparison between the two panels of Fig. 3 reveals that theM1 andĒ2 contributions are different for the extrapolated crystallites and for the bulk crystal, with the difference being exactly compensated by the v contributions that are only present in the latter. A similar situation occurs for the ground-state orbital magnetization of a crystal, whose bulk expression contains a subtle Berry-curvature term without an obvious counterpart in the molecular theory [69,70]. For a two-dimensional insulator with a single valence band, that term reads −(e/2ℏ) k ε k (∂ x A y − ∂ y A x ) or, after integrating by parts, (e/2) k (v x A y − v y A x ). As in Fig. 3, this additional band-dispersion contribution must be included in the bulk calculation to recover the net orbital magnetization of a large flake [69,70].
To conclude, let us illustrate the different low-frequency behaviors of the rotatory power in insulating and conducting states of our model. We evaluate ρ(ω, τ ) for the bulk model from Eqs. (26) and (28), setting ℏ/τ = 2 × 10 −3 . The frequency range is chosen as 0 ≤ ℏω ≤ 10 −2 , and the calculation is carried out at zero temperature for ε F = 0.0 (insulator) and ε F = 1.0 (metal). In both cases, a uniform mesh of 100×100×100 k points  Fig. 2 into three types of origin-independent contributions: intrinsic magnetic dipole (M1), intrisic electric quadrupole (Ē2), and band dispersive (v). The latter is only present in the bulk calculation on the right, and it must be included to obtain the same total result as in the extrapolated crystallite calculations on the left. is used to sample the Brillouin zone; to improve the convergence of the calculation in the metallic case, the Fermi-surface terms in Eq. (26) are evaluated as Fermi-sea integrals by performing an integration by parts.
The resulting ρ(ω) profiles, plotted in Fig. 4 as solid lines, display the behavior dictated by Eq. (29). The insulator displays a simple ρ ∝ ω 2 decay with a negligible influence from the scattering time τ . Instead, for the metal one can distinguish two different parabolic regimes (dashed lines) delimited by ωτ ∼ 1, in accordance with Eq. (30). This distinctive low-frequency profile of the optical rotatory dispersion in conducting crystals awaits experimental verification.

Summary and discussion
In summary, we have developed a band-theoretical description of optical spatial dispersion in insulating and conducting crystals. The novelty with respect to previous formulations resides in the fact that the current induced by the optical field is given in a physicallytransparent form, as a sum of contributions that are individually origin independent, and which remain invariant under single-band gauge transformations of the Bloch eigenstates. Although we have focused on the optical conductivity σ ab,c (ω) at first order in q, higherorder responses can in principle be treated in a similar manner.
For a crystal composed of nonoverlapping units, our formula for σ ab,c (ω) reduces to the standard multipole-theory expression for molecules, but with the transition moments (d, m, q) of Eq. (1) replaced by their intrinsic (origin-independent) counterparts (d,m,q) given by Eq. (2). Away from the molecular limit, σ ab,c (ω) changes in two ways. First, the intrinsic transition moments between delocalized Bloch eigenstates are no longer given by Eq. (2); instead, one should use either the quantum-geometric expressions in Eq. (23), or their sum-over-states counterparts in Eq. (24). The second change is that σ ab,c (ω) acquires additional band-dispersion contributions associated with electron transfer between crystal cells; this is in line with the modern theories of electric polarization and orbital magnetization in crystals [39]. There were two key aspects to our derivation. The first was the use of covariant Bloch-state derivatives to expand the optical conductivity in powers of q; this allowed us to eliminate spurious noncovariant terms, and to isolate the physically relevant matrix elementsd,m, andq. The second was the choice of the nonsingular form of the Kubo formula in Eq. (11), rather than Eq. (7), as the starting point for the expansion in q. An order-by-order analysis of the equivalence between the Kubo formulas (7) and (11) led us to identify a hierarchy of optical sum rules. In particular, we found that the well-known rotatory-strength sum rule from molecular physics remains valid for crystals, thanks to a topological argument involving the k-space Berry curvature.
Our work opens up new prospects for realistic ab initio calculations of spatiallydispersive optical responses in crystals, and of static magnetoelectric and quadrupolar responses as well. An implementation based on the sum-over-states formulas for (d,m,q) in Eq. (24) has already been carried out in a concurrent work done in coordination with the present one [71], and we also envision evaluating the k-derivative formulas in Eq. (23) using Wannier interpolation [72].
In closing, we mention a possible connection with the theory of the orbital Hall effect. It was recently proposed [73][74][75] to evaluate the orbital Hall conductivity using an orbitalcurrent operator whose matrix elements in the Bloch-eigenstate basis are proportional to l m z,orb ml v a ln + v a ml m z,orb ln . Here, m orb ln is the bulk generalization of Eq. (1b), given by the same expression as in Eq. (23b) form orb ln , but with the covariant derivatives therein replaced by ordinary derivatives. The two definitions are related bȳ and therefore they agree for l = n only. For l ̸ = n, the two terms on the right-hand side of Eq. (61) are not separately gauge covariant, and the lack of gauge covariance of m orb ln makes the orbital Hall conductivity gauge dependent. This suggest that one should generally usē m orb ln instead of m orb ln when evaluating the orbital Hall conductivity. In other words, one is allowed to work with m orb ln only in the parallel-transport gauge, where A l = A n = 0.