Multipole theory of optical spatial dispersion in crystals

Óscar Pozo1, Ivo Souza1,2

1 Centro de Física de Materiales, Universidad del País Vasco, 20018 San
Sebastián, Spain

2 Ikerbasque Foundation, 48013 Bilbao, Spain

## Abstract

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.

###### Contents

  1. 1 Introduction
  2. 2 Basic definitions and relations
  3. 3 The bulk formula for $\sigma_{ab,c}(\omega)$
    1. 3.1 Kubo formula for the effective conductivity $\sigma_{ab}(\omega,{\bf q})$
    2. 3.2 Expansion in ${\bf q}$ of the band-energy terms
    3. 3.3 Expansion in ${\bf q}$ of the optical matrix elements
      1. 3.3.1 Nondegenerate bands
      2. 3.3.2 Degenerate bands
    4. 3.4 Conversion to length (multipole) form
    5. 3.5 Final expression
      1. 3.5.1 Antisymmetric (time-even) part
      2. 3.5.2 Symmetric (time-odd) part
    6. 3.6 Magnetoelectric and quadrupolar responses
  4. 4 Molecular limit
  5. 5 Sum rules
    1. 5.1 Equivalence between the two forms of the Kubo formula
    2. 5.2 Optical sum rules
  6. 6 A tight-binding example
  7. 7 Summary and discussion

##  1 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 ${\bf 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 spatially-dispersive 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-
over-states 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
[12, 11] (see Refs. [13, 14, 15] for reviews).

By comparison, there have been relatively few attempts to formulate bulk
theories of optical spatial dispersion [16, 17, 18], particularly within the
one-electron band picture [19, 20, 21, 22, 23, 24, 25, 26]. As a result, only
a small number of ab initio calculations of natural optical activity [20, 21,
27, 28, 29] and of nonreciprocal directional dichroism [30, 31] have been
carried out for crystals. Such effects provide valuable information about
broken structural and magnetic symmetries, and their study in novel bulk [32,
30, 31] and quasi-two-dimensional [33, 34, 35, 36, 37] materials calls for
improved theoretical descriptions.

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 independent-particle 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 ${\bf 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 independent-particle
approximation, they take the form

$\displaystyle{\bf d}_{nl}$ $\displaystyle=-e\langle\phi_{n}|{\bf
r}|\phi_{l}\rangle\,,$ (1a) $\displaystyle{\bf m}_{nl}$
$\displaystyle=-\frac{e}{2}\langle\phi_{n}|{\bf r}\times{\bf
v}|\phi_{l}\rangle\,,$ (1b) $\displaystyle q^{ab}_{nl}$
$\displaystyle=-e\langle\phi_{n}|r_{a}r_{b}|\phi_{l}\rangle\,,$ (1c)

where $\phi_{n}({\bf r})$ and $\phi_{l}({\bf 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 $\phi_{n}({\bf r})$ are replaced by Bloch eigenstates $\psi_{n{\bf
k}}({\bf r})=e^{i{\bf k}\cdot{\bf r}}u_{n{\bf k}}({\bf r})$, given that the
matrix elements in Eq. (1) involve the nonperiodic position operator ${\bf
r}$. For E1 transitions between nondegenerate bands $n$ and $l$, there is a
well-known prescription, namely ${\bf d}_{nl}({\bf k})=-e\langle u_{n{\bf
k}}|i{\boldsymbol{\partial}}_{\bf k}u_{l{\bf k}}\rangle$ [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 [9, 8].

In our independent-particle formulation, the optical conductivity at first
order in ${\bf q}$ is written in terms of “intrinsic” multipole transition
moments $\bar{\text{E}}1$, $\bar{\text{M}}1$, and $\bar{\text{E}}2$ that are
origin independent and well defined for both molecules and crystals. For
molecules, these modified transition moments take the form

$\displaystyle\bar{\bf d}_{nl}$ $\displaystyle=-e\langle\phi_{n}|{\bf
r}-\left(\bar{\bf r}_{n}+\bar{\bf r}_{l}\right)/2|\phi_{l}\rangle\,,$ (2a)
$\displaystyle\bar{\bf m}_{nl}$
$\displaystyle=-\frac{e}{2}\langle\phi_{n}|\left[{\bf r}-\left(\bar{\bf
r}_{n}+\bar{\bf r}_{l}\right)/2\right]\times{\bf v}|\phi_{l}\rangle\,,$ (2b)
$\displaystyle\bar{q}^{ab}_{nl}$
$\displaystyle=-e\langle\phi_{n}|\big{[}r_{a}-\big{(}\bar{r}^{a}_{n}+\bar{r}^{a}_{l}\big{)}/2\big{]}\big{[}r_{b}-\big{(}\bar{r}^{b}_{n}+\bar{r}^{b}_{l}\big{)}/2\big{]}|\phi_{l}\rangle\,.$
(2c)

As they are defined relative to an intrinsic origin located halfway between
the centers $\bar{\bf r}_{n}=\langle\phi_{n}|{\bf r}|\phi_{n}\rangle$ and
$\bar{\bf r}_{l}=\langle\phi_{l}|{\bf r}|\phi_{l}\rangle$ of the two orbitals,
the matrices $\bar{d}$, $\bar{m}$, and $\bar{q}$ are manifestly origin
independent. For crystals, we find that they acquire “quantum-geometric” forms
when expressed in terms of the cell-periodic Bloch eigenstates: $\bar{d}$ is
given by the interband Berry connection, while $\bar{m}$ and $\bar{q}$ 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 ${\bf q}$. The derivation is split into several steps. We start in
Sec. 3.1 from the Kubo formula for the frequency- and wavevector-dependent
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 ${\bf q}$,
treating the ${\bf 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 ${\bf 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.

##  2 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 ${\bf A}(t,{\bf r})={\rm Re}\left[{\bf
A}(\omega,{\bf q})e^{i\left({\bf q}\cdot{\bf r}-\omega t\right)}\right]$. To
linear order in the field amplitude, the induced current density reads

$j_{a}(\omega,{\bf q})=\Pi_{ab}(\omega,{\bf q})A_{b}(\omega,{\bf
q})=\frac{1}{i\omega}\Pi_{ab}(\omega,{\bf q})E_{b}(\omega,{\bf q})\,,$ (3)

so that one may define an effective conductivity as $\sigma_{ab}(\omega,{\bf
q})=(1/i\omega)\Pi_{ab}(\omega,{\bf q})$ [1, 3, 2]. If the spatial dispersion
is weak, the effective conductivity can be expanded as

$\sigma_{ab}(\omega,\mathbf{q})=\sigma_{ab}(\omega,{\mathbf{0}})+\sigma_{ab,c}(\omega)q_{c}+\mathcal{O}(q^{2})\
,$ (4)

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,

$\displaystyle\sigma_{ab,c}^{\mathrm{S}}$ $\displaystyle={\rm
Re}\,\sigma_{ab,c}^{\mathrm{H}}+i\,{\rm Im}\,\sigma_{ab,c}^{\mathrm{AH}}\ ,$
(5a) $\displaystyle\sigma_{ab,c}^{\mathrm{A}}$ $\displaystyle={\rm
Re}\,\sigma_{ab,c}^{\mathrm{AH}}+i\,{\rm Im}\,\sigma_{ab,c}^{\mathrm{H}}\ .$
(5b)

The H and AH parts are absorptive and reactive, respectively [1, 2]. Regarding
the S and A parts, they transform differently when time reversal $\mathcal{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
$\mathcal{T}$ even, and the S part is $\mathcal{T}$ odd; the former describes
natural optical activity, and the latter describes spatially-dispersive
magneto-optical effects. The entire $\sigma_{ab,c}$ tensor is odd under
spatial inversion $\cal{P}$, and hence it vanishes in centrosymmetric systems.
If both $\cal{P}$ and $\mathcal{T}$ are broken but the combined
$\cal{P}\mathcal{T}$ symmetry is present, $\sigma^{\text{A}}_{ab,c}$ vanishes
but $\sigma^{\text{S}}_{ab,c}$ can be nonzero.

##  3 The bulk formula for $\sigma_{ab,c}(\omega)$

###  3.1 Kubo formula for the effective conductivity $\sigma_{ab}(\omega,{\bf
q})$

We now specialize to a three-dimensional crystal described by a single-
particle Pauli Hamiltonian $\mathcal{H}$ [38] with a local external potential,
and introduce the electromagnetic perturbation via an interaction Hamiltonian
$\mathcal{H}_{\rm I}$ expressed in the velocity gauge [42]. To linear order in
the vector potential ${\bf A}(t,{\bf r})$, the interaction Hamiltonian can be
written as [23]

$\mathcal{H}_{\rm I}^{(\eta)}(t,{\bf r})=\frac{e}{2}\left[\tilde{\bf
A}^{(\eta)}(t,{\bf r})\cdot{\bf v}+{\bf v}\cdot\tilde{\bf A}^{(\eta)}(t,{\bf
r})\right]+\frac{e}{m_{e}}{\bf S}\cdot\boldsymbol{\partial}_{\bf
r}\times\tilde{\bf A}^{(\eta)}(t,{\bf r})\,,$ (6)

where $m_{e}$ is the electron mass, ${\bf v}=(1/i\hbar)[\mathcal{H},{\bf r}]$
is the unperturbed velocity operator, and ${\bf S}$ is the spin operator. In
addition, we have defined $\tilde{\bf A}^{(\eta)}(t,{\bf r})=e^{\eta t}{\bf
A}(t,{\bf r})$, where the parameter $\eta$ 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],

$\sigma^{(\eta)}_{ab}(\omega,{\bf
q})=\delta_{ab}\dfrac{ie^{2}N}{m_{e}\omega}+\dfrac{ie^{2}}{\hbar\omega}\sum_{n,l}\int_{\bf
k}\dfrac{f_{ln{\bf k}}({\bf q})}{\omega_{ln{\bf k}}({\bf
q})-\omega-i\eta}\,\mathcal{M}_{nl{\bf k}}^{ab}({\bf q})\ .$ (7)

Here $\int_{\bf k}=\int_{\mathrm{BZ}}d{\bf k}/(2\pi)^{3}$, $\omega_{ln{\bf
k}}({\bf q})=\omega_{l}({\bf k}+{\bf q}/2)-\omega_{n}({\bf k}-{\bf q}/2)$,
where $\hbar\omega_{n}({\bf k})=\varepsilon_{n}({\bf k})$ is the band energy,
and $f_{ln{\bf k}}({\bf q})=f_{l}({\bf k}+{\bf q}/2)-f_{n}({\bf k}-{\bf
q}/2)$, where $f_{n}({\bf k})=f[\omega_{n}({\bf 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

$\mathcal{M}_{nl{\bf k}}^{ab}({\bf q})=\left[I^{a}_{ln{\bf k}}({\bf
q})\right]^{*}I^{b}_{ln{\bf k}}({\bf q})\ ,$ (8)

where ${\bf I}_{ln{\bf k}}({\bf q})$ is a sum of orbital and spin
contributions [23],

$\displaystyle{\bf I}_{ln{\bf k}}^{\mathrm{orb}}({\bf q})$
$\displaystyle=\langle u_{l}({\bf k}+{\bf q}/2)|{\bf v}({\bf k})|u_{n}({\bf
k}-{\bf q}/2)\rangle\ ,$ (9a) $\displaystyle{\bf I}_{ln{\bf
k}}^{\mathrm{spin}}({\bf q})$ $\displaystyle=\dfrac{ig_{s}}{2m_{e}}\langle
u_{l}({\bf k}+{\bf q}/2)|{\bf S}|u_{n}({\bf k}-{\bf q}/2)\rangle\times{\bf q}\
.$ (9b)

In Eq. (9a), ${\bf v}({\bf k})=(1/\hbar)\boldsymbol{\partial}_{\bf k}H({\bf
k})$ with $H({\bf k})=e^{-i{\bf k}\cdot{\bf r}}\mathcal{H}e^{i{\bf k}\cdot{\bf
r}}$, and in Eq. (9b), $g_{s}\approx 2$ is the spin $g$-factor of the
electron. Henceforth, the index ${\bf k}$ will be omitted for brevity.

The $1/\omega$ prefactors in Eq. (7), inherited from Eq. (3), make it singular
at $\omega=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_{\eta\rightarrow 0^{+}}(x-i\eta)^{-1}=1/x+i\pi\delta(x)$. Next, notice
that the reactive part can be rewritten by invoking the Kramers-Krönig
relation

$\sigma_{ab}^{\mathrm{AH}}(\omega_{0},{\bf
q})=-\dfrac{i}{\pi}\mathrm{P}\int_{-\infty}^{\infty}d\omega\,\dfrac{\sigma_{ab}^{\mathrm{H}}(\omega,{\bf
q})}{\omega-\omega_{0}}\ ,$ (10)

while in the absorptive part the factor $1/\omega$ can be replaced with
$1/\omega_{ln}$ thanks to the delta function. Finally, recombine the two parts
to obtain

$\sigma_{ab}(\tilde{\omega},{\bf q})=\dfrac{e^{2}}{\hbar}\sum_{n,l}\int_{\bf
k}\mathcal{E}_{ln}(\tilde{\omega},{\bf q})\mathcal{M}^{ab}_{nl}({\bf q})\,,$
(11)

where we have defined $\tilde{\omega}=\omega+i\eta$ and

$\mathcal{E}_{ln}(\tilde{\omega},{\bf q})=\dfrac{f_{ln}({\bf
q})}{\omega_{ln}({\bf q})}\dfrac{i}{\omega_{ln}({\bf q})-\tilde{\omega}}\,.$
(12)

At zeroth order in ${\bf q}$ Eq. (11) reduces to Eq. (25) of Ref. [46], and at
first order it becomes, in the notation of Eq. (4),

$\sigma_{ab,c}(\tilde{\omega})=\frac{e^{2}}{\hbar}\sum_{n,l}\int_{\bf
k}\left[\mathcal{E}^{,c}_{ln}(\tilde{\omega})\mathcal{M}^{ab}_{nl}({\mathbf{0}})+\mathcal{E}_{ln}(\tilde{\omega},{\mathbf{0}})\mathcal{M}^{ab,c}_{nl}\right]\,.$
(13)

In Sec. 5, we discuss how the equivalence between the Kubo formulas (7) and
(11) at zeroth and first orders in ${\bf 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
$\mathcal{E}$ and $\mathcal{M}$ matrices appearing in Eq. (13). We start in
Sec. 3.2 with the expansion of $\mathcal{E}$, and then devote Secs. 3.3 and
3.4 to the expansion and subsequent manipulations of $\mathcal{M}$, which is
where our treatment differs more substantially from that of previous works.

The terms in the resulting expression for $\sigma_{ab,c}(\tilde{\omega})$ 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
$\mathcal{E}^{,c}_{ln}(\tilde{\omega})$ involves the band velocities ${\bf
v}_{n}=\boldsymbol{\partial}_{\bf k}\omega_{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].

###  3.2 Expansion in ${\bf q}$ of the band-energy terms

When expanding Eq. (12) in powers of ${\bf q}$, the intraband ($l=n$) and
interband ($l\not=n$) parts must be treated separately. To first order in
${\bf q}$, one finds

$\displaystyle\mathcal{E}_{nn}(\tilde{\omega},{\mathbf{0}})$
$\displaystyle=-\dfrac{i}{\tilde{\omega}}f^{\prime}_{n}\,,$ (14a)
$\displaystyle\mathcal{E}_{ln}(\tilde{\omega},{\mathbf{0}})$
$\displaystyle=i\dfrac{f_{ln}}{\omega_{ln}}(\omega_{ln}+\tilde{\omega})Z_{ln}(\tilde{\omega})\,,$
(14b) $\displaystyle\mathcal{E}^{,c}_{nn}(\tilde{\omega})$
$\displaystyle=-\dfrac{i}{\tilde{\omega}^{2}}f^{\prime}_{n}v_{n}^{c}\,,$ (14c)
$\displaystyle\mathcal{E}^{,c}_{ln}(\tilde{\omega})$
$\displaystyle=\dfrac{i}{2}\dfrac{v_{l}^{c}f^{\prime}_{l}+v_{n}^{c}f^{\prime}_{n}}{\omega_{ln}}(\omega_{ln}+\tilde{\omega})Z_{ln}(\tilde{\omega})$
$\displaystyle-
if_{ln}\dfrac{Z_{ln}^{2}(\tilde{\omega})}{\omega_{ln}^{2}}\left[\left(\omega_{ln}+\tilde{\omega}\right)^{2}\left(\omega_{ln}-\tilde{\omega}/2\right)\right](v_{l}^{c}+v_{n}^{c})\,,$
(14d)

where $f^{\prime}_{n}=\partial f_{n}/\partial\omega_{n}$,
$f_{ln}=f_{l}-f_{n}$, $\omega_{ln}=\omega_{l}-\omega_{n}$, and
$Z_{ln}(\tilde{\omega})=1/(\omega_{ln}^{2}-\tilde{\omega}^{2})$. For the
intraband identities, we used

$\frac{f_{nn}({\bf q})}{\omega_{nn}({\bf q})}=f_{n}^{\prime}+{\cal
O}(q^{2})\,,$ (15)

which follows from $\omega_{nn}({\bf q})=(\partial_{a}\omega_{n})q_{a}+{\cal
O}(q^{3})$ and $f_{nn}({\bf q})=(\partial_{a}f_{n})q_{a}+{\cal O}(q^{3})$,
where $\partial_{a}=\partial_{k_{a}}$.

###  3.3 Expansion in ${\bf q}$ of the optical matrix elements

####  3.3.1 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}\rangle\ \rightarrow\ e^{-i\beta_{n}}|u_{n}\rangle$, where $\beta_{n}$
is a real function of ${\bf k}$.

As the $\mathcal{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 $\sigma_{ab,c}(\tilde{\omega})$. 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 $\sigma_{ab,c}(\tilde{\omega})$ in terms of
origin-independent quantities.

The coefficient $\mathcal{M}^{ab}({\mathbf{0}})$ appearing in the first term
of Eq. (13) is trivially gauge invariant, as it involves a single term,

$\mathcal{M}_{nl}^{ab}({\mathbf{0}})=v^{a}_{nl}v^{b}_{ln}\,,$ (16)

which is a product of gauge-covariant velocity matrix elements (clearly, those
matrix elements are also origin independent). Instead, the coefficient
$\mathcal{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 $\langle
u_{n}|v_{a}|\partial_{c}u_{l}\rangle$, because Bloch-state derivatives
transform noncovariantly as $|\boldsymbol{\partial}_{\bf k}u_{n}\rangle\
\rightarrow\ e^{-i\beta_{n}}\left(|\boldsymbol{\partial}_{\bf
k}u_{n}\rangle-i(\boldsymbol{\partial}_{\bf k}\beta_{n})|u_{n}\rangle\right)$.
This can be fixed by writing $|\boldsymbol{\partial}_{\bf k}u_{n}\rangle$ as
$|{\bf D}_{\bf k}u_{n}\rangle-i\boldsymbol{\mathcal{A}}_{n}|u_{n}\rangle$,
where $|{\bf D}_{\bf k}u_{n}\rangle$ is the covariant derivative [39] and
$\boldsymbol{\mathcal{A}}_{n}=\langle u_{n}|i\boldsymbol{\partial}_{\bf
k}u_{n}\rangle$ is the intraband Berry connection. The terms containing Berry
connections cancel out, leaving

$\displaystyle\mathcal{M}^{ab,c}_{nl}$
$\displaystyle=\dfrac{1}{2}\Big{(}v_{nl}^{a}\langle
D_{c}u_{l}|v_{b}|u_{n}\rangle+\langle u_{n}|v_{a}|D_{c}u_{l}\rangle
v_{ln}^{b}-\langle D_{c}u_{n}|v_{a}|u_{l}\rangle v_{ln}^{b}-v_{nl}^{a}\langle
u_{l}|v_{b}|D_{c}u_{n}\rangle\Big{)}$
$\displaystyle+\dfrac{ig_{s}}{2m_{e}}\left(\epsilon_{acd}S_{nl}^{d}v_{ln}^{b}-\epsilon_{bcd}S_{ln}^{d}v_{nl}^{a}\right)\,,$
(17)

where every term is a gauge-invariant product of gauge-covariant matrix
elements, just like in Eq. (16).

####  3.3.2 Degenerate bands

Gyrotropic birefringence and nonreciprocal directional dichroism occur in
antiferromagnetic crystals such as Cr2O3 [3, 4], where the energy bands are
doubly degenerate at every ${\bf k}$ as a result of the combined
$\cal{P}\mathcal{T}$ symmetry [48]. To treat such cases, we introduce
degeneracy indices $\lambda$ and $\nu$ 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,
$\mathcal{M}^{ab}_{nl}=\sum_{\lambda,\nu}\,(I^{a}_{l\lambda,n\nu})^{*}I^{b}_{l\lambda,n\nu}$.
The reasoning leading up to Eq. (17) follows through, provided that the
covariant derivative is generalized as $|{\bf D}_{\bf
k}u_{n\nu}\rangle=|\boldsymbol{\partial}_{\bf
k}u_{n\nu}\rangle+i\sum_{\nu^{\prime}}\,|u_{n\nu^{\prime}}\rangle\boldsymbol{\mathcal{A}}_{n}^{\nu^{\prime}\nu}$,
where $\boldsymbol{\mathcal{A}}_{n}^{\nu^{\prime}\nu}=\langle
u_{n\nu^{\prime}}|i\boldsymbol{\partial}_{\bf k}u_{n\nu}\rangle$ [39]. The
object $|{\bf D}_{\bf k}u_{n\nu}\rangle$ transforms covariantly under
multiband gauge transformations of the form
$|u_{n\nu}\rangle\rightarrow\sum_{\nu^{\prime}}\,|u_{n\nu^{\prime}}\rangle
U_{n}^{\nu^{\prime}\nu}$, where $U_{n}$ is a ${\bf k}$-dependent unitary
matrix in the degeneracy indices. To alleviate the notation, from now on we
will assume nondegenerate bands.

###  3.4 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
$\sigma_{ab,c}(\tilde{\omega})$ 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 [49, 9], this
is achieved by means of identities such as $\langle\phi_{l}|{\bf
v}|\phi_{n}\rangle=i\omega_{ln}\langle\phi_{l}|{\bf r}|\phi_{n}\rangle$.

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

$(\boldsymbol{\partial}_{\bf k}H)|u_{n}\rangle=(\boldsymbol{\partial}_{\bf
k}\varepsilon_{n})|u_{n}\rangle-(H-\varepsilon_{n})|{\bf D}_{\bf
k}u_{n}\rangle\ ,$ (18)

which can be obtained by differentiating
$H|u_{n}\rangle=\varepsilon_{n}|u_{n}\rangle$, and then writing
$|\boldsymbol{\partial}_{\bf k}u_{n}\rangle$ as $|{\bf D}_{\bf
k}u_{n}\rangle-i\boldsymbol{\mathcal{A}}_{n}|u_{n}\rangle$. Contracting with
$\langle u_{l}|$ gives ${\bf v}_{ln}=\delta_{ln}{\bf
v}_{n}+i\omega_{ln}\boldsymbol{{\bf A}}_{ln}$, where ${\bf v}_{n}$ is the band
velocity, and

${\bf A}_{ln}=\langle u_{l}|i\mathbf{D}u_{n}\rangle=(1-\delta_{ln})\langle
u_{l}|i\boldsymbol{\partial}_{\bf k}u_{n}\rangle=\begin{cases}\frac{{\bf
v}_{ln}}{i\omega_{ln}}&\text{if $l\not=n$}\,,\\\ 0&\text{if $l=n$}\end{cases}$
(19)

is the interband Berry connection.

With Eqs. (18) and (19), one can split Eq. (17) for $\mathcal{M}^{ab,c}_{nl}$
as

$\mathcal{M}^{ab,c}_{nl}=i\left(v_{nl}^{a}B_{ln}^{bc}-v_{ln}^{b}B_{nl}^{ac}\right)+\dfrac{\omega_{ln}}{2}\left[(v_{n}^{a}+v_{l}^{a})A^{b}_{ln}A^{c}_{nl}+(v_{n}^{b}+v_{l}^{b})A^{a}_{nl}A^{c}_{ln}\right]\
,$ (20)

where we have defined

$B_{ln}^{bc}=\dfrac{1}{2i\hbar}\left(\langle
D_{b}u_{l}|H-\varepsilon_{l}|D_{c}u_{n}\rangle-\langle
D_{c}u_{l}|H-\varepsilon_{n}|D_{b}u_{n}\rangle\right)-\dfrac{g_{s}}{2m_{e}}\epsilon_{bcd}S_{ln}^{d}\
.$ (21)

The first term in Eq. (20) is molecular for $l\not=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 is
$\bar{d}^{a}=-eA^{a}$, while the intrinsic magnetic-dipole and electric-
quadrupole matrices are related to the antisymmetric and symmetric parts of
$B^{bc}$ as follows,

$\displaystyle\bar{m}^{a}_{ln}$
$\displaystyle=\frac{e}{2}\epsilon_{abc}B^{bc}_{ln}\,,$ (22a)
$\displaystyle\bar{q}^{bc}_{ln}$
$\displaystyle=\frac{ie}{\omega_{ln}}\left(B^{bc}_{ln}+B^{cb}_{ln}\right)\,.$
(22b)

Thus,

$\displaystyle\bar{\bf d}_{nl}$ $\displaystyle=-e\langle u_{n}|i{\bf D}_{\bf
k}u_{l}\rangle\,,$ (23a) $\displaystyle\bar{\bf m}_{ln}$
$\displaystyle=\frac{e}{2i\hbar}\langle{\bf D}_{\bf
k}u_{l}|\times\left(H-\frac{\varepsilon_{l}+\varepsilon_{n}}{2}\right)|{\bf
D}_{\bf k}u_{n}\rangle-\dfrac{eg_{s}}{2m_{e}}{\bf S}_{ln}\,,$ (23b)
$\displaystyle\bar{q}^{bc}_{ln}$ $\displaystyle=-\dfrac{e}{2}\left(\langle
D_{b}u_{l}|D_{c}u_{n}\rangle+\langle D_{c}u_{l}|D_{b}u_{n}\rangle\right)\,,$
(23c)

where $\bar{\bf m}_{ln}$ comprises orbital and spin contributions. By
expanding the covariant derivatives and then setting $l=n$, one finds that
$\bar{\bf d}_{nn}={\mathbf{0}}$ [see Eq. (19)], and that $\bar{\bf m}_{nn}$
and $-\bar{q}^{bc}_{nn}/e$ are respectively the intrinsic magnetic moment
${\bf m}_{n}$ [40] and the quantum metric $g^{bc}_{n}$ [41] of a Bloch
eigenstate.

Using Eq. (19) for ${\bf A}_{ln}$ together with the completeness relation, Eq.
(23) can be recast in a more convenient form for numerical work,

$\displaystyle\bar{\bf d}_{nl}$ $\displaystyle=\begin{cases}ie\frac{{\bf
v}_{nl}}{\omega_{nl}}&\text{if $l\not=n$}\\\ 0&\text{if $l=n$}\end{cases}\,,$
(24a) $\displaystyle\bar{\bf m}^{\text{orb}}_{ln}$
$\displaystyle=\frac{e}{4i}\sum_{p\not=l,n}\left(\frac{1}{\omega_{pl}}+\frac{1}{\omega_{pn}}\right){\bf
v}_{lp}\times{\bf v}_{pn}\,,$ (24b) $\displaystyle\bar{q}^{bc}_{ln}$
$\displaystyle=-\frac{e}{2}\sum_{p\not=l,n}\left[\frac{v^{b}_{lp}v^{c}_{pn}}{\omega_{pl}\omega_{pn}}+(b\leftrightarrow
c)\right]\,.$ (24c)

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. (2) will be established in Sec. 4.

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 $-\bar{q}^{bc}_{ln}/e$ is distinct from the “band-resolved
quantum metric”
$g^{bc}_{ln}=\left(A^{b}_{ln}A^{c}_{nl}+A^{c}_{ln}A^{b}_{nl}\right)/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}_{n}={\rm
Re}\,\langle\partial_{b}u_{n}|\partial_{c}u_{n}\rangle-{\cal A}^{b}_{n}{\cal
A}^{c}_{n}$. Instead, $-\bar{q}^{bc}_{ln}/e$ is gauge covariant for $l\not=n$,
and for $l=n$ it reduces to $g^{bc}_{n}$.

###  3.5 Final expression

We have now gathered all the needed ingredients to evaluate Eq. (13) for
$\sigma_{ab,c}(\tilde{\omega})$, namely the expansion coefficients of
$\mathcal{E}$ in Eq. (14), and those of $\mathcal{M}$ in Eqs. (16) and (20).
In Eqs. (26) and (31) below, we break down the resulting expression into
antisymmetric ($\mathcal{T}$-even) and symmetric ($\mathcal{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}=\omega_{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

$A^{a}_{nl}=-\frac{1}{e}\bar{d}^{a}_{nl}\,,\quad
B^{bc}_{ln}=\frac{1}{e}\bar{m}^{a}_{ln}\epsilon_{abc}+\frac{\omega_{ln}}{2ie}\bar{q}^{bc}_{ln}\,.$
(25)

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 $\cal{P}$, $\mathcal{T}$, or
$\cal{P}\mathcal{T}$ 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.

####  3.5.1 Antisymmetric (time-even) part

The antisymmetric part of $\sigma_{ab,c}(\tilde{\omega})$ takes the form

$\displaystyle\dfrac{\hbar}{e^{2}}\,\sigma_{ab,c}^{\text{A}}(\tilde{\omega})$
$\displaystyle=\tilde{\omega}\sum_{n,l}\int_{\bf
k}Z_{ln}(\tilde{\omega})\bigg{\\{}-f_{ln}{\rm
Im}\left[A^{a}_{nl}B^{bc}_{ln}-(a\leftrightarrow b)\right]$
$\displaystyle+f_{ln}\left[\dfrac{1}{2}\left(v^{a}_{n}+v^{a}_{l}\right){\rm
Im}\left(A^{b}_{nl}A^{c}_{ln}\right)-(a\leftrightarrow b)\right]$
$\displaystyle+f_{ln}\left(3\omega_{ln}^{2}-\tilde{\omega}^{2}\right)Z_{ln}(\tilde{\omega}){\rm
Im}\left(A^{a}_{nl}A^{b}_{ln}\right)\dfrac{1}{2}\left(v^{c}_{n}+v^{c}_{l}\right)$
$\displaystyle-f^{\prime}_{n}\omega_{ln}{\rm
Im}\left(A^{a}_{nl}A^{b}_{ln}\right)v^{c}_{n}\bigg{\\}}$
$\displaystyle+\frac{1}{\tilde{\omega}}\sum_{n}\int_{\bf
k}f^{\prime}_{n}\left(v^{a}_{n}B^{bc}_{nn}-v^{b}_{n}B^{ac}_{nn}\right)\,.$
(26)

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
$\sigma_{ab,c}^{\text{A}}(\tilde{\omega})$ came from differentiating the
$\mathcal{E}$ matrix, that is, from the first term of Eq. (13) for
$\sigma_{ab,c}(\tilde{\omega})$. It went unnoticed in those works that the
other term in that equation – where one differentiates the $\mathcal{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
$\left(\epsilon_{acd}K_{bd}-\epsilon_{bcd}K_{ad}\right)/(e\tilde{\omega})$,
where

$K_{ab}=-\sum_{n}\int_{\bf
k}f^{\prime}_{n}v^{a}_{n}m^{b}_{n}=\sum_{n}\int_{\bf
k}f_{n}\partial_{a}m^{b}_{n}\,.$ (27)

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]

$\rho(\omega,\tau)=\frac{\omega}{2c^{2}\epsilon_{0}}{\rm
Re}\left[\sigma^{\text{A}}_{xy,z}(\omega+i\tau^{-1})\right]\,,$ (28)

where $\epsilon_{0}$ is the vacuum permittivity and $c$ is the speed of light.
To deal with absorption, the positive infinitesimal $\eta$ in
$\tilde{\omega}=\omega+i\eta$ has been reinterpreted heuristically as a
phenomenological scattering rate $\tau^{-1}$ [6, 46, 55]. For frequencies and
scattering rates well below the threshold for interband transitions,
$\omega,\tau^{-1}\ll\omega_{\text{gap}}$, Eq. (26) yields

$\rho(\omega,\tau)=\dfrac{(\omega\tau)^{2}}{1+(\omega\tau)^{2}}a+b\omega^{2}\,.$
(29)

The coefficient $b$ comes from the interband terms which have $\tilde{\omega}$
prefactors, and $a=-(e/c^{2}\epsilon_{0}\hbar)K_{xx}$ (with $K_{xx}=K_{yy}$)
comes from the intraband term with a $1/\tilde{\omega}$ prefactor. In
insulators the coefficient $a$ vanishes, and hence the rotatory power displays
the familiar $\omega^{2}$ dependence at low frequencies [6]; in conductors
that coefficient is nonzero, and one can distinguish two different regimes as
follows,

$\rho(\omega,\tau)\simeq\begin{cases}(\tau^{2}a+b)\omega^{2}&\text{if
$\omega\tau\ll 1$}\ ,\\\ a+b\omega^{2}&\text{if $\omega\tau\gg 1$}\
.\end{cases}$ (30)

In Sec. 6, we will illustrate these low-frequency profiles for a concrete
tight-binding model.

####  3.5.2 Symmetric (time-odd) part

The symmetric part of $\sigma_{ab,c}(\tilde{\omega})$ reads

$\displaystyle\frac{\hbar}{e^{2}}\,\sigma_{ab,c}^{\text{S}}(\tilde{\omega})$
$\displaystyle=i\sum_{n,l}\int_{\bf
k}Z_{ln}(\tilde{\omega})\bigg{\\{}f_{ln}\omega_{ln}{\rm
Re}\left[A^{a}_{nl}B^{bc}_{ln}+(a\leftrightarrow b)\right]$
$\displaystyle+f_{ln}\omega_{ln}\left[\frac{1}{2}\left(v^{a}_{n}+v^{a}_{l}\right){\rm
Re}\left(A^{b}_{nl}A^{c}_{ln}\right)+(a\leftrightarrow b)\right]$
$\displaystyle-f_{ln}\omega^{3}_{ln}Z_{ln}(\tilde{\omega}){\rm
Re}\left(A^{a}_{nl}A^{b}_{ln}\right)\left(v^{c}_{n}+v^{c}_{l}\right)$
$\displaystyle+f^{\prime}_{n}\omega^{2}_{ln}{\rm
Re}\left(A^{a}_{nl}A^{b}_{ln}\right)v^{c}_{n}\bigg{\\}}$
$\displaystyle-\frac{i}{\tilde{\omega}^{2}}\sum_{n}\int_{\bf
k}f^{\prime}_{n}v^{a}_{n}v^{b}_{n}v^{c}_{n}\,.$ (31)

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 $\omega^{2}_{ln}Z_{ln}(\tilde{\omega})$
as $1+\tilde{\omega}^{2}Z_{ln}(\tilde{\omega})$ and noting that ${\rm
Re}\sum_{l}A^{a}_{nl}A^{b}_{ln}$ is the quantum metric
$g^{ab}_{n}=-\bar{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,

$i\sum_{n}\int_{\bf
k}f^{\prime}_{n}g^{ab}_{n}v^{c}_{n}+i\tilde{\omega}^{2}\sum_{n,l}\int_{\bf
k}f^{\prime}_{n}Z_{ln}(\tilde{\omega}){\rm
Re}\left(A^{a}_{nl}A^{b}_{ln}\right)v^{c}_{n}\,.$ (32)

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.

###  3.6 Magnetoelectric and quadrupolar responses

As recognized in Ref. [3], the physical basis for $\sigma_{ab,c}(\omega)$ is
provided by quadrupolar and magnetoelectric couplings. The quadrupolar
couplings are described by a $\mathcal{T}$-odd totally-symmetric tensor
$\gamma_{abc}(\omega)$, and the magnetoelectric ones by a $\mathcal{T}$-odd
traceless tensor $\tilde{\alpha}_{ab}(\omega)$ together with a
$\mathcal{T}$-even tensor $\check{\alpha}_{ab}(\omega)$. Those tensors are
defined as [3, 22]

$\displaystyle\gamma_{abc}(\omega)$
$\displaystyle=\frac{1}{3i}\left[\sigma^{\text{S}}_{ab,c}(\omega)+\sigma^{\text{S}}_{bc,a}(\omega)+\sigma^{\text{S}}_{ca,b}(\omega)\right]\,,$
(33a) $\displaystyle\tilde{\alpha}_{ab}(\omega)$
$\displaystyle=\frac{1}{3i}\sigma^{\text{S}}_{ac,d}(\omega)\epsilon_{cdb}\,,$
(33b) $\displaystyle\check{\alpha}_{ab}(\omega)$
$\displaystyle=\frac{1}{4i}\epsilon_{bcd}\left[\sigma^{\text{A}}_{cd,a}(\omega)-2\sigma^{\text{A}}_{ac,d}(\omega)\right]\,,$
(33c)

so that

$\displaystyle\sigma^{\text{S}}_{ab,c}(\omega)$
$\displaystyle=i\left[\epsilon_{acd}\tilde{\alpha}_{bd}(\omega)+\epsilon_{bcd}\tilde{\alpha}_{ad}(\omega)\right]+i\gamma_{abc}(\omega)\,,$
(34a) $\displaystyle\sigma^{\text{A}}_{ab,c}(\omega)$
$\displaystyle=i\left[\epsilon_{acd}\check{\alpha}_{bd}(\omega)-\epsilon_{bcd}\check{\alpha}_{ad}(\omega)\right]\,.$
(34b)

Inserting Eq. (34) in Eq. (3) for the induced current density and then using
the Maxwell-Faraday equation yields the constitutive relation

$j^{\omega,{\bf q}}_{a}=\left({\boldsymbol{\nabla}}\times{\bf M}^{\omega,{\bf
q}}\right)_{a}+\left[\tilde{\alpha}_{ab}(\omega)-\check{\alpha}_{ab}(\omega)\right]\partial_{t}B^{\omega,{\bf
q}}_{b}+\gamma_{abc}(\omega)\nabla_{b}E^{\omega,{\bf q}}_{c}\,,$ (35)

where ${\bf X}^{\omega,{\bf q}}$ denotes ${\bf X}(\omega,{\bf q})e^{i({\bf
q}\cdot{\bf r}-i\omega t)}$, and $M_{a}^{\omega,{\bf
q}}=\left[\tilde{\alpha}_{ba}(\omega)+\check{\alpha}_{ba}(\omega)\right]E^{\omega,{\bf
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 $\mathcal{T}$-even magnetoelectric tensor $\check{\alpha}_{ab}$ describes
a “kinetic” magnetoelectric effect in gyrotropic conductors [56], while the
$\mathcal{T}$-odd tensor $\tilde{\alpha}_{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,
$\partial_{t}P_{a}^{\omega,{\bf
q}}=\tilde{\alpha}_{ab}(\omega)\partial_{t}B_{b}^{\omega,{\bf q}}$;
integrating the adiabatic current in the quasi-static limit, one obtains the
familiar form $P_{a}=\tilde{\alpha}_{ab}B_{b}$ of the inverse magnetoelectric
effect [1]. The full $\mathcal{T}$-odd magnetoelectric tensor includes an
additional trace (“axion”) piece $(\theta e^{2}/2\pi h)\delta_{ab}$; the axion
angle $\theta$ 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 $\sigma^{\text{S}}_{ab,c}(\omega)$ and
$\sigma^{\text{A}}_{ab,c}(\omega)$, evaluated at $\omega=0$. This has been
done in previous works for selected terms only: in Ref. [22], the Fermi-sea
orbital contribution to $\tilde{\alpha}_{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
$\check{\alpha}_{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 $\mathcal{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.

##  4 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
moments $\bar{d}$, $\bar{m}$, and $\bar{q}$ reduce to those in Eq. (2). We
then show how the formulas for $\sigma^{\text{A}}_{ab,c}(\tilde{\omega})$ and
$\sigma^{\text{S}}_{ab,c}(\tilde{\omega})$ 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_{n{\bf k}}({\bf r})\overset{\boldsymbol{.}}{=}e^{-i{\bf
k}\cdot\boldsymbol{\xi}({\bf r})}\phi_{n}[\boldsymbol{\xi}({\bf r})]\ ,$ (36)

where $\boldsymbol{\xi}({\bf r})={\bf r}-{\bf R}({\bf r})$ is the intracell
coordinate, ${\bf R}({\bf r})$ is the lattice vector that folds the absolute
coordinate ${\bf r}$ into the home unit cell, $\phi_{n}({\bf r})$ is
vanishingly small outside that cell, and $\overset{\boldsymbol{.}}{=}$ 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,

$\boldsymbol{\mathcal{A}}_{n}\overset{\boldsymbol{.}}{=}\int_{\text{cell}}d{\bf
r}\ \phi_{n}^{*}({\bf r})e^{i{\bf k}\cdot{\bf r}}\ i\boldsymbol{\partial}_{\bf
k}\left[e^{-i{\bf k}\cdot{\bf r}}\phi_{n}({\bf r})\right]=\bar{\bf r}_{n}\,,$
(37)

and the covariant derivative of a Bloch state reduces to

${\bf D}_{\bf k}u_{n{\bf k}}({\bf r})\overset{\boldsymbol{.}}{=}-ie^{-i{\bf
k}\cdot{\bf r}}\left({\bf r}-\bar{\bf r}_{n}\right)\phi_{n}({\bf r})\,,$ (38)

for ${\bf r}$ in the home cell. Using this identity in Eq. (23) for $\bar{d}$,
$\bar{m}$, and $\bar{q}$, we recover after some manipulations the expressions
in Eq. (2). In the case of $\bar{m}$, it is necessary to invoke the operator
identity $[v_{a},r_{b}]=[v_{b},r_{a}]$ to rewrite $({\bf r}\times{\bf v}-{\bf
v}\times{\bf r})/2$ as ${\bf r}\times{\bf v}$.

With these relations in hand, we can address the molecular limit of Eqs. (26)
and (31) for $\sigma^{\text{A}}_{ab,c}(\tilde{\omega})$ and
$\sigma^{\text{S}}_{ab,c}(\tilde{\omega})$. 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 ${\bf k}$, we can set $\int_{\bf
k}\overset{\boldsymbol{.}}{=}1/V_{c}$ in those terms ($V_{c}$ is the cell
volume) to find

$\displaystyle V_{c}\sigma^{\text{A}}_{ab,c}$
$\displaystyle\overset{\boldsymbol{.}}{=}\bar{G}^{\prime}_{ad}\epsilon_{dbc}+\frac{\tilde{\omega}}{2}\bar{a}_{abc}-(a\leftrightarrow
b)\,,$ (39a) $\displaystyle iV_{c}\sigma^{\text{S}}_{ab,c}$
$\displaystyle\overset{\boldsymbol{.}}{=}-\bar{G}_{ad}\epsilon_{dbc}+\frac{\tilde{\omega}}{2}\bar{a}^{\prime}_{abc}+(a\leftrightarrow
b)\,,$ (39b)

where we dropped the $\tilde{\omega}$ dependence for brevity, and defined the
(extensive) molecular tensors

$\displaystyle\bar{G}_{ab}$
$\displaystyle=\frac{1}{\hbar}\sum_{n,l}f_{nl}\omega_{ln}Z_{ln}{\rm
Re}\left(\bar{d}_{nl}^{a}\bar{m}_{ln}^{b}\right)\ ,$ (40a)
$\displaystyle\bar{G}_{ab}^{\prime}$
$\displaystyle=-\frac{1}{\hbar}\sum_{n,l}f_{nl}\tilde{\omega}Z_{ln}{\rm
Im}\left(\bar{d}_{nl}^{a}\bar{m}_{ln}^{b}\right)\ ,$ (40b)
$\displaystyle\bar{a}_{abc}$
$\displaystyle=\frac{1}{\hbar}\sum_{n,l}f_{nl}\omega_{ln}Z_{ln}{\rm
Re}\left(\bar{d}_{nl}^{a}\bar{q}_{ln}^{bc}\right)\ ,$ (40c)
$\displaystyle\bar{a}_{abc}^{\prime}$
$\displaystyle=-\frac{1}{\hbar}\sum_{n,l}f_{nl}\dfrac{\omega_{ln}^{2}}{\tilde{\omega}}Z_{ln}{\rm
Im}\left(\bar{d}_{nl}^{a}\bar{q}_{ln}^{bc}\right)\ .$ (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,
$\bar{G}$, $\bar{G}^{\prime}$ and $\bar{a}$ become single-particle versions of
the multipolar susceptibility tensors $G$, $G^{\prime}$ 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 with $\bar{d}$,
$\bar{m}$, and $\bar{q}$. It is not immediately clear that the same is true
for $\bar{a}^{\prime}$, since Eq. (40d) contains a factor of
$\omega^{2}_{ln}/\tilde{\omega}$ in place of the $\tilde{\omega}$ factor
appearing in Eq. (2.84) of Ref. [7] for $a^{\prime}$. However, those factors
are interchangeable in the expression for $a^{\prime}$, as can be seen in the
manner described around around Eqs. (2.75–2.78) of Ref. [7]. Thus,
$(\bar{G},\bar{G}^{\prime},\bar{a},\bar{a}^{\prime})$ are origin-independent
versions of the molecular tensors $(G,G^{\prime},a,a^{\prime})$ entering the
standard multipole theory.

Next, let us consider the propagation of light inside our idealized molecular
crystal. For a given propagation direction $\hat{\bf n}$, we define
(intensive) optical-activity and gyrotropic-birefringence tensors as
$\beta^{\text{A}}_{ab}(\hat{\bf n})=-\sigma^{\text{A}}_{ab,c}\hat{n}_{c}$ and
$\beta^{\text{S}}_{ab}(\hat{\bf n})=i\sigma^{\text{S}}_{ab,c}\hat{n}_{c}$,
respectively. Using Eq. (39), we obtain Eqs. (5.8) and (5.9) of Ref. [7] for
those tensors, but with $(\bar{G},\bar{G}^{\prime},\bar{a},\bar{a}^{\prime})$
in place of $(G,G^{\prime},a,a^{\prime})$. 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, $(\bar{d},\bar{m},\bar{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 ${\bf 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 magnetic-dipole and
electric-quadrupole terms [6, 7, 8, 9].

##  5 Sum rules

In Sec. 3.1, we wrote two alternative Kubo formulas for
$\sigma_{ab}(\omega,{\bf q})$, namely Eqs. (7) and (11). The former displays
apparent $1/\omega$ divergences at $\omega=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 ${\bf q}$, and relate those identities to the oscillator-
and rotatory-strength sum rules.

###  5.1 Equivalence between the two forms of the Kubo formula

Let us denote as $(ie^{2}/\omega)\Delta_{ab}(\mathbf{q})$ the difference
between the reactive parts of the Kubo formulas (7) and (11). Writing
$1/[x(a-x)]$ as $(1/a)[1/x+1/(a-x)]$, we find

$\Delta_{ab}({\bf
q})=\delta_{ab}\dfrac{N}{m_{e}}+\dfrac{1}{\hbar}\sum_{n,l}\int_{\bf
k}\dfrac{f_{ln}({\bf q})}{\omega_{ln}({\bf q})}M_{nl}^{ab}({\bf q})\ ,$ (41)

and using

$\displaystyle\omega_{ln}({\bf q})$ $\displaystyle=-\omega_{nl}(-{\bf q})\,,$
(42a) $\displaystyle f_{ln}({\bf q})$ $\displaystyle=-f_{nl}(-{\bf q})\,,$
(42b) $\displaystyle M^{ab}_{nl}({\bf q})$
$\displaystyle=\left[M^{ab}_{ln}(-{\bf q})\right]^{*}\,,$ (42c) $\displaystyle
M^{ab}_{nl}({\bf q})$ $\displaystyle=\left[M^{ba}_{nl}({\bf q})\right]^{*}\,,$
(42d)

we obtain

$\displaystyle{\rm Re}\,\Delta_{ab}({\bf q})$ $\displaystyle={\rm
Re}\,\Delta_{ba}({\bf q})={\rm Re}\,\Delta_{ab}(-{\bf q})\,,$ (43a)
$\displaystyle{\rm Im}\,\Delta_{ab}({\bf q})$ $\displaystyle=-{\rm
Im}\,\Delta_{ba}({\bf q})=-{\rm Im}\,\Delta_{ab}(-{\bf q})\,.$ (43b)

For the two Kubo formulas to be equivalent, $\Delta_{ab}({\bf 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 $\Delta_{ab}({\bf
q})$ at zeroth and first order in ${\bf q}$, we expand it as

$\Delta_{ab}({\bf
q})=\Delta_{ab}({\mathbf{0}})+\Delta_{ab,c}q_{c}+\mathcal{O}(q^{2})\ .$ (44)

Let us start with the zeroth-order term in the expansion. Writing the electron
density $N$ in the first term of Eq. (41) as $\sum_{n}\int_{\bf 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

$\Delta_{ab}({\mathbf{0}})=\sum_{n}\int_{\bf
k}f_{n}\left[\dfrac{\delta_{ab}}{m_{e}}+2\sum_{l\neq n}\dfrac{{\rm
Re}\left(v_{nl}^{a}v_{ln}^{b}\right)}{\varepsilon_{n}-\varepsilon_{l}}-\dfrac{1}{\hbar^{2}}\dfrac{\partial^{2}\varepsilon_{n}}{\partial
k_{a}\partial k_{b}}\right]=0\ .$ (45)

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 $\omega=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 $\Delta_{ab}({\bf q})$ with the $\omega\rightarrow 0$ limit of
$\Pi_{ab}(\omega,{\bf q})=i\omega\sigma_{ab}(\omega,{\bf q})$ [see Eq. (3)],
evaluated using the Kubo formula (7). This gives $-e^{2}\Delta_{ab}({\bf
q})=\Pi_{ab}(0,{\bf q})$, and so
$-e^{2}\Delta_{ab,c}=\left.\partial_{q_{c}}\Pi_{ab}(0,{\bf q})\right|_{{\bf
q}={\mathbf{0}}}$. From the analysis of this quantity in Ref. [23] (see Sec.
III.A.1 of its Supplemental Material), we conclude that

$\Delta_{ab,c}=-\dfrac{i}{\hbar}\epsilon_{abc}\sum_{n}\int_{\bf k}f_{n}\,{\bf
v}_{n}\cdot\boldsymbol{\Omega}_{n}=0\ ,$ (46)

where $\boldsymbol{\Omega}_{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 $\sigma_{ab,c}(\omega)$
that contains a term diverging as $1/\omega$. 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 $\Delta_{ab}(\mathbf{q})$ dictates the high-frequency
behavior of the optical conductivity as follows. Suppose there is a frequency
$\omega_{\text{max}}$ above which the system does not absorb [6]; setting
$\omega\gg\omega_{\text{max}}$ in Eqs. (11) and (12) and comparing with Eq.
(41), we can deduce that

$\sigma_{ab}(\omega\gg\omega_{\mathrm{max}},{\bf
q})=\delta_{ab}\dfrac{ie^{2}N}{\omega
m_{e}}-\frac{ie^{2}}{\omega}\Delta_{ab}({\bf q})\,.$ (47)

Thus, at high frequencies the optical conductivity reduces to the diamagnetic
term; and since that term is independent of ${\bf q}$, we conclude that
$\sigma_{ab}(\omega,{\mathbf{0}})$ decays as $1/\omega$.

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 $\Delta_{ab,c}$, as we will now show.

###  5.2 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
$\int_{0}^{\infty}f(\omega)d\omega$ as $\langle f(\omega)\rangle$ and using
Eq. (5) yields

$\left<\sigma^{\text{H}}_{ab}(\omega,{\bf q})\right>=\left<{\rm
Re}\,\sigma^{\text{S}}_{ab}(\omega,{\bf q})\right>+i\left<{\rm
Im}\,\sigma^{\text{A}}_{ab}(\omega,{\bf q})\right>\,.$ (48)

To evaluate this quantity, we begin by taking the Hermitian part of Eq. (11)
for $\eta\rightarrow 0^{+}$,

$\sigma^{\text{H}}_{ab}(\omega,{\bf q})=-\frac{\pi
e^{2}}{\hbar}\sum_{n,l}\int_{\bf k}\frac{f_{ln}({\bf q})}{\omega_{ln}({\bf
q})}\mathcal{M}^{ab}_{nl}({\bf q})\delta\left[\omega-\omega_{ln}({\bf
q})\right]\,.$ (49)

Making the substitution

$f_{ln}({\bf q})=f_{l}({\bf k}+{\bf q}/2)\left[1-f_{n}({\bf k}-{\bf
q}/2)\right]-f_{n}({\bf k}-{\bf q}/2)\left[1-f_{l}({\bf k}+{\bf
q}/2)\right]\,,$ (50)

and noting that at zero temperature only the second term contributes to Eq.
(49) when $\omega>0$ and ${\bf q}\approx{\mathbf{0}}$, we obtain

$\left<{\rm Re}\,\sigma^{\text{S}}_{ab}(\omega,{\bf q})\right>+i\left<{\rm
Im}\,\sigma^{\text{A}}_{ab}(\omega,{\bf q})\right>=R_{ab}({\bf
q})+i\,I_{ab}({\bf q})\,,$ (51)

where we have defined

$\displaystyle R_{ab}({\bf q})$ $\displaystyle=\dfrac{\pi
e^{2}}{\hbar}\sum_{n,l}\int_{\bf k}f_{n}({\bf k}-{\bf q}/2)\left[1-f_{l}({\bf
k}+{\bf q}/2)\right]\dfrac{{\rm Re}\left[\mathcal{M}_{nl}^{ab}({\bf
q})\right]}{\omega_{ln}({\bf q})}\,,$ (52a) $\displaystyle I_{ab}({\bf q})$
$\displaystyle=\dfrac{\pi e^{2}}{\hbar}\sum_{n,l}\int_{\bf k}f_{n}({\bf
k}-{\bf q}/2)\left[1-f_{l}({\bf k}+{\bf q}/2)\right]\dfrac{{\rm
Im}\left[\mathcal{M}_{nl}^{ab}({\bf q})\right]}{\omega_{ln}({\bf q})}\,.$
(52b)

Let us split $R_{ab}({\bf q})$ and $I_{ab}({\bf q})$ in Eq. (51) into even and
odd parts in ${\bf q}$. Using the identities in Eq. (42), one finds that the
even part of $R_{ab}({\bf q})$ plus the odd part of $iI_{ab}({\bf q})$ is
proportional to the second term in Eq. (41) for $\Delta_{ab}({\bf q})$.
Therefore,

$\displaystyle\left<{\rm Re}\,\sigma^{\text{S}}_{ab}(\omega,{\bf q})\right>$
$\displaystyle=\frac{\pi e^{2}}{2}\left[\delta_{ab}N/m_{e}-{\rm
Re}\,\Delta_{ab}({\bf q})\right]+\frac{1}{2}\left[R_{ab}({\bf q})-R_{ab}(-{\bf
q})\right]\,,$ (53a) $\displaystyle\left<{\rm
Im}\,\sigma^{\text{A}}_{ab}(\omega,{\bf q})\right>$ $\displaystyle=-\frac{\pi
e^{2}}{2}{\rm Im}\,\Delta_{ab}({\bf q})+\frac{1}{2}\left[I_{ab}({\bf
q})+I_{ab}(-{\bf q})\right]\,,$ (53b)

where we keep track of the vanishing quantity $\Delta_{ab}({\bf q})$.

The expansion of Eq. (53) in powers of ${\bf q}$ generates a series of sum
rules. Since, according to Eq. (43), the terms ${\rm Re}\,\Delta_{ab}$ and
${\rm Im}\,\Delta_{ab}$ only contribute (formally) at even and odd orders in
${\bf q}$, respectively, and since the reverse is true for the second terms in
Eqs. (53a) and (53b), we obtain

$\displaystyle\left<{\rm
Re}\,\sigma^{\text{S}}_{ab}(\omega,{\mathbf{0}})\right>$
$\displaystyle=\frac{\pi
e^{2}}{2}\left[\delta_{ab}N/m_{e}-\Delta_{ab}({\mathbf{0}})\right]\,,$ (54a)
$\displaystyle\left<{\rm
Im}\,\sigma^{\text{A}}_{ab}(\omega,{\mathbf{0}})\right>$
$\displaystyle=I_{ab}({\mathbf{0}})\,,$ (54b) $\displaystyle\left<{\rm
Re}\,\sigma^{\text{S}}_{ab,c}(\omega)\right>$ $\displaystyle=R_{ab,c}\,,$
(54c) $\displaystyle\left<{\rm Im}\,\sigma^{\text{A}}_{ab,c}(\omega)\right>$
$\displaystyle=-\frac{\pi e^{2}}{2}\Delta_{ab,c}\,,$ (54d)

to linear order in ${\bf q}$. Below, we consider each of these identities in
turn.

Once we set $\Delta_{ab}({\mathbf{0}})=0$ in accordance with Eq. (45), Eq.
(54a) becomes the oscillator-strength sum rule [1]

$\left<{\rm
Re}\,\sigma^{\text{S}}_{ab}(\omega,{\mathbf{0}})\right>=\frac{\omega^{2}_{\text{p}}}{8}\delta_{ab}\,,$
(55)

where $\omega_{\text{p}}=(4\pi e^{2}N/m_{e})^{1/2}$ is the plasma frequency.
As already mentioned, for tight-binding models the diamagnetic term in Eq. (7)
changes form while Eq. (11) remains unchanged, which leads to a modified
oscillator-strength sum rule [47].

Equation (54b) is the rotatory-strength sum rule for magnetic circular
dichroism. At ${\bf q}={\mathbf{0}}$, the intraband part of Eq. (52b) vanishes
because $\mathcal{M}_{nn}^{ab}({\mathbf{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
${\bf q}$.

Finally, by setting $\Delta_{ab,c}=0$ in Eq. (54d) in accordance with Eq.
(46), we arrive at the rotatory-strength sum rule for natural circular
dichroism,

$\left<{\rm Im}\,\sigma^{\text{A}}_{ab,c}(\omega)\right>=0\,.$ (56)

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
$\Delta_{ab}({\bf q})$ – at zero and first order in ${\bf q}$, respectively.
More generally, the expansion of Eq. (53) in powers of ${\bf q}$ yields at
each order two optical sum rules, one of which relies on the vanishing of
$\Delta_{ab}({\bf q})$ at that order.

##  6 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.

Figure 1: (a) The tight-binding model of Eq. (57). The crystallographic
vectors are $\mathbf{a}_{1}=(\sqrt{3}a,0,0)$,
$\mathbf{a}_{2}=(\sqrt{3}a/2,3a/2,0)$, and $\mathbf{a}_{3}=(0,0,c)$, with $a$
the distance between nearest-neighbor sites on the same layer, and $c$ the
interlayer distance. The on-site energies and hoppings are indicated
schematically. As shown on the right, the complex intralayer hoppings from
$+\Delta$ sites to $-\Delta$ sites have amplitude $it_{1}$; the reverse
hoppings (not shown) have amplitude $-it_{1}$. (b) Band structure of the model
for the Hamiltonian parameters given in the main text.

As a simple model of a bulk crystal with nonzero $\sigma^{\text{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 $\sigma^{\text{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

$\mathcal{H}=\Delta\sum_{i}\xi_{i}c_{i}^{\dagger}c_{i}+it_{1}\sum_{\langle
i,j\rangle}\xi_{j}c_{i}^{\dagger}c_{j}+\dfrac{i\lambda_{1}}{a}\sum_{\langle
i,j\rangle}c_{i}^{\dagger}\left({\boldsymbol{\sigma}}\cdot{\boldsymbol{\delta}}_{ij}\right)c_{j}+\dfrac{i\lambda_{2}}{a}\sum_{[i,j]}c_{i}^{\dagger}\left({\boldsymbol{\sigma}}\cdot{\bf
d}_{ij}\right)c_{j}\ .$ (57)

The first term is a staggered on-site potential, with $\xi_{i}=\pm 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, ${\boldsymbol{\sigma}}$ is the vector of Pauli
matrices and ${\boldsymbol{\delta}}_{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 ${\bf 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$, $\Delta=0.5$, $\lambda_{1}=-0.06$, and $\lambda_{2}=0.05$.

Exploiting the translational symmetry of the crystal, we replace the site
indices $\\{i\\}$ with $\\{{\bf R}i\\}$, where the lattice vector ${\bf R}$
labels the cell, and $i$ is now an intracell site index. The Hamiltonian
matrix elements are denoted by $\mathcal{H}_{ij}({\bf
R})=\langle\phi_{{\mathbf{0}}i}|\mathcal{H}|\phi_{{\bf R}j}\rangle$, where
$\phi_{{\bf R}j}({\bf r})=\varphi_{j}({\bf r}-{\bf
R}-{\boldsymbol{\tau}}_{j})$ is a basis orbital centered at ${\bf
R}-{\boldsymbol{\tau}}_{j}$ [39]. The tight-binding Hamiltonian in ${\bf k}$
space is constructed as

$H^{\bf k}_{ij}=\sum_{\bf R}e^{i{\bf k}\cdot({\bf
R}+{\boldsymbol{\tau}}_{j}-{\boldsymbol{\tau}}_{i})}\mathcal{H}_{ij}({\bf
R})\,,$ (58)

leading to the eigenvalue equation $H_{\bf k}\cdot C_{n{\bf
k}}=\varepsilon_{n{\bf k}}C_{n{\bf k}}$.

Figure 2:  Numerical results for the nonzero components of $\sigma_{ab,c}$ for
the model of Eq. (57) and Fig. 1, with the two lowest bands treated as
occupied. The results are plotted as a function of frequency up to
$\hbar\omega=0.3$, which is well below the threshold for interband absorption
($\varepsilon_{\text{gap}}\approx 0.53$). Filled circles: extrapolation from
finite-size crystallites. Solid lines: bulk crystal. The tensor
$\sigma_{ab,c}$ has been divided by $e^{2}/\hbar$ to make it dimensionless.

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 $\sigma^{\text{A}}_{ab,c}$ and
$\sigma^{\text{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]

${\bf v}_{nl}({\bf k})=\frac{1}{\hbar}C_{n{\bf
k}}^{\dagger}\cdot\left(\boldsymbol{\partial}_{\bf k}H_{\bf k}\right)\cdot
C_{l{\bf k}}\,.$ (59)

As the system is insulating, the only nonzero contributions to $\sigma_{ab,c}$
come from the Fermi-sea 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 $\eta=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],
$\sigma_{yz,x}^{\mathrm{A}}$, $\sigma_{xy,z}^{\mathrm{A}}$,
$\sigma_{xx,y}^{\mathrm{S}}$ and $\sigma_{xz,y}^{\mathrm{S}}$, only the first
three are actually nonzero when the Fermi level $\varepsilon_{\text{F}}$ lies
in the gap. Converged results, obtained by sampling the Brillouin zone on a
uniform mesh of $50\times 50\times 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 $\sigma^{\text{A}}_{ab,c}$ and $\sigma^{\text{S}}_{ab,c}$ from Eq.
(39) under open boundary conditions. Calculations are performed for samples
with $L+1$ unit cells in each crystallographic direction, with $L$ ranging
from 4 to 12\. The results are extrapolated to $L\rightarrow\infty$ by fitting
them to the function

$f(L)=f_{0}+f_{1}/L+f_{2}/L^{2}+f_{3}/L^{3}\,,$ (60)

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 $\sigma^{\text{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
containing $\bar{m}$ are denoted as $\bar{\text{M}}1$, and those containing
$\bar{q}$ are denoted as $\bar{\text{E}}2$. For the bulk crystal (right
panel), there are, in addition to $\bar{\text{M}}1$ and $\bar{\text{E}}2$
contributions from the first term in Eq. (26), band-dispersion contributions
from the second and third terms, which are denoted as ${\bf v}$.

A comparison between the two panels of Fig. 3 reveals that the
$\bar{\text{M}}1$ and $\bar{\text{E}}2$ contributions are different for the
extrapolated crystallites and for the bulk crystal, with the difference being
exactly compensated by the ${\bf 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 [70,
69]. For a two-dimensional insulator with a single valence band, that term
reads $-(e/2\hbar)\int_{\bf k}\varepsilon_{\bf
k}(\partial_{x}{\mathcal{A}}_{y}-\partial_{y}{\mathcal{A}}_{x})$ or, after
integrating by parts, $(e/2)\int_{\bf
k}(v_{x}{\mathcal{A}}_{y}-v_{y}{\mathcal{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 [70,
69].

To conclude, let us illustrate the different low-frequency behaviors of the
rotatory power in insulating and conducting states of our model. We evaluate
$\rho(\omega,\tau)$ for the bulk model from Eqs. (26) and (28), setting
$\hbar/\tau=2\times 10^{-3}$. The frequency range is chosen as
$0\leq\hbar\omega\leq 10^{-2}$, and the calculation is carried out at zero
temperature for $\varepsilon_{\text{F}}=0.0$ (insulator) and
$\varepsilon_{\text{F}}=1.0$ (metal). In both cases, a uniform mesh of
$100\times 100\times 100$ $k$ points 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 $\rho(\omega)$ profiles, plotted in Fig. 4 as solid lines,
display the behavior dictated by Eq. (29). The insulator displays a simple
$\rho\propto\omega^{2}$ decay with a negligible influence from the scattering
time $\tau$. Instead, for the metal one can distinguish two different
parabolic regimes (dashed lines) delimited by $\omega\tau\sim 1$, in
accordance with Eq. (30). This distinctive low-frequency profile of the
optical rotatory dispersion in conducting crystals awaits experimental
verification.

Figure 3: Decomposition of the $\sigma_{xy,z}^{\mathrm{A}}$ curves in Fig. 2
into three types of origin-independent contributions: intrinsic magnetic
dipole ($\bar{\text{M}}1$), intrisic electric quadrupole ($\bar{\text{E}}2$),
and band dispersive (${\bf 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.

##  7 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 physically-transparent 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 $\sigma_{ab,c}(\omega)$
at first order in ${\bf q}$, higher-order responses can in principle be
treated in a similar manner.

Figure 4: Low-frequency optical rotatory dispersion of the model of Eq. (57)
and Fig. 1 in an insulating state ($\varepsilon_{\text{F}}=0.0$), and in a
metallic state ($\varepsilon_{\text{F}}=1.0$). Numerical results based on Eqs.
(26) and (28) are shown as solid lines. The dashed lines are guides to the
eye, indicating the distinct parabolic behaviors in the metallic state for
$\omega\tau\ll 1$ and for $\omega\tau\gg 1$, as per Eq. (30), and the faint
vertical line indicates the crossover frequency $\omega=1/\tau$. The rotatory
power has units of radians per length.

For a crystal composed of nonoverlapping units, our formula for
$\sigma_{ab,c}(\omega)$ 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
$(\bar{d},\bar{m},\bar{q})$ given by Eq. (2). Away from the molecular limit,
$\sigma_{ab,c}(\omega)$ 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
$\sigma_{ab,c}(\omega)$ 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 ${\bf q}$; this allowed us to eliminate spurious noncovariant terms, and to
isolate the physically relevant matrix elements $\bar{d}$, $\bar{m}$, and
$\bar{q}$. 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 ${\bf 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 ${\bf k}$-space Berry curvature.

Our work opens up new prospects for realistic ab initio calculations of
spatially-dispersive optical responses in crystals, and of static
magnetoelectric and quadrupolar responses as well. An implementation based on
the sum-over-states formulas for $(\bar{d},\bar{m},\bar{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 ${\bf 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 orbital-current operator whose matrix elements in
the Bloch-eigenstate basis are proportional to
$\sum_{l}\big{(}m^{z,\text{orb}}_{ml}v^{a}_{ln}+v^{a}_{ml}m^{z,\text{orb}}_{ln}\big{)}$.
Here, ${\bf m}^{\text{orb}}_{ln}$ is the bulk generalization of Eq. (1b),
given by the same expression as in Eq. (23b) for $\bar{\bf
m}^{\text{orb}}_{ln}$, but with the covariant derivatives therein replaced by
ordinary derivatives. The two definitions are related by

$\bar{\bf m}^{\text{orb}}_{ln}={\bf
m}^{\text{orb}}_{ln}-\frac{e}{4i}\omega_{ln}\left(\boldsymbol{\mathcal{A}}_{l}+\boldsymbol{\mathcal{A}}_{n}\right)\times{\bf
A}_{ln}\,,$ (61)

and therefore they agree for $l=n$ only. For $l\not=n$, the two terms on the
right-hand side of Eq. (61) are not separately gauge covariant, and the lack
of gauge covariance of ${\bf m}^{\text{orb}}_{ln}$ makes the orbital Hall
conductivity gauge dependent. This suggest that one should generally use
$\bar{\bf m}^{\text{orb}}_{ln}$ instead of ${\bf m}^{\text{orb}}_{ln}$ when
evaluating the orbital Hall conductivity. In other words, one is allowed to
work with ${\bf m}^{\text{orb}}_{ln}$ only in the parallel-transport gauge,
where
$\boldsymbol{\mathcal{A}}_{l}=\boldsymbol{\mathcal{A}}_{n}={\mathbf{0}}$.

## Acknowledgements

This work was supported by the IKUR Strategy under the collaboration agreement
between the Ikerbasque Foundation and the Material Physics Center on behalf of
the Department of Education of the Basque Government, and by Grant No.
PID2021-129035NB-I00 funded by MCIN/AEI/10.13039/501100011033. The authors
wish to thank Cheol-Hwan Park, Stepan Tsirkin, Xiaoming Wang, Fernando de
Juan, Raffaele Resta, and Massimiliano Stengel for valuable discussions and
comments on the manuscript.

## References

  * [1] L. D. Landau and E. M. Lifshitz,  _Electrodynamics of Continuous Media_ ,  Elsevier, Amsterdam, Netherlands,  ISBN 9780080302751 (1984).
  * [2] D. B. Melrose and R. C. McPhedran,  _Electromagnetic Processes in Dispersive Media_ ,  Cambridge University Press, Cambridge, UK,  ISBN 9780511600036 (1991),  10.1017/CBO9780511600036.
  * [3] R. M. Hornreich and S. Shtrikman,  _Theory of gyrotropic birefringence_ ,  Phys. Rev. 171, 1065 (1968),  10.1103/PhysRev.171.1065.
  * [4] B. Graham and R. E. Raab,  _Magnetic effects in antiferromagnetic crystals in the electric quadrupole-magnetic dipole approximation_ ,  Phil. Mag. B 66, 269 (1992),  10.1080/13642819208224589.
  * [5] E. U. Condon,  _Theories of optical rotatory power_ ,  Rev. Mod. Phys. 9, 432 (1937),  10.1103/RevModPhys.9.432.
  * [6] L. D. Barron,  _Molecular Light Scattering and Optical Activity_ ,  Cambridge University Press, Cambridge, UK,  ISBN 9780511535468 (2004),  10.1017/CBO9780511535468.
  * [7] R. E. Raab and O. L. De Lange,  _Multipole Theory in Electromagnetism: Classical, quantum, and symmetry aspects, with applications_ ,  Oxford University Press, Oxford, UK,  ISBN 9780198567271 (2005),  10.1093/acprof:oso/9780198567271.001.0001.
  * [8] P. Norman, K. Ruud and T. Saue,  _Principles and Practices of Molecular Properties: Theory, modeling, and simulations_ ,  John Wiley & Sons, Hoboken, USA,  ISBN 9780470725627 (2018),  10.1002/9781118794821.
  * [9] T. B. Pedersen and A. E. Hansen,  _Ab initio calculation and display of the rotatory strength tensor in the random phase approximation. Method and model studies_ ,  Chem. Phys. Lett. 246, 1 (1995),  10.1016/0009-2614(95)01036-9.
  * [10] P. L. Polavarapu,  _Ab initio molecular optical rotations and absolute configurations_ ,  Mol. Phys. 91, 551 (1997),  10.1080/002689797171436.
  * [11] D. Varsano, L. A. Espinosa-Lean, X. Andrade, M. A. L. Marque, R. di Felice and A. Rubio,  _Towards a gauge invariant method for molecular chiroptical properties in TDDFT_ ,  Phys. Chem. Chem. Phys. 11, 4481 (2009),  10.1039/b903200b.
  * [12] K. Yabana and G. F. Bertsch,  _Application of the time-dependent local density approximation to optical activity_ ,  Phys. Rev. A 60, 1271 (1999),  10.1103/PhysRevA.60.1271.
  * [13] P. L. Polavarapu,  _Optical rotation: Recent advances in determining the absolute configuration_ ,  Chirality 14, 768 (2002),  10.1002/chir.10145.
  * [14] J. Autschbach, L. Nitsch-Velasquez and M. Rudolph,  _Time-Dependent Density Functional Response Theory for Electronic Chiroptical Properties of Chiral Molecules_ ,  Top. Curr. Chem. 298, 1 (2011),  10.1007/128_2010_72.
  * [15] J. Mattiat and S. Luber,  _Recent Progress in the Simulation of Chiral Systems with Real Time Propagation Methods_ ,  Helv. Chim. Acta 104, e2100154 (2021),  10.1002/hlca.202100154.
  * [16] W. Moffitt,  _Optical Rotatory Dispersion of Helical Polymers_ ,  J. Chem. Phys. 25, 467 (1956),  10.1063/1.1742946.
  * [17] T. Kato, I. Tsujikawa and T. Murao,  _Optical Activity Caused by Exciton Dispersion_ ,  J. Phys. Soc. Jpn 34, 763 (1973),  10.1143/JPSJ.34.763.
  * [18] V. M. Agranovich and V. L. Ginzburg,  _Crystal Optics with Spatial Dispersion, and Excitons_ ,  Springer, Berlin, Heidelberg, Germany,  ISBN 9783662024089 (1984),  10.1007/978-3-662-02406-5.
  * [19] K. Natori,  _Band theory of the optical activity of crystals_ ,  J. Phys. Soc. Japan 39, 1013 (1975),  10.1143/JPSJ.39.1013.
  * [20] H. Zhong, Z. H. Levine, D. C. Allan and J. W. Wilkins,  _Band-theoretic calculations of the optical-activity tensor of $\alpha$-quartz and trigonal Se_,  Phys. Rev. B 48, 1384 (1993),  10.1103/PhysRevB.48.1384.
  * [21] L. Jönsson, Z. H. Levine and J. W. Wilkins,  _Large local-field corrections in optical rotatory power of quartz and selenium_ ,  Phys. Rev. Lett. 76, 1372 (1996),  10.1103/PhysRevLett.76.1372.
  * [22] A. Malashevich and I. Souza,  _Band theory of spatial dispersion in magnetoelectrics_ ,  Phys. Rev. B 82, 245118 (2010),  10.1103/PhysRevB.82.245118.
  * [23] S. Zhong, J. E. Moore and I. Souza,  _Gyrotropic Magnetic Effect and the Magnetic Moment on the Fermi Surface_ ,  Phys. Rev. Lett. 116, 077201 (2016),  10.1103/PhysRevLett.116.077201.
  * [24] J. Ma and D. A. Pesin,  _Chiral magnetic effect and natural optical activity in metals with or without Weyl points_ ,  Phys. Rev. B 92, 235205 (2015),  10.1103/PhysRevB.92.235205.
  * [25] Y. Gao and D. Xiao,  _Nonreciprocal directional dichroism induced by the quantum metric dipole_ ,  Phys. Rev. Lett. 122, 227402 (2019),  10.1103/PhysRevLett.122.227402.
  * [26] A. H. Duff and J. E. Sipe,  _Magnetoelectric polarizability and optical activity: Spin and frequency dependence_ ,  Phys. Rev. B 106, 085413 (2022),  10.1103/PhysRevB.106.085413.
  * [27] F. Hidalgo, A. Sánchez-Castillo and C. Noguez,  _Efficient first-principles method for calculating the circular dichroism of nanostructures_ ,  Phys. Rev. B 79, 075438 (2009),  10.1103/PhysRevB.79.075438.
  * [28] S. S. Tsirkin, P. A. Puente and I. Souza,  _Gyrotropic effects in trigonal tellurium studied from first principles_ ,  Phys. Rev. B 97, 035158 (2018),  10.1103/PhysRevB.97.035158.
  * [29] M. Rérat and B. Kirtman,  _First-Principles Calculation of the Optical Rotatory Power of Periodic Systems: Application on $\alpha$-Quartz, Tartaric Acid Crystal, and Chiral $(n,m)$-Carbon Nanotubes_,  J. Chem. Theory Comput. 17, 4063 (2021),  10.1021/acs.jctc.1c00243.
  * [30] M. O. Yokosuk, H.-S. Kim, K. D. Hughey, J. Kim, A. V. Stier, K. R. O. Neal, J. Yang, S. A. Crooker, K. Haule, S.-W. Cheong, D. Vanderbilt and J. L. Musfeldt,  _Nonreciprocal directional dichroism of a chiral magnet in the visible range_ ,  npj Quantum Materials 5, 20 (2020),  10.1038/s41535-020-0224-6.
  * [31] K. Park, M. O. Yokosuk, M. Goryca, J. J. Yang, S. A. Crooker, S.-W. Cheong, K. Haule, D. Vanderbilt, H.-S. Kim and J. L. Musfeldt,  _Nonreciprocal directional dichroism at telecom wavelengths_ ,  npj Quantum Materials 7, 38 (2022),  10.1038/s41535-022-00438-6.
  * [32] T. Arima,  _Magneto-electric optics in non-centrosymmetric ferromagnets_ ,  J. Phys. Condens. Matter 20, 434211 (2008),  10.1088/0953-8984/20/43/434211.
  * [33] C.-J. Kim, A. Sánchez-Castillo, Z. Ziegler, Y. Ogawa, C. Noguez and J. Park,  _Chiral atomically thin films_ ,  Nature Nanotech. 11, 520 (2016),  10.1038/nnano.2016.3.
  * [34] E. S. Morell, L. Chico and L. Brey,  _Twisting Dirac fermions: circular dichroism in bilayer graphene_ ,  2D Materials 4, 035015 (2017),  10.1088/2053-1583/aa7eb6.
  * [35] T. Stauber, T. Low and G. Gómez-Santos,  _Chiral response of twisted bilayer graphene_ ,  Phys. Rev. Lett. 120, 046801 (2018),  10.1103/PhysRevLett.120.046801.
  * [36] K. Chang, Z. Zheng, J. E. Sipe and J. L. Cheng,  _Theory of optical activity in doped systems with application to twisted bilayer graphene_ ,  Phys. Rev. B 106, 245405 (2022),  10.1103/PhysRevB.106.245405.
  * [37] Z. Han, F. Wang, J. Sun, X. Wang and Z. Tang,  _Recent advances in ultrathin chiral metasurfaces by twisted stacking_ ,  Adv. Mater. 35, 2206141 (2023),  10.1002/adma.202206141.
  * [38] E. I. Blount,  _Formalisms of band theory_ ,  in _Solid state physics_ , Academic Press, Cambridge, USA,  ISBN 9780126077131 (1962),  10.1016/S0081-1947(08)60459-2.
  * [39] D. Vanderbilt,  _Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators_ ,  Cambridge University Press, Cambridge, UK,  ISBN 9781316662205 (2018),  10.1017/9781316662205 .
  * [40] D. Xiao, M.-C. Chang and Q. Niu,  _Berry phase effects on electronic properties_ ,  Rev. Mod. Phys. 82, 1959 (2010),  10.1103/RevModPhys.82.1959.
  * [41] J. P. Provost and G. Vallee,  _Riemannian Structure on Manifolds of Quantum States_ ,  Commun. Math. Phys. 76, 289 (1980),  10.1007/BF02193559.
  * [42] K. Rzazewski and R. W. Boyd,  _Equivalence of interaction Hamiltonians in the electric dipole approximation_ ,  J. Mod. Opt. 51, 1137 (2004),  10.1080/09500340408230412.
  * [43] W. A. Harrison,  _Solid State Theory_ ,  Dover Publications, New York, USA,  ISBN 9780486639482 (1980).
  * [44] M. Dressel and G. Grüner,  _Electrodynamics of Solids: Optical properties of electrons in matter_ ,  Cambridge University Press, Cambridge, UK,  ISBN 9780511606168 (2002),  10.1017/CBO9780511606168.
  * [45] H. Bruus and K. Flensberg,  _Many-body quantum theory in condensed matter physics_ ,  Oxford University Press, Oxford, UK,  ISBN 9780198566355 (2004).
  * [46] P. B. Allen,  _Electron transport_ ,  in _Conceptual Foundations of Materials_ ,  Elsevier, Amsterdam, Netherlands,  ISBN 9780444509765 (2006),  10.1016/S1572-0934(06)02006-3.
  * [47] M. Graf and P. Vogl,  _Electromagnetic fields and dielectric response in empirical tight-binding theory_ ,  Phys. Rev. B 51, 4940 (1995),  10.1103/PhysRevB.51.4940.
  * [48] Y. Yafet,  _$g$ factors and spin-lattice relaxation of conduction electrons_,  in _Solid State Physics_ ,  Academic Press, Cambridge, USA,  ISBN 9780126077148 (1963),  10.1016/S0081-1947(08)60259-3.
  * [49] J. Snir and J. Schellman,  _Optical Activity of Oriented Helices. Quadrupole Contributions_ ,  J. Phys. Chem. 77, 2653 (1973),  10.1021/j100632a012.
  * [50] Q. Ma, A. G. Grushin and K. S. Burch,  _Topology and geometry under the nonlinear electromagnetic spotlight_ ,  Nat. Mater. 20, 1601 (2021),  10.1038/s41563-021-00992-7.
  * [51] Y. Gao, Y. Zhang and D. Xiao,  _Tunable Layer Circular Photogalvanic Effect in Twisted Bilayers_ ,  Phys. Rev. Lett. 124, 077401 (2020),  10.1103/PhysRevLett.124.077401.
  * [52] H. Watanabe and Y. Yanase,  _Chiral Photocurrent in Parity-Violating Magnet and Enhanced Response in Topological Antiferromagnet_ ,  Phys. Rev. X 11, 011001 (2021),  10.1103/PhysRevX.11.011001.
  * [53] J. Ahn, G.-Y. Guo, N. Nagaosa and A. Vishwanath,  _Riemannian geometry of resonant optical responses_ ,  Nature Phys. 18, 290 (2022),  10.1038/s41567-021-01465-z.
  * [54] C. Şahin, J. Rou, J. Ma and D. A. Pesin,  _Pancharatnam-Berry phase and kinetic magnetoelectric effect in trigonal tellurium_ ,  Phys. Rev. B 97, 205206 (2018),  10.1103/PhysRevB.97.205206.
  * [55] D. J. Passos, G. B. Ventura, J. M. V. P. Lopes, J. M. B. L. d. Santos and N. M. R. Peres,  _Nonlinear optical responses of crystalline systems: Results from a velocity gauge analysis_ ,  Phys. Rev. B 97, 235446 (2018),  10.1103/PhysRevB.97.235446.
  * [56] L. S. Levitov, Y. V. Nazarov and G. M. Eliashberg,  _Magnetoelectric effects in conductors with mirror isomer symmetry_ ,  Sov. Phys. JETP 61, 133 (1985).
  * [57] A. Malashevich, I. Souza, S. Coh and D. Vanderbilt,  _Theory of orbital magnetoelectric response_ ,  New J. Phys. 12(5), 053032 (2010),  10.1088/1367-2630/12/5/053032.
  * [58] T. Yoda, T. Yokoyama and S. Murakami,  _Current-induced orbital and spin magnetizations in crystals with helical structure_ ,  Sci. Rep. 5, 12024 (2015),  10.1038/srep12024.
  * [59] M. F. Lapa and T. L. Hughes,  _Semiclassical wave packet dynamics in nonuniform electric fields_ ,  Phys. Rev. B 99, 121111 (2019),  10.1103/PhysRevB.99.121111.
  * [60] C. Xiao, H. Liu, J. Zhao, S. A. Yang and Q. Niu,  _Thermoelectric generation of orbital magnetization in metals_ ,  Phys. Rev. B 103, 045401 (2021),  10.1103/PhysRevB.103.045401.
  * [61] C. Xiao, Y. Ren and B. Xiong,  _Adiabatically induced orbital magnetization_ ,  Phys. Rev. B 103, 115432 (2021),  10.1103/PhysRevB.103.115432.
  * [62] W. Kohn,  _Theory of the insulating state_ ,  Phys. Rev. 133, A171 (1964),  10.1103/PhysRev.133.A171.
  * [63] G. Weinreich,  _Solids: Elementary Theory for Advanced Students_ ,  John Wiley & Sons, New York, USA,  ISBN 9780471925903 (1965).
  * [64] J. E. Sipe and E. Ghahramani,  _Nonlinear optical response of semiconductors in the independent-particle approximation_ ,  Phys. Rev. B 48, 11705 (1993),  10.1103/PhysRevB.48.11705.
  * [65] I. Souza and D. Vanderbilt,  _Dichroic $f$-sum rule and the orbital magnetization of crystals_,  Phys. Rev. B 77, 054438 (2008),  10.1103/PhysRevB.77.054438.
  * [66] W. Yao, D. Xiao and Q. Niu,  _Valley-dependent optoelectronics from inversion symmetry breaking_ ,  Phys. Rev. B 77, 235406 (2008),  10.1103/PhysRevB.77.235406.
  * [67] R. E. Newnham,  _Properties of materials: anisotropy, symmetry, structure_ ,  Oxford University Press, Oxford, UK,  ISBN 9780198520764 (2005),  10.1093/oso/9780198520757.001.0001.
  * [68] S. V. Gallego, J. Etxebarria, L. Elcoro, E. S. Tasci and J. M. Perez-Mato,  _Automatic calculation of symmetry-adapted tensors in magnetic and non-magnetic materials: a new tool of the Bilbao Crystallographic Server_ ,  Acta Crystallographica Section A 75, 438 (2019),  10.1107/S2053273319001748.
  * [69] D. Ceresoli, T. Thonhauser, D. Vanderbilt and R. Resta,  _Orbital magnetization in crystalline solids: Multi-band insulators, Chern insulators, and metals_ ,  Phys. Rev. B 74, 024408 (2006),  10.1103/PhysRevB.74.024408.
  * [70] T. Thonhauser, D. Ceresoli, D. Vanderbilt and R. Resta,  _Orbital magnetization in periodic insulators_ ,  Phys. Rev. Lett. 95, 137205 (2005),  10.1103/PhysRevLett.95.137205.
  * [71] X. Wang and Y. Yan,  _Optical activity of solids from first principles_ ,  Phys. Rev. B 107, 045201 (2023),  10.1103/PhysRevB.107.045201.
  * [72] X. Wang, J. R. Yates, I. Souza and D. Vanderbilt,  _Ab initio calculation of the anomalous Hall conductivity by Wannier interpolation_ ,  Phys. Rev. B 74, 195118 (2006),  10.1103/PhysRevB.74.195118.
  * [73] S. Bhowal and G. Vignale,  _Orbital Hall effect as an alternative to valley Hall effect in gapped graphene_ ,  Phys. Rev. B 103, 195309 (2021),  10.1103/PhysRevB.103.195309.
  * [74] T. P. Cysne, S. Bhowal, G. Vignale and T. G. Rappoport,  _Orbital Hall effect in bilayer transition metal dichalcogenides: From the intra-atomic approximation to the Bloch states orbital magnetic moment approach_ ,  Phys. Rev. B 105, 195421 (2022),  10.1103/PhysRevB.105.195421.
  * [75] A. Pezo, D. García Ovalle and A. Manchon,  _Orbital Hall effect in crystals: Interatomic versus intra-atomic contributions_ ,  Phys. Rev. B 106, 104414 (2022),  10.1103/PhysRevB.106.104414.