# Natural optical activity from density-functional perturbation theory

Asier Zabalo  Institut de Ciència de Materials de Barcelona (ICMAB-CSIC),
Campus UAB, 08193 Bellaterra, Spain  Massimiliano Stengel  Institut de Ciència
de Materials de Barcelona (ICMAB-CSIC), Campus UAB, 08193 Bellaterra, Spain
ICREA-Institució Catalana de Recerca i Estudis Avançats, 08010 Barcelona,
Spain

###### Abstract

We present an accurate and computationally efficient first-principles
methodology to calculate the natural optical activity. Our approach is based
on the long-wave density-functional perturbation theory and includes self-
consistent field (SCF) terms naturally in the formalism, which are found to be
of crucial importance. The final result is expressed exclusively in terms of
response functions to uniform field perturbations and avoids troublesome
summations over empty states. Our strategy is validated by computing the
natural optical activity tensor in representative chiral crystals (trigonal
Se, $\alpha$-HgS and $\alpha$-SiO2) and molecules (C4H4O2), finding excellent
agreement with experiment and previous theoretical calculations.

Natural optical activity (NOA) refers to the first-order spatial dispersion of
the macroscopic dielectric tensor [1]. Empirically, it manifests as optical
rotation (OR), which is a property of certain structures to rotate the plane
of the polarization of light that travels through them [2, 3]; at difference
with the Faraday effect, NOA is reciprocal and doesn’t require magnetism to be
present [4]. It was first measured in quartz crystals back in 1811 by Arago,
and historically, most of the studied optically active materials turned out to
be chiral. In fact, chirality is a sufficient but not necessary condition for
NOA to be present, as optically active achiral systems also exist [5]. Since
its discovery, natural optical activity has been attracting increasing
research interest, and reliable experimental measurements now exist for many
materials, both in molecular [6, 7, 8, 9, 10] and crystalline form [11, 12,
13, 14, 15, 16, 17].

Parallel to the experiments, there have been considerable advances in the
theoretical understanding of optical rotation as well [2, 11, 18, 19, 20]. Ab-
initio methods like Hartee-Fock (HF), [9] coupled-cluster (CC) [21] and
density functional theory (DFT) [6, 22, 7] have recently become popular in the
context of NOA. While most of the available literature is about small
molecules, notable attempts at calculating optical activity in solids do
exist. It is worth mentioning, for example, the pioneering works by Zhong,
Levine, Allan and Wilkins, [19, 20] based on a numerical long-wavelength
expansion of the electromagnetic response function. Later, Malashevich and
Souza [10] and Pozo and Souza [24] derived analytical expressions for the NOA,
thus reviving the interest in the field; their formalism has been implemented
very recently within an ab initio context [25]. The agreement between theory
and experiment achieved in these works is quite good, e.g., for trigonal Se
[24, 10], $\alpha$-quartz [19, 20] and trigonal Te [26].

In spite of the remarkable progress, however, a systematic, first-principles-
based and computationally efficient methodology to compute the NOA has not
been established yet. The first issue concerns the treatment of the self-
consistent (SCF) fields. These were accounted for in Ref. [27] and found to be
of crucial importance, but the numerical differentiations with respect to the
wave vector ${\bf q}$ that were used therein have limited a widespread
application of their method. The existing analytical expressions [10, 24] for
the NOA are, in principle, better suited to an _ab initio_ implementation
[25], but the SCF contributions are systematically neglected therein. Another
disadvantage with the existing techniques lies in that they require cumbersome
sums over empty states; this introduces an additional potential source of
error, as the convergence with respect to the number of bands tends to be
slow. There are additional technical subtleties that have not been considered
in the context of the NOA, for example regarding the correct treatment of the
current-density response in presence of nonlocal pseudopotentials [28]. It is
unquestionable that the current limitations rule out the study of many systems
of outstanding interest (e.g., electrotoroidic compounds [29, 30]), which are
hardly accessible to the currently available schemes.

Here we present, within the framework of first-principles long-wave density
functional perturbation theory (DFPT), a method to calculate the natural
optical activity that overcomes the aforementioned limitations and is equally
valid for molecules and extended solids. Building on Ref. [7], we express the
natural optical activity tensor as the first-order spatial dispersion (i.e.,
derivative with respect to the wave vector $\mathbf{q}$) of the macroscopic
dielectric function. Crucially, the capabilities of the recently implemented
long-wave module [32] of abinit [33, 34] allow for an efficient calculation by
combining response functions that are already available in the code (e.g.,
$\mathbf{k}$-derivatives, electric and orbital magnetic field perturbations).
This way, summations over excited states are entirely avoided, and the effect
of local fields is automatically included without the need for an ad-hoc
treatment. We validate our methodology by computing the NOA tensor for well
known chiral structures, including trigonal crystals (Se, $\alpha$-HgS and
$\alpha$-SiO2) and the C4H4O2 molecule. Our numerical results show fast
convergence with respect to the main computational parameters, and are in
excellent agreement with experiment and earlier theoretical calculations.

Our starting point is the double Fourier transform in frequency $\omega$ and
wave vector $\mathbf{q}$ of the permittivity function,
$\epsilon_{\alpha\beta}(\omega,\mathbf{q})$. By expanding
$\epsilon_{\alpha\beta}(\omega,\mathbf{q})$ in powers of the wave vector
$\mathbf{q}$, around $\mathbf{q=0}$, we obtain

$\epsilon_{\alpha\beta}(\omega,\mathbf{q})=\epsilon_{\alpha\beta}(\omega,\mathbf{q=0})+iq_{\gamma}\eta_{\alpha\beta\gamma}(\omega)+\dots,$
(1)

where $\eta_{\alpha\beta\gamma}(\omega)$ is the natural optical activity
tensor [1]. (From now on, we adopt Einstein summation conventions for the
Cartesian indices $\alpha\beta\gamma$.) In absence of dissipation (i.e., in
the transparent regime), $\epsilon_{\alpha\beta}(\omega,\mathbf{q})$ is a
$3\times 3$ Hermitian matrix, which at ${\bf q=0}$ becomes real symmetric in
crystals with time-reversal (TR) symmetry. The frequency-dependent natural
optical activity tensor is then also real and satisfies
$\eta_{\alpha\beta\gamma}(\omega)=-\eta_{\beta\alpha\gamma}(\omega)$, which
means that only 9 of the 27 components of $\eta_{\alpha\beta\gamma}$ are
independent. As a consequence, $\eta_{\alpha\beta\gamma}$ is often rearranged
into the second-rank gyration or gyrotropic tensor, $g_{\alpha\beta}$, [1]

$g_{\alpha\beta}(\omega)=\frac{1}{2}\epsilon_{\gamma\delta\alpha}\eta_{\gamma\delta\beta}(\omega),$
(2)

where $\epsilon_{\gamma\alpha\delta}$ is the Levi-Civita symbol. Assuming a
crystal structure with the point group 32 (trigonal Se, $\alpha$-HgS and
$\alpha$-SiO2 belong to this crystal class), and considering that the optical
axis is oriented along the $z$ Cartesian direction [19],

$\boldsymbol{g}(\omega)=\begin{pmatrix}g_{11}(\omega)&0&0\\\
0&g_{11}(\omega)&0\\\ 0&0&g_{33}(\omega)\end{pmatrix},$ (3)

where $g_{11}=\eta_{231}$ and $g_{33}=\eta_{123}$. The optical rotatory power
$\rho$ is then given by [19]

$\rho(\omega)=\frac{\omega^{2}}{2c^{2}}g_{33}(\omega),$ (4)

where $c$ is the speed of light. In this work, we shall focus on the
$\omega\rightarrow 0$ limit, where the components of both $\boldsymbol{g}$ and
$\boldsymbol{\eta}$ tend to a finite constant,

$\eta_{\alpha\beta\gamma}=\eta_{\alpha\beta\gamma}(\omega\rightarrow 0),\qquad
g_{\alpha\beta}=g_{\alpha\beta}(\omega\rightarrow 0).$ (5)

At leading order in the frequency, this yields a rotatory power of

$\rho(\omega)\simeq(\hbar\omega)^{2}\bar{\rho},\qquad\bar{\rho}=\frac{g_{33}}{2(\hbar
c)^{2}},$ (6)

where $\hbar$ is the reduced Planck constant. The constant $\bar{\rho}$ is
usually expressed in the units of deg/[mm (eV)2] and can be directly compared
to experimental measurements.

To make further progress, we shall express the dielectric function in the low-
frequency limit as a second derivative of the ground state energy with respect
to two spatially modulated electric fields ($\boldsymbol{\mathcal{E}}$) [35]

$\epsilon_{\alpha\beta}(\mathbf{q})=\delta_{\alpha\beta}-\frac{4\pi}{\Omega}E_{\mathbf{q}}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}.$
(7)

This allows us to write the natural optical activity tensor as the first
derivative of $\epsilon_{\alpha\beta}(\mathbf{q})$ with respect to
$q_{\gamma}$,

$\eta_{\alpha\beta\gamma}=-\frac{4\pi}{\Omega}\text{Im}E_{\gamma}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}},\qquad
E_{\gamma}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}=\frac{\partial
E_{\mathbf{q}}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}}{\partial
q_{\gamma}}\Big{|}_{\bf q=0},$ (8)

where $\Omega$ is the volume of the unit cell. By virtue of the “$2n+1$”
theorem [7], $E_{\gamma}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}$ can be
written in terms of uniform-field response functions, which are already
available in public first-principles packages like abinit. More specifically,
we find

$E_{\gamma}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}=E_{\text{elst},\gamma}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}+2s\int_{\text{BZ}}[d^{3}k]E_{\mathbf{k},\gamma}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}},$
(9)

where $s=2$ is the spin multiplicity, and the shorthand notation
$[d^{3}k]={\Omega}/{(2\pi)^{3}}d^{3}k$ is used for the Brillouin-zone (BZ)
integral. (We assume that the system under study is a time-reversal (TR)
symmetric insulator.) The electrostatic (elst) term is defined as

$E_{\text{elst},\gamma}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}=\int_{\Omega}\int
n^{\mathcal{E}_{\alpha}}(\mathbf{r})K_{\gamma}(\mathbf{r,r^{\prime}})n^{\mathcal{E}_{\beta}}d^{3}rd^{3}r^{\prime},$
(10)

where $n^{\mathcal{E}_{\beta}}$ is the first-order charge density response to
$\mathcal{E}_{\beta}$, and $K_{\gamma}(\mathbf{r,r^{\prime}})$ is the first
$\mathbf{q}$-derivative of the Hartree exchange and correlation (Hxc) kernel.
The wave function term of Eq. (9), in turn, can be written as

$\begin{split}E_{\mathbf{k},\gamma}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}=&\mathcal{X}_{\mathbf{k}}^{\mathcal{E}_{\alpha}k_{\gamma}\mathcal{E}_{\beta}}+\mathcal{Y}_{\mathbf{k}}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}k_{\gamma}}+\mathcal{Y}_{\mathbf{k}}^{k_{\gamma}\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}\\\
&+\mathcal{W}^{\alpha,\beta\gamma}_{\bf
k}+\left(\mathcal{W}^{\beta,\alpha\gamma}_{\bf k}\right)^{*},\end{split}$ (11)

We shall explain Eq. (11) term by term in the following.

For three generic perturbations, $\lambda_{1}$, $\lambda_{2}$ and
$\lambda_{3}$, the calligraphic symbols in the first line are defined as

$\mathcal{X}_{\mathbf{k}}^{\lambda_{1}\lambda_{2}\lambda_{3}}=\sum_{m}\bra*{u_{m\mathbf{k}}^{\lambda_{1}}}\hat{\mathcal{H}}^{\lambda_{2}}_{\mathbf{k}}\ket*{u_{m\mathbf{k}}^{\lambda_{3}}}$
(12)

and

$\mathcal{Y}_{\mathbf{k}}^{\lambda_{1}\lambda_{2}\lambda_{3}}=-\sum_{mn}\bra*{u_{m\mathbf{k}}^{\lambda_{1}}}\ket*{u_{n\mathbf{k}}^{\lambda_{3}}}\bra*{u_{n\mathbf{k}}^{(0)}}\hat{\mathcal{H}}_{\mathbf{k}}^{\lambda_{2}}\ket*{u_{m\mathbf{k}}^{(0)}}.$
(13)

(The band indices $m,n$ run over the occupied states only.) Here,
$\ket*{u_{m\mathbf{k}}^{\lambda}}$ are the first-order wave functions and the
first-order calligraphic Hamiltonian is given by
$\hat{\mathcal{H}_{\mathbf{k}}^{\lambda}}=\hat{H}^{\lambda}_{\mathbf{k}}+\hat{V}^{\lambda}$,
where $\hat{H}^{\lambda}_{\mathbf{k}}$ is the external perturbation and
$\hat{V}^{\lambda}$ is the self-consistent field (SCF) potential response.
Note that
$\hat{\mathcal{H}}_{\mathbf{k}}^{k_{\gamma}}=\hat{{H}}_{\mathbf{k}}^{k_{\gamma}}$
as there is no SCF contribution to the derivative in $\mathbf{k}$-space, and
$\hat{\mathcal{H}}_{\mathbf{k}}^{\mathcal{E}_{\alpha}}=\hat{V}^{\mathcal{E}_{\alpha}}$
in the above equations since the “external potential” is a purely cross-gap
operator in the electric-field case [7].

The third line is defined as

$\mathcal{W}^{\alpha,\beta\gamma}_{\bf
k}=\sum_{m}i\bra*{u_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}\ket*{u_{m\mathbf{k},\gamma}^{A_{\beta}}},$
(14)

where $\ket*{u_{m\mathbf{k},\gamma}^{A_{\beta}}}$ indicates the wave function
response to an electromagnetic vector potential at first order in the
modulation vector ${\bf q}$. (See Sec. IV of Ref. [36] for more details.) We
can write $\mathcal{W}$ as a sum of two contributions that are, respectively,
symmetric ($\mathcal{S}_{\mathbf{k}}^{\alpha,\beta\gamma}$) and an
antisymmetric ($\mathcal{A}_{\mathbf{k}}^{\alpha,\beta\gamma}$) with respect
to $\beta\leftrightarrow\gamma$ exchange,

$\mathcal{W}^{\alpha,\beta\gamma}_{\bf
k}=\mathcal{S}_{\mathbf{k}}^{\alpha,\beta\gamma}+\mathcal{A}_{\mathbf{k}}^{\alpha,\beta\gamma}.$
(15)

These objects are given by

$\mathcal{S}_{\mathbf{k}}^{\alpha,\beta\gamma}=\frac{i}{2}\sum_{m}\bra*{u_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}\ket*{\partial_{\beta\gamma}^{2}u_{m\mathbf{k}}^{(0)}}$
(16)

and

$\mathcal{A}_{\mathbf{k}}^{\alpha,\beta\gamma}=\frac{1}{2}\sum_{m}\epsilon_{\delta\beta\gamma}\bra*{u_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}\ket*{u_{m\mathbf{k}}^{B_{\delta}}}.$
(17)

In Eq. (16), $\partial^{2}_{\beta\gamma}$ represents a second derivative in
$\mathbf{k}$-space. The
$\ket*{\partial_{\beta\gamma}^{2}u_{m\mathbf{k}}^{(0)}}$ functions in
$\mathcal{S}$ are the well known $d^{2}/dk_{\beta}dk_{\gamma}$ wave functions
[7, 8, 38]; whereas in Eq. (17), $\ket*{u_{m\mathbf{k}}^{B_{\delta}}}$ is the
wave function response to a uniform orbital magnetic field, $B_{\delta}$, as
defined in Refs. [9, 40].

For finite systems, the above theory nicely recovers the established formulas
that are used in quantum chemistry calculations (more details can be found in
Sec. VI of Ref. [36]). Our formulation, however, presents many crucial
advantages. First, Eq. (9) has been derived within a DFPT framework, and hence
avoids the cumbersome summations over unoccupied states that are required by
other methods. Second, all contributions to Eq. (11) are individually
independent of the choice of the origin, and equally valid for both molecules
and extended crystals; this implies that our formulas are free of cancellation
errors due to incomplete basis sets. Third, all the aforementioned terms are
independent of the choice of the wave function gauge by construction, as they
are all expressed as parametric derivatives (with respect to ${\bf q}$) of
multiband gauge-invariant quantities. Fourth, the treatment of the current-
density response in presence of nonlocal pseudopotentials complies with the
prescriptions of Ref. [28]. Finally, and most importantly, SCF terms naturally
appear in our formalism, both directly in $E_{\text{elst}}$ and $\mathcal{Y}$
(both terms vanish if local fields are neglected), and indirectly in the other
terms via the first-order wave functions
$\ket*{u_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}$ (see Sec. V. of Ref. [36]).

A natural question to ask at this point is whether Eq. (11) is unique, or
whether there other combinations of the same ingredients that yield the same
result. Two inequivalent definitions of
$E_{\mathbf{k},\gamma}^{\mathcal{E}_{\alpha}^{*}\mathcal{E}_{\beta}}$ can, at
most, differ by a vanishing Brillouin-zone integral; so the question boils
down to asking whether we can combine the individual pieces in Eq. (11) in
such a way that the result is the total $\mathbf{k}$-derivative of some
function $f({\bf k})$. An obvious choice for $f({\bf k})$ consists in
identifying it with the $\mathbf{k}$-derivative of the macroscopic dielectric
tensor. Indeed, by applying the “$2n+1$” theorem to the stationary expression
[41, 42, 35] for $E_{\mathbf{k},{\bf
q=0}}^{\mathcal{E}_{\alpha}^{*}\mathcal{E}_{\beta}}$, we find

$\begin{split}\frac{\partial
E_{\mathbf{k,q}}^{\mathcal{E}_{\alpha}^{*}\mathcal{E}_{\beta}}}{\partial
k_{\gamma}}\Big{|}_{\mathbf{q=0}}=&\mathcal{X}_{\mathbf{k}}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}k_{\gamma}}+\mathcal{X}_{\mathbf{k}}^{k_{\gamma}\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}+\mathcal{X}_{\mathbf{k}}^{\mathcal{E}_{\alpha}k_{\gamma}\mathcal{E}_{\beta}}\\\
&+\mathcal{Y}_{\mathbf{k}}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}k_{\gamma}}+\mathcal{Y}_{\mathbf{k}}^{\mathcal{E}_{\alpha}k_{\gamma}\mathcal{E}_{\beta}}+\mathcal{Y}_{\mathbf{k}}^{k_{\gamma}\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}\\\
&+2\mathcal{S}_{\mathbf{k}}^{\alpha,\beta\gamma}+2\Big{(}\mathcal{S}_{\mathbf{k}}^{\beta,\alpha\gamma}\Big{)}^{*}.\end{split}$
(18)

Then, by subtracting the latter expression from Eq. (11), we obtain another
equally valid formula for the NOA,

$\begin{split}\Big{[}E_{\mathbf{k},\gamma}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}\Big{]}^{\prime}=&-\Big{(}\mathcal{X}_{\mathbf{k}}^{\mathcal{E}_{\alpha}\mathcal{E}_{\beta}k_{\gamma}}+\mathcal{X}_{\mathbf{k}}^{k_{\gamma}\mathcal{E}_{\alpha}\mathcal{E}_{\beta}}+\mathcal{Y}_{\mathbf{k}}^{\mathcal{E}_{\alpha}k_{\gamma}\mathcal{E}_{\beta}}\Big{)}\\\
&-\mathcal{W}^{\alpha,\gamma\beta}_{\bf
k}-\left(\mathcal{W}^{\beta,\gamma\alpha}_{\bf k}\right)^{*}.\end{split}$ (19)

Numerical tests confirm the consistency of Eq. (11) and (19) to a very high
degree of accuracy. We therefore conclude that Eq. (9) is not unique; on the
contrary, there are infinite possible definition of the gyrotropy tensor,
differing from our Eq. (9) by a dimensionless constant times Eq. (18).

This arbitrariness can be regarded a direct consequence of the
_electromagnetic_ (EM) gauge freedom. Indeed, the last lines in both Eq. (11)
and Eq. (19) have the physical meaning of Berry curvatures in the parameter
space spanned by a uniform magnetic field (${\bf B}$) and an electric field.
Such curvatures are, as we said, insensitive to the choice of the coordinate
origin and the wave function gauge. This result was achieved by expressing the
${\bf B}$-field response function in a cell-periodic form, consistent with the
density-operator theory of Essin et al. [9]. Notwithstanding these undeniable
advantages, the aforementioned Berry curvatures retain an inherent dependence
on the EM-gauge [43]. More specifically, the symbol
$\mathcal{W}^{\alpha,\beta\gamma}$ is expressed in a Landau gauge where the
$\beta$ component of the ${\bf A}$-field increases linearly along $\gamma$; so
when going from Eq. (11) to Eq. (19) we have essentially switched between two
Landau gauges in the last term, and collected the leftovers in the form of
$\mathcal{X}$ and $\mathcal{Y}$. (It is, of course, possible to define a third
variant of Eq. (11), where the contribution of $\mathcal{S}$ cancels out, at
the expense of having a slightly longer list of $\mathcal{X}$\- and
$\mathcal{Y}$-symbols.) Ideally, one would like to exploit this freedom to
obtain a physically intuitive separation between well-defined (and possibly
individually measurable) physical effects; whether such a choice exists is an
interesting open question, which we shall defer to a later work.

Table 1: Calculated independent components of the gyration tensor (in bohr) and the optical rotatory power $\bar{\rho}$ defined in Eq. (6) (in deg/[mm (eV)2] units). Values in brackets are computed neglecting the SCF terms. | $g_{11}$ | $g_{33}$ | $\bar{\rho}$
---|---|---|---
Se | -1.307 (-1.547) | -1.913 (-0.458) | -74.5 (-17.8)
$\alpha$-HgS | 0.775 (0.554) | -1.861 (-1.274) | -72.5 (-49.6)
$\alpha$-SiO2 | -0.071 (-0.001) | 0.125 (0.019) | 4.9 (0.7)

Our first principles calculations are performed with the open-source abinit
[34, 33] package. (Details of the computational parameters are provided in
Sec. I of Ref. [36].) Overall, our approach displays a remarkably fast
convergence with respect to the main computational parameters (plane-wave
energy cutoff and number of $\mathbf{k}$ points, see Ref. [36], Sec. III). In
Table 1 we show the converged numerical values for the independent components
of the gyration tensor and the optical rotatory power in our test set of
trigonal crystals: trigonal Se, $\alpha$-HgS and $\alpha$-SiO2 (numerical
values in brackets are obtained neglecting SCF terms). Our results are in
fairly good agreement with literature values, even if a scissor operator was
applied in Ref. [27] to correct the LDA band gap. (Trigonal Se is an
interesting exception, in that Ref. [27] reports an opposite sign to ours for
the non-SCF value of the $g_{33}$ component; although the reason for this
discrepancy is unclear, we remain confident in the accuracy of our results, as
other values nicely agree with ours in both magnitude and sign.) Overall, our
results confirm the crucial importance of local-field SCF contributions,
consistent with the conclusions of Ref. [27].

Table 2: Comparison between LDA and GGA for the independent components of the gyration tensor for Se, $\alpha$-HgS and $\alpha$-SiO2, for different structures. In the Structure column, “exp” refers to the experimental structure, while Se (LDA) means that the structure was relaxed with LDA, for example. Structure |  $g_{11}$ (bohr) |  $g_{33}$ (bohr)
---|---|---
| LDA | GGA | LDA | GGA
Se(exp) | -1.306 | -1.301 | -1.910 | -1.329
Se(GGA) | -1.408 | -1.431 | -1.802 | -1.216
$\alpha$-HgS (LDA) | 0.775 | 0.663 | -1.861 | -1.645
$\alpha$-HgS (GGA) | -0.716 | -0.692 | -0.065 | -0.065
$\alpha$-SiO2 (LDA) | -0.071 | -0.071 | 0.125 | 0.125
$\alpha$-SiO2 (GGA) | -0.085 | -0.085 | 0.168 | 0.167

Given the large impact of SCF fields on the results, we decided to repeat our
calculations within the PBE [6] parametrization of the GGA. The corresponding
values are reported in Table 2. (Further details can be found in Ref. [36],
Sec. II.) Interestingly, for a given crystal structure the choice between LDA
and GGA seems to have a relatively small influence on the calculated
coefficients, except for the $g_{33}$ component of Se where such deviation
reaches $\sim$50%. Conversely, the structural parameters do appear to have a
significant impact on the final result. To account for this fact, we have
tested various models for the crystal structure, either using the experimental
(exp) one, or relaxed to mechanical equilibrium (either within LDA or GGA).
Our analysis shows that the fundamental gap depends on the volume of the unit
cell, and such a dependence has a strong impact on the calculated
$\boldsymbol{g}$-tensor components. For example, in the LDA equilibrium
structure of Se the electronic band gap is so small that we were unable to
converge $g_{11}$ and $g_{33}$ to meaningful values. While GGA displays the
usual overcorrection of the equilibrium volume, it yields results that are in
much closer agreement with the experiment. It is also interesting to note that
the NOA, unlike other linear-response properties (e.g. the dielectric tensor),
has a nontrivial dependence on the structure (and hence on the amplitude of
the gap). The final result originates from the mutual cancellation of several
terms, not all of which are expected to diverge in the metallic limit. This
means that some components of ${\boldsymbol{g}}$ may change rather
dramatically with structure, while others remain essentially unaltered (see
Sec. II of Ref. [36] for more details).

We now focus on the isolated molecule C4H4O2.

Table 3: Calculated independent components of the gyration tensor times the volume of the simulation cell ($\Omega$) for C4H4O2. Values are given in Hartree atomic units. | $\Omega g_{11}$ | $\Omega g_{22}$ | $\Omega g_{33}$ | $\frac{\Omega}{2}(g_{12}+g_{21})$ | $\beta$
---|---|---|---|---|---
With SCF | -69.69 | -68.12 | -33.98 | -267.32 | -2.28
Without SCF | -72.52 | -56.18 | 144.90 | -629.35 | 0.21

Table 3 shows our computed gyration tensor (multiplied by the volume of the
simulation cell $\Omega$), with and without SCF terms; as in crystals, the
latter have a huge impact on some components. We also report the optical
rotatory parameter $\beta$, which in molecular systems relates to the rotatory
power $\alpha(\omega)$ via [45, 46]

$\alpha(\omega)=\frac{N_{A}\omega^{2}}{Mc^{2}}\beta,\qquad\beta=\frac{\Omega}{4\pi}\frac{1}{2}\sum_{a}\frac{1}{3}g_{aa}.$
(20)

Here $N_{A}$ is the Avogadro number and $M$ is the molar mass of the molecule.
Our computed value of $\beta$ almost exactly matches the value of
$\beta=-2.29$ that was reported in Ref. [47]. Although such a level of
agreement gives us confidence in the correctness of our implementation, it may
be to some extent coincidental, given the differences in our respective
approximations and computational schemes.

In summary, we have presented a formulation of optical dispersion within the
framework of density-functional perturbation theory. Our methodology brings
the first-principles calculation of the gyration tensor to the same level of
accuracy and computational ease as standard linear-response properties, e.g.,
the dielectric tensor. We have also discussed some formal aspects of the
theory, e.g., the non-uniqueness of Eq. (9), which we relate to the gauge
freedom of electromagnetism. As an outlook, a natural step forward consists in
generalizing our method to finite frequencies, and to magnetic materials with
broken time-reversal symmetry; progress along these lines will be presented in
a forthcoming publication.

###### Acknowledgements.

We acknowledge support from Ministerio de Ciencia e Innovación (MICINN-Spain)
through Grant No. PID2019-108573GB-C22; from Severo Ochoa FUNFUTURE center of
excellence (CEX2019-000917-S); from Generalitat de Catalunya (Grant No. 2021
SGR 01519); and from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program (Grant Agreement No.
724529).

## References

  * Landau and Lifshitz [1984] L. D. Landau and E. M. Lifshitz,  _Electrodynamics of continuous media, Landau and Lifshitz course of theoretical physics, vol. 8_ (Pergamon Press, 1984).
  * Condon [1937] E. U. Condon, Rev. Mod. Phys. 9, 432 (1937).
  * Mason and Mason [1982] S. F. Mason and S. F. Mason,  _Molecular optical activity and the chiral discriminations_ (Cambridge University Press, 1982).
  * Melrose and McPhedran [1991] D. B. Melrose and R. C. McPhedran,  _Electromagnetic processes in dispersive media_ (Cambridge University Press, 1991).
  * Halasyamani and Poeppelmeier [1998] P. S. Halasyamani and K. R. Poeppelmeier, Chem. Mater. 10, 2753 (1998).
  * Stephens  _et al._ [2001] P. J. Stephens, F. J. Devlin, J. R. Cheeseman, and M. J. Frisch, J. Phys. Chem. A 105, 5356 (2001).
  * Autschbach  _et al._ [2002a] J. Autschbach, T. Ziegler, S. J. A. van Gisbergen, and E. J. Baerends, J. Chem. Phys. 116, 6930 (2002a).
  * Autschbach  _et al._ [2002b] J. Autschbach, S. Patchkovskii, T. Ziegler, S. J. A. van Gisbergen, and E. Jan Baerends, J. Chem. Phys. 117, 581 (2002b).
  * Polavarapu [1997] P. L. Polavarapu, Mol. Phys. 91, 551 (1997).
  * Mort and Autschbach [2005] B. C. Mort and J. Autschbach, J. Phys. Chem. A 109, 8617 (2005).
  * Jerphagnon and Chemla [1976] J. Jerphagnon and D. S. Chemla, J. Chem. Phys. 65, 1522 (1976).
  * Ades and Champness [1975] S. Ades and C. H. Champness, J. Opt. Soc. Am. 65, 217 (1975).
  * Stadnicka  _et al._ [1985] K. Stadnicka, A. Glazer, and J. Moxon, J. Appl. Cryst. 18, 237 (1985).
  * [14] S. Fukuda, T. Shiosaki, and A. Kawabata, Phys. Status Solidi B 68, K107.
  * Nomura [1960] K. C. Nomura, Phys. Rev. Lett. 5, 500 (1960).
  * Day [1971] G. W. Day, Appl. Phys. Lett. 18, 347 (1971).
  * Devarajan and Glazer [1986] V. Devarajan and A. M. Glazer, Acta Crystallogr. A 42, 560 (1986).
  * Barron [2009] L. D. Barron,  _Molecular light scattering and optical activity_ (Cambridge University Press, 2009).
  * Zhong  _et al._ [1993] H. Zhong, Z. H. Levine, D. C. Allan, and J. W. Wilkins, Phys. Rev. B 48, 1384 (1993).
  * Zhong  _et al._ [1992] H. Zhong, Z. H. Levine, D. C. Allan, and J. W. Wilkins, Phys. Rev. Lett. 69, 379 (1992).
  * Tam  _et al._ [2004] M. C. Tam, N. J. Russ, and T. D. Crawford, J. Chem. Phys. 121, 3550 (2004).
  * Norman  _et al._ [2004] P. Norman, K. Ruud, and T. Helgaker, J. Chem. Phys. 120, 5027 (2004).
  * Malashevich and Souza [2010] A. Malashevich and I. Souza, Phys. Rev. B 82, 245118 (2010).
  * Pozo and Souza [2022] Ó. Pozo and I. Souza, arXiv preprint arXiv:2211.10183 (2022).
  * Wang and Yan [2023] X. Wang and Y. Yan, Phys. Rev. B 107, 045201 (2023).
  * Tsirkin  _et al._ [2018] S. S. Tsirkin, P. A. Puente, and I. Souza, Phys. Rev. B 97, 035158 (2018).
  * Jönsson  _et al._ [1996] L. Jönsson, Z. H. Levine, and J. W. Wilkins, Phys. Rev. Lett. 76, 1372 (1996).
  * Dreyer  _et al._ [2018] C. E. Dreyer, M. Stengel, and D. Vanderbilt, Phys. Rev. B 98, 075153 (2018).
  * Prosandeev  _et al._ [2013] S. Prosandeev, A. Malashevich, Z. Gui, L. Louis, R. Walter, I. Souza, and L. Bellaiche, Phys. Rev. B 87, 195111 (2013).
  * Naumov  _et al._ [2004] I. I. Naumov, L. Bellaiche, and H. Fu, Nature 432, 737 (2004).
  * Royo and Stengel [2019] M. Royo and M. Stengel, Phys. Rev. X 9, 021050 (2019).
  * Romero  _et al._ [2020] A. H. Romero, D. C. Allan, B. Amadon, G. Antonius, T. Applencourt, L. Baguet, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, _et al._ , J. Chem. Phys. 152, 124102 (2020).
  * Gonze  _et al._ [2020] X. Gonze, B. Amadon, G. Antonius, F. Arnardi, L. Baguet, J.-M. Beuken, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, _et al._ , Comput. Phys. Commun. 248, 107042 (2020).
  * Gonze  _et al._ [2009] X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken, F. Bottin, P. Boulanger, F. Bruneval, D. Caliste, R. Caracas, M. Côté, _et al._ , Comput. Phys. Commun. 180, 2582 (2009).
  * Gonze and Lee [1997] X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997).
  * [36] See Supplemental Material at http://link for more information about theoretical derivations, computational details and results.
  * Zabalo  _et al._ [2022] A. Zabalo, C. E. Dreyer, and M. Stengel, Phys. Rev. B 105, 094305 (2022).
  * Miwa [2011] K. Miwa, Phys. Rev. B 84, 094304 (2011).
  * Essin  _et al._ [2010] A. M. Essin, A. M. Turner, J. E. Moore, and D. Vanderbilt, Phys. Rev. B 81, 205104 (2010).
  * Gonze and Zwanziger [2011] X. Gonze and J. W. Zwanziger, Phys. Rev. B 84, 064445 (2011).
  * Gonze [1995a] X. Gonze, Phys. Rev. A 52, 1096 (1995a).
  * Gonze [1995b] X. Gonze, Phys. Rev. A 52, 1086 (1995b).
  * Resta [2023] R. Resta, J. Chem. Phys. 158, 024105 (2023).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * GRI [2001] Chem. Phys. Lett. 339, 380 (2001).
  * Autschbach  _et al._ [2002c] J. Autschbach, S. Patchkovskii, T. Ziegler, S. J. A. van Gisbergen, and E. Jan Baerends, J. Chem. Phys. 117, 581 (2002c).
  * Pelloni and Lazzeretti [2014] S. Pelloni and P. Lazzeretti, J. Chem. Phys. 140, 074105 (2014).

´

##  I Supplementary notes for “Natural optical activity from density-
functional perturbation theory”

Asier Zabalo1 and Massimiliano Stengel1,2

1Institut de Ciència de Materials de Barcelona (ICMAB-CSIC), Campus UAB, 08193
Bellaterra, Spain
2ICREA-Institució Catalana de Recerca i Estudis Avançats, 08010 Barcelona,
Spain

##  II I. Computational details

Our numerical results of the main text are obtained using the DFT and DFPT
implementations of the abinit package with the Perdew-Wang [1] parametrization
of the LDA. We use norm-conserving pseudopotentials from the Pseudo Dojo [2]
website and we regenerate them without exchange-correlation nonlinear core
corrections using the ONCVPSP [3] software. Spin-orbit coupling (SOC) is
neglected in our first-principles calculations. All the materials studied in
this work, trigonal Se and $\alpha$-HgS and $\alpha$-SiO2, belong to the point
group 32, but they may crystalize in two enantiomorphic structures with space
groups $P3_{1}21$ and $P3_{2}21$. These structures with opposite handedness
have the same rotatory power in magnitude, but with opposite sign. In this
work, we shall consider the $P3_{1}21$ space group structure for the three
crystals under study. The crystal structure is either set to the experimental
one (Ref. [4] for Se), or relaxed to mechanical equilibrium until the forces
are smaller than $10^{-6}$ Ha/bohr. We use a plane-wave cutoff of 50 Ha for Se
and $\alpha$-SiO2 and 40 Ha for $\alpha$-HgS. We find that the results are
remarkably sensitive to this choice, and on the exchange-correlation
functional that is used in the structural relaxation; details are provided in
the next Section.

Regarding C4H4O2, a box with sides of $a=35$ bohr ($\Omega=a^{3}$) is used in
order to simulate an isolated molecule, with a plane-wave energy cutoff of 50
Ha. The non-relaxed coordinates are taken from Ref. [5].

Table S1: Cartesian coordinates (in bohr) of the C4H4O2 molecule, from Ref. [5]. C4H4O2 | x | y | z
---|---|---|---
O1 | 0.643397 | 1.236357 | -2.077495
O2 | -0.643397 | -1.236357 | -2.077495
C1 | -0.120062 | 1.373532 | 2.331775
C2 | 0.120062 | 2.555671 | 0.116615
C3 | 0.120062 | -1.373532 | 2.331775
C4 | -0.120062 | -2.555671 | 0.116615
H1 | 0.117087 | 4.574339 | -0.221634
H2 | -0.496036 | 2.427108 | 4.040414
H3 | -0.117087 | -4.574339 | -0.221634
H4 | 0.496036 | -2.427108 | 4.040414

##  III II. Effect of the structural parameters and exchange and correlation
functionals (LDA vs GGA)

The aim of this section is to examine how the exchange and correlation
functional affects the structural characteristics of the solids being
analyzed. We shall see how going from LDA to GGA changes the structural
parameters and the electronic band structure of Se, $\alpha$-HgS and
$\alpha$-SiO2 and see how this affects the dielectric tensor, and hence, the
natural optical activity. Table S2 shows the structural parameters for
different cases. As expected, GGA gives a larger unit cell when the forces are
relaxed and, in general, a larger electronic band gap (the only exception here
is $\alpha$-SiO2).

Table S2: Structural parameters and electronic band gap of Se, $\alpha$-HgS and $\alpha$-SiO2, with LDA and GGA. In the Structure column, “exp” refers to the experimental structure, while Se (LDA) means that the structure was relaxed with LDA, for example. The Functional column specifies the exchange and correlation functional used for computing the electronic band structure. Structure | Functional |  $a$ (bohr) |  $c$ (bohr) |  $\Omega$ (bohr3) | band gap (eV)
---|---|---|---|---|---
Se (exp) | LDA | 8.201 | 9.354 | 544.894 | 0.850
Se (exp) | GGA | 8.201 | 9.354 | 544.894 | 0.948
Se (LDA) | LDA | 7.431 | 9.695 | 463.605 | 0.225
Se (GGA) | GGA | 8.463 | 9.564 | 593.187 | 0.961
$\alpha$-HgS (LDA) | LDA | 7.629 | 17.530 | 883.498 | 0.995
$\alpha$-HgS (GGA) | GGA | 8.434 | 18.430 | 1135.414 | 1.829
$\alpha$-SiO2 (LDA) | LDA | 9.174 | 10.110 | 736.974 | 5.987
$\alpha$-SiO2 (GGA) | GGA | 9.489 | 10.411 | 811.792 | 5.768

Fig. S1, S2 and S3 show the electronic band structure of Se, $\alpha$-HgS and
$\alpha$-SiO2, respectively. For Se, we study four cases: the non-relaxed
structure and the relaxed structure, both with LDA and GGA. When the structure
is relaxed with LDA, the band gap drastically decreases and, in fact, it
becomes so small (around 0.2 eV) that even standard DFPT calculations, e.g.
the quadrupoles, fail to converge. In general, we observe that the larger the
band gap, the faster the convergence is, which might explain the slightly
faster convergence of the natural optical activity with GGA, in contrast to
LDA, that is shown in Sec. III. The most intriguing discovery arises when
examining the impact of structural parameters on the dielectric tensor, as
shown in Table S3. Surprisingly, the choice of the exchange and correlation
functional has minimal effect on the independent components of the dielectric
tensor. However, the structural parameters play a significant role here,
particularly evident in the case of $\alpha$-HgS. Although the effect is less
dramatic for Se and $\alpha$-SiO2, it becomes more pronounced when calculating
the natural optical activity tensor, which is just the first-order derivative
of the dielectric tensor with respect to $\mathbf{q}$ (see Table II of main
text).

Figure S1: Band structure of trigonal Se, computed with LDA and GGA.

Figure S2: Band structure of $\alpha$-HgS, computed with LDA and GGA.

Figure S3: Band structure of $\alpha$-SiO2, computed with LDA and GGA. Table S3: Comparison between LDA and GGA for the independent components of the dielectric permittivity tensor for Se, $\alpha$-HgS and $\alpha$-SiO2, for different structures, in Hartree atomic units. In the Structure column, “exp” refers to the experimental structure, while Se (LDA) means that the structure was relaxed with LDA, for example. Structure | $\epsilon_{11}$ | $\epsilon_{33}$
---|---|---
| LDA | GGA | LDA | GGA
Se (exp) | 9.05817 | 8.44465 | 14.66883 | 13.89905
Se (GGA) | 9.19622 | 8.53039 | 14.77391 | 13.92631
$\alpha$-HgS (LDA) | 10.07620 | 9.40749 | 12.10559 | 11.23652
$\alpha$-HgS (GGA) | 5.89176 | 5.56827 | 8.21707 | 7.72728
$\alpha$-SiO2 (LDA) | 2.48932 | 2.45577 | 2.52075 | 2.48670
$\alpha$-SiO2 (GGA) | 2.42803 | 2.39744 | 2.45630 | 2.42505

To gain further insight into the effect of the structural parameters, in Table
S4 we show the individual contributions to the NOA tensor for different
structures. In general, the effect of the structure is clearly noticeable in
the individual terms, which can be easily appreciated for $\alpha$-HgS.
However, even if the individual contributions differ to a large degree (see
the $\mathcal{Y}$ and $\mathcal{W}$ terms in Se), an accidental cancellation
of several terms can lead to a similar final result for the NOA tensor
components.

Table S4: Contribution to the $g_{33}$ tensor component (in bohr) from different terms of Eq. (10) and Eq. (11) of the main text. We are showing the results obtained with the LDA exchange and correlation functional for the computation of the NOA, in all cases. Structure | $E_{\text{elst}}$ | $\mathcal{X}$ | $\mathcal{Y}$ | $\mathcal{W}$ | $g_{33}$
---|---|---|---|---|---
Se (exp) | -0.164 | -0.836 | -1.884 | 0.973 | -1.910
Se (GGA) | -0.169 | -0.729 | -2.482 | 1.578 | -1.802
$\alpha$-HgS (LDA) | -0.095 | -1.106 | -0.801 | 0.149 | -1.861
$\alpha$-HgS (GGA) | -0.204 | -0.054 | -0.479 | 0.672 | -0.065
$\alpha$-SiO2 (LDA) | 0.010 | 0.020 | 0.099 | -0.005 | 0.125
$\alpha$-SiO2 (GGA) | 0.016 | 0.029 | 0.118 | 0.006 | 0.168

##  IV III. Convergence study of the natural optical activity tensor

As a further step for the validation of our method, we plot the independent
components of the gyration tensor, $g_{11}$ and $g_{33}$, for trigonal Se,
$\alpha$-HgS and $\alpha$-SiO2 . Fig. S4, S5 and S6 show the obtained
numerical results as a function of the plane-wave cutoff energy and the
$\mathbf{k}$ points mesh resolution. For comparison purposes, we also compute
the natural optical activity with the PBE [6] parametrization of the GGA. (For
trigonal Se, the experimental structure is used both for LDA and GGA; whereas
for $\alpha$-HgS and $\alpha$-SiO2, the structure is relaxed until the forces
are smaller than $10^{-6}$ Ha/bohr, either with LDA or GGA.)

###  IV.1 Trigonal Se

Figure S4: Convergence of the independent components of the gyration tensor,
$g_{11}$ and $g_{33}$, of non-relaxed trigonal Se, with respect to the plane-
wave cutoff and the density of the $\mathbf{k}$-point mesh. The top panel
shows the obtained results with LDA, whereas the bottom panel shows the
results obtained with GGA.

###  IV.2 $\alpha$-HgS

Figure S5: Convergence of the independent components of the gyration tensor,
$g_{11}$ and $g_{33}$, of $\alpha$-HgS, with respect to the plane-wave cutoff
and the density of the $\mathbf{k}$-point mesh. The top panel shows the
obtained results with LDA, whereas the bottom panel shows the results obtained
with GGA.

###  IV.3 $\alpha$-SiO2

Figure S6: Convergence of the independent components of the gyration tensor,
$g_{11}$ and $g_{33}$, of $\alpha$-SiO2, with respect to the plane-wave cutoff
and the density of the $\mathbf{k}$-point mesh. The top panel shows the
obtained results with LDA, whereas the bottom panel shows the results obtained
with GGA.

##  V IV. Response to a long-wavelength magnetic vector potential field

The wave-function response to a long-wavelength magnetic vector potential
field is given by the following Sternheimer equation [7, 8]

$(\hat{H}_{\mathbf{k}}^{(0)}+a\hat{P}_{\mathbf{k}}-\epsilon_{\mathbf{k}}^{(0)})\ket*{u_{m\mathbf{k},\gamma}^{A_{\beta}}}=-\hat{Q}_{\mathbf{k}}\hat{O}_{\mathbf{k}}^{\alpha\gamma}\ket*{u_{m\mathbf{k}}^{(0)}},$
(21)

where $\hat{H}^{(0)}_{\mathbf{k}}$ is the ground state Hamiltonian,
$\epsilon_{m\mathbf{k}}^{(0)}$ is an unperturbed energy eigenvalue and
$\hat{Q}_{\mathbf{k}}=1-\hat{P}_{\mathbf{k}}$ is the conduction-band projector
with
$\hat{P}_{\mathbf{k}}=\sum_{n}\ket*{u_{n\mathbf{k}}^{(0)}}\bra*{u_{n\mathbf{k}}^{(0)}}$.
In Eq. (21), $a$ is a parameter with the dimension of energy that ensures
stability [7]. The perturbing operator $\hat{O}_{\mathbf{k}}^{\alpha\gamma}$
is given by

$\hat{O}_{\mathbf{k}}^{\alpha\gamma}=\partial_{\gamma}\hat{H}_{\mathbf{k}}^{(0)}\partial_{\alpha}\hat{P}_{\mathbf{k}}-\partial_{\gamma}\hat{P}_{\mathbf{k}}\partial_{\alpha}\hat{H}_{\mathbf{k}}^{(0)}+\frac{1}{2}\partial^{2}_{\alpha\gamma}\hat{H}_{\mathbf{k}}^{(0)},$
(22)

where
$\partial^{2}_{\gamma\alpha}\equiv\partial^{2}/\partial_{k_{\alpha}}\partial_{k_{\gamma}}$
is a second derivative in $\mathbf{k}$ space. (Notice that Eq. (22) gives the
correct response to a long-wavelength magnetic vector potential field. There
is a typo in Ref. [8], which makes Eq. (22) differ from the expression given
in Ref. [8] by a half factor.)

The perturbing operator $\hat{O}_{\mathbf{k}}^{\beta\gamma}$ can be expressed
as a sum of a symmetric $\hat{O}_{\mathbf{k}}^{S,\beta\gamma}$ and an
antisymmetric $\hat{O}_{\mathbf{k}}^{A,\beta\gamma}$ part with respect to
$\beta\leftrightarrow\gamma$ exchange, in a way that the symmetric part
recovers the usual $d^{2}/dkdk$ wave functions,
$\ket*{\partial^{2}_{\beta\gamma}u_{m\mathbf{k}}^{(0)}}$, and the
antisymmetric part recovers the wave function response to a uniform orbital
$\mathbf{B}$ field, $\ket*{u_{m\mathbf{k}}^{B_{\delta}}}$, such that

$\hat{O}_{\mathbf{k}}^{\beta\gamma}=\frac{1}{2}\big{(}\hat{O}_{\mathbf{k}}^{S,\beta\gamma}+\hat{O}_{\mathbf{k}}^{A,\beta\gamma}\big{)},$
(23)

where $\hat{O}_{\mathbf{k}}^{S,\beta\gamma}$ is the perturbing operator that
enters the Sternheimer equation for the
$\ket*{\partial^{2}_{\beta\gamma}u_{m\mathbf{k}}^{(0)}}$ wave functions,

$\hat{O}_{\mathbf{k}}^{S,\beta\gamma}=\partial_{\gamma}\hat{H}^{(0)}_{\mathbf{k}}\partial_{\beta}\hat{P}_{\mathbf{k}}+\partial_{\beta}\hat{H}_{\mathbf{k}}^{(0)}\partial_{\gamma}\hat{P}_{\mathbf{k}}-\partial_{\gamma}\hat{P}_{\mathbf{k}}\partial_{\beta}\hat{H}^{(0)}_{\mathbf{k}}-\partial_{\beta}\hat{P}_{\mathbf{k}}\partial_{\gamma}\hat{H}^{(0)}_{\mathbf{k}}+\partial^{2}_{\beta\gamma}\hat{H}^{(0)}_{\mathbf{k}},$
(24)

and $\hat{O}_{\mathbf{k}}^{A,\beta\gamma}$ is intimately related to the
perturbing operator that enters the Sternheimer equation for the wave function
response to a uniform orbital magnetic field,
$\hat{O}_{\mathbf{k}}^{B_{\delta}}$, as defined by Essin et al. [9],

$\hat{O}_{\mathbf{k}}^{B_{\delta}}=\frac{i}{2}\epsilon_{\delta\beta\gamma}\hat{O}_{\mathbf{k}}^{A,\beta\gamma},$
(25)

where

$\hat{O}_{\mathbf{k}}^{A,\beta\gamma}=\partial_{\gamma}\hat{H}_{\mathbf{k}}^{(0)}\partial_{\beta}\hat{P}_{\mathbf{k}}-\partial_{\beta}\hat{H}_{\mathbf{k}}^{(0)}\partial_{\gamma}\hat{P}_{\mathbf{k}}-\partial_{\gamma}\hat{P}_{\mathbf{k}}\partial_{\beta}\hat{H}_{\mathbf{k}}^{(0)}+\partial_{\beta}\hat{P}_{\mathbf{k}}\partial_{\gamma}\hat{H}_{\mathbf{k}}^{(0)}.$
(26)

Using the above relations one can readily verify that Eq. (15), (16) and (17)
of the main text are nicely fulfilled.

##  VI V. Natural optical activity without SCF terms

We shall use tilded symbols to indicate that SCF terms are excluded. Without
self-consistency, the natural optical activity tensor reads as

$\tilde{\eta}_{\alpha\beta\gamma}=-\frac{4\pi}{\Omega}\text{Im}\tilde{E}_{\gamma}^{\mathcal{E}_{\alpha}^{*}\mathcal{E}_{\beta}},$
(27)

where

$\tilde{E}_{\gamma}^{\mathcal{E}_{\alpha}^{*}\mathcal{E}_{\beta}}=2s\int_{\text{BZ}}[d^{3}k]\tilde{E}_{\mathbf{k},\gamma}^{\mathcal{E}_{\alpha}^{*}\mathcal{E}_{\beta}}.$
(28)

The wave function term is given by

$\begin{split}\tilde{E}_{\mathbf{k},\gamma}^{\mathcal{E}_{\alpha}^{*}\mathcal{E}_{\beta}}&=\tilde{\mathcal{X}}_{\mathbf{k}}^{\mathcal{E}_{\alpha}k_{\gamma}\mathcal{E}_{\beta}}+\tilde{\mathcal{W}}_{\mathbf{k}}^{\alpha,\beta\gamma}+\big{[}\tilde{\mathcal{W}}_{\mathbf{k}}^{\beta,\alpha\gamma}\big{]}^{*}\\\
&=\sum_{m}\Big{(}\bra*{\tilde{u}_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}\hat{H}_{\mathbf{k}}^{k_{\gamma}}\ket*{\tilde{u}_{m\mathbf{k}}^{\mathcal{E}_{\beta}}}+\bra*{\tilde{u}_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}\ket*{iu_{m\mathbf{k},\gamma}^{A_{\beta}}}+\bra*{iu_{m\mathbf{k},\gamma}^{A_{\alpha}}}\ket*{\tilde{u}_{m\mathbf{k}}^{\mathcal{E}_{\beta}}}\Big{)}.\end{split}$
(29)

In general, when self-consistency is taken into account, the first-order wave
functions $\ket*{u_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}$ are obtained by
solving the following Sternheimer equation,

$(\hat{H}^{(0)}_{\mathbf{k}}+a\hat{P}_{\mathbf{k}}-\epsilon_{m\mathbf{k}}^{(0)})\ket*{u_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}=-\hat{Q}_{\mathbf{k}}\hat{\mathcal{H}}_{\mathbf{k}}^{\mathcal{E}_{\alpha}}\ket*{u_{m\mathbf{k}}^{(0)}},$
(30)

where
$\hat{\mathcal{H}}_{\mathbf{k}}^{\mathcal{E}_{\alpha}}=i\partial_{\alpha}\hat{P}_{\mathbf{k}}+\hat{V}^{\mathcal{E}_{\alpha}}$,
which includes self-consistency via the usual SCF potential
$\hat{V}^{\mathcal{E}_{\alpha}}$. The tilded first-order wave functions
$\ket*{\tilde{u}_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}$ are obtained from Eq.
(30), once the $\hat{V}^{\mathcal{E}_{\alpha}}$ term is excluded,

$(\hat{H}_{\mathbf{k}}^{(0)}+a\hat{P}_{\mathbf{k}}-\epsilon_{m\mathbf{k}}^{(0)})\ket*{\tilde{u}_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}=-i\ket*{u_{m\mathbf{k}}^{k_{\alpha}}}.$
(31)

##  VII VI. Recovering the formulas of “Band theory of spatial dispersion in
magnetoelectrics” for finite samples

This section is devoted to establish a formal link between our long-wave DFPT-
based expression for the natural optical activity and the formulas obtained in
Ref. [10] for finite systems. First of all, we note that our tilded (see Sec.
V) natural optical activity tensor, $\tilde{\eta}_{\alpha\beta\gamma}$, is
directly related to the $\mathcal{T}$-even part (where $\mathcal{T}$
represents time-reversal symmetry) of the first $\mathbf{q}$ gradient of the
effective conductivity tensor computed in Ref. [10],
$\sigma_{\alpha\beta\gamma}$,

$\tilde{\eta}_{\alpha\beta\gamma}=\frac{4\pi}{\omega}\sigma^{\text{A}}_{\alpha\beta\gamma}(\omega)\Big{|}_{\omega\longrightarrow
0}.$ (32)

Our starting point is our Eq. (27). In full generality, one can express the
first-order wave functions as follows,

$\ket*{\tilde{u}_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}=\sum_{m,l}^{o,e}\ket*{u_{l\mathbf{k}}^{(0)}}\frac{\bra*{u_{l\mathbf{k}}^{(0)}}\hat{H}_{\mathbf{k}}^{\mathcal{E}_{\alpha}}\ket*{u_{m\mathbf{k}}^{(0)}}}{\epsilon_{m\mathbf{k}}^{(0)}-\epsilon_{l\mathbf{k}}^{(0)}},$
(33)

where the $m$ and $l$ indices run over the occupied ($o$) and empty ($e$)
bands, respectively. For finite samples

$\hat{H}_{\mathbf{k}}^{\mathcal{E}_{\alpha}}\longrightarrow
r_{\alpha},\quad\ket*{u_{n\mathbf{k}}^{(0)}}\longrightarrow\ket*{n},\quad\epsilon_{n\mathbf{k}}^{(0)}\longrightarrow\epsilon_{n}.$
(34)

(Also, and to simplify the notation, we shall use
$\hat{H}^{(0)}\rightarrow\hat{H}$.) Along the same lines, it is useful to
recall that

$\partial_{\gamma}\hat{H}^{(0)}_{\mathbf{k}}\longrightarrow
i[\hat{H},r_{\gamma}],\quad\partial_{\gamma}\hat{P}_{\mathbf{k}}\longrightarrow
i[\hat{P},r_{\gamma}].$ (35)

###  VII.1 Computing
$\tilde{\mathcal{X}}_{\mathbf{k}}^{\mathcal{E}_{\alpha}k_{\gamma}\mathcal{E}_{\beta}}$

$\begin{split}\tilde{\mathcal{X}}_{\mathbf{k}}^{\mathcal{E}_{\alpha}k_{\gamma}\mathcal{E}_{\beta}}&=\sum_{m}^{o}\bra*{\tilde{u}_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}\hat{H}_{\mathbf{k}}^{k_{\gamma}}\ket*{\tilde{u}_{m\mathbf{k}}^{\mathcal{E}_{\beta}}}\\\
&=i\sum_{m,ls}^{o,e}\Big{[}\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}[\hat{H},r_{\gamma}]\ket*{s}\bra*{s}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})(\epsilon_{s}-\epsilon_{m})}\Big{]}\\\
&=i\sum_{m,ls}^{o,e}\Big{[}\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}\hat{H}r_{\gamma}\ket*{s}\bra*{s}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})(\epsilon_{s}-\epsilon_{m})}-\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}r_{\gamma}\hat{H}\ket*{s}\bra*{s}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})(\epsilon_{s}-\epsilon_{m})}\Big{]}\\\
&=i\sum_{m,ls}^{o,e}\Big{[}(\epsilon_{l}-\epsilon_{s})\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}r_{\gamma}\ket*{s}\bra*{s}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})(\epsilon_{s}-\epsilon_{m})}\Big{]}\\\
&=i\sum_{m,ls}^{o,e}\Big{[}(\epsilon_{l}-\epsilon_{m})\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}r_{\gamma}\ket*{s}\bra*{s}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})(\epsilon_{s}-\epsilon_{m})}+(\epsilon_{m}-\epsilon_{s})\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}r_{\gamma}\ket*{s}\bra*{s}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})(\epsilon_{s}-\epsilon_{m})}\Big{]}\\\
&=i\sum_{m,ls}^{o,e}\Big{[}\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}r_{\gamma}\ket*{s}\bra*{s}r_{\beta}\ket*{m}}{(\epsilon_{s}-\epsilon_{m})}-\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}r_{\gamma}\ket*{s}\bra*{s}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})}\Big{]}\\\
&=i\sum_{m,l}^{o,e}\Big{[}\frac{\bra*{m}r_{\alpha}\hat{Q}r_{\gamma}\ket*{l}\bra*{l}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})}-\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}r_{\gamma}\hat{Q}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})}\Big{]},\end{split}$
(36)

where $m$ runs over the occupied states $(o)$ only and $s$ and $l$ run over
the empty states $(e)$. We have used
$\epsilon_{l}-\epsilon_{s}=(\epsilon_{l}-\epsilon_{m})+(\epsilon_{m}-\epsilon_{s})$
in the fifth line and $\hat{Q}=\sum_{l}\ket*{l}\bra*{l}$.

###  VII.2 Computing $\tilde{\mathcal{W}}_{\mathbf{k}}^{\alpha,\beta\gamma}$

$\begin{split}\tilde{\mathcal{W}}_{\mathbf{k}}^{\alpha,\beta\gamma}&=i\sum_{m}^{o}\bra*{\tilde{u}_{m\mathbf{k}}^{\mathcal{E}_{\alpha}}}\ket*{u^{A_{\beta}}_{m\mathbf{k},\gamma}}\\\
&=i\sum_{m,l}^{o,e}\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}\ket*{u_{m\mathbf{k},\gamma}^{A_{\beta}}}}{\epsilon_{m}-\epsilon_{l}}\\\
&=i\sum_{m,l}^{o,e}\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}\hat{O}^{\beta\gamma}\ket*{m}}{(\epsilon_{m}-\epsilon_{l})^{2}},\end{split}$
(37)

where we have used Eq. (21) in order to write

$\bra*{l}\ket*{u_{m\mathbf{k},\gamma}^{A_{\beta}}}\longrightarrow\frac{\bra*{l}\hat{O}^{\beta\gamma}\ket*{m}}{\epsilon_{m}-\epsilon_{l}}.$
(38)

We shall now focus on computing

$\bra*{l}\hat{O}^{\beta\gamma}\ket*{m}=\frac{1}{2}\bra*{l}\big{(}\hat{O}^{\text{S},\beta\gamma}+\hat{O}^{\text{A},\beta\gamma}\big{)}\ket*{m}.$
(39)

####  VII.2.1 Symmetric part

$\begin{split}\bra*{l}\hat{O}^{\text{S},\beta\gamma}\ket*{m}=&\underbrace{-\bra*{l}[\hat{H},r_{\gamma}][\hat{P},r_{\beta}]\ket*{m}-\bra*{l}[\hat{H},r_{\beta}][\hat{P},r_{\gamma}]\ket*{m}-\bra*{l}[[\hat{H},r_{\beta}],r_{\gamma}]\ket*{m}}_{S_{1}}\\\
&\underbrace{+\bra*{l}[\hat{P},r_{\beta}][\hat{H},r_{\gamma}]\ket*{m}+\bra*{l}[\hat{P},r_{\gamma}][\hat{H},r_{\beta}]\ket*{m}}_{S_{2}}.\end{split}$
(40)

Let us first focus on the $S_{1}$ term,

$\begin{split}S_{1}=&\underbrace{-\bra*{l}[\hat{H},r_{\gamma}]r_{\beta}\ket*{m}+\bra*{l}[\hat{H},r_{\beta}]r_{\gamma}\ket*{m}-\bra*{l}[[\hat{H},r_{\beta}],r_{\gamma}]\ket*{m}}_{S_{1}^{(a)}}\underbrace{-\bra*{l}[\hat{H},r_{\gamma}]\hat{P}r_{\beta}\ket*{m}-\bra*{l}[\hat{H},r_{\beta}]\hat{P}r_{\gamma}\ket*{m}}_{S_{2}^{(b)}}.\end{split}$
(41)

The $S_{1}^{(a)}$ term can be simplified as follows

$\begin{split}S_{1}^{(a)}=&\bra*{l}\hat{H}r_{\beta}r_{\gamma}\ket*{m}-\bra*{l}r_{\beta}\hat{H}r_{\gamma}\ket*{m}+\bra*{l}\hat{H}r_{\beta}r_{\gamma}\ket*{m}-\bra*{l}r_{\gamma}\hat{H}r_{\beta}\ket*{m}-\bra*{l}[\hat{H},r_{\beta}]r_{\gamma}\ket*{m}+\bra*{l}r_{\gamma}[\hat{H},r_{\beta}]\ket*{m}\\\
=&\bra*{l}\hat{H}r_{\gamma}r_{\beta}\ket*{m}-\bra*{l}r_{\gamma}r_{\beta}\ket*{m}\\\
=&\bra*{l}[\hat{H},r_{\beta}r_{\gamma}]\ket*{m}.\end{split}$ (42)

In the same way, we can simplify the $S_{1}^{(b)}$ term, where only to terms
survive, as we have a conduction-band state on the left,

$\begin{split}S_{1}^{(b)}&=\bra*{l}[\hat{P},r_{\beta}][\hat{H}.r_{\gamma}]\ket*{m}+\bra*{l}[\hat{P},r_{\gamma}][\hat{H},r_{\beta}]\ket*{m}\\\
&=-\bra*{l}r_{\beta}\hat{P}[\hat{H},r_{\gamma}]\ket*{m}-\bra*{l}r_{\gamma}\hat{P}[\hat{H},r_{\beta}]\ket*{m}.\end{split}$
(43)

Collecting all terms, $S_{1}^{(a)}+S_{1}^{(b)}+S_{2}$, we obtain

$\begin{split}\bra*{l}\hat{O}^{\text{S},\beta\gamma}\ket*{m}=&\bra*{l}[\hat{H},r_{\beta}r_{\gamma}]\ket*{m}-\bra*{l}[\hat{H},r_{\gamma}]\hat{P}r_{\beta}\ket*{m}-\bra*{l}[\hat{H},r_{\beta}]\hat{P}r_{\gamma}\ket*{m}-\bra*{l}r_{\beta}\hat{P}[\hat{H},r_{\gamma}]\ket*{m}-\bra*{l}r_{\gamma}\hat{P}[\hat{H},r_{\beta}]\ket*{m}\\\
=&\bra*{l}[\hat{H},r_{\beta}r_{\gamma}]\ket*{m}-\bra*{l}\hat{H}r_{\gamma}\hat{P}r_{\beta}\ket*{m}+\cancel{\bra*{l}r_{\gamma}\hat{H}\hat{P}r_{\beta}\ket*{m}}-\bra*{l}\hat{H}r_{\beta}\hat{P}r_{\gamma}\ket*{m}+\cancel{\bra*{l}r_{\beta}\hat{H}\hat{P}r_{\gamma}\ket*{m}}\\\
&-\cancel{\bra*{l}r_{\beta}\hat{P}\hat{H}r_{\gamma}\ket*{m}}+\bra*{l}r_{\beta}\hat{P}r_{\gamma}\hat{H}\ket*{m}-\cancel{r_{\gamma}\hat{P}\hat{H}r_{\beta}\ket*{m}}+\bra*{l}r_{\gamma}\hat{P}r_{\beta}\hat{H}\ket*{m}\\\
=&\bra*{l}[\hat{H},r_{\beta}r_{\gamma}]\ket*{m}-\bra*{l}\hat{H}r_{\gamma}\hat{P}r_{\beta}\ket*{m}-\bra*{l}\hat{H}r_{\beta}\hat{P}r_{\gamma}\ket*{m}+\bra*{l}r_{\beta}\hat{P}r_{\gamma}\hat{H}\ket*{m}+\bra*{l}r_{\gamma}\hat{P}r_{\beta}\hat{H}\ket*{m}.\end{split}$
(44)

Thus the symmetric part yields a purely geometric contribution to the wave
function response,

$\bra*{l}\hat{O}^{\text{S},\beta\gamma}\ket*{m}=(\epsilon_{l}-\epsilon_{m})\bra*{l}\big{(}r_{\beta}r_{\gamma}-r_{\beta}\hat{P}r_{\gamma}-r_{\gamma}\hat{P}r_{\beta}\big{)}\ket*{m}.$
(45)

####  VII.2.2 Antisymmetric part

$\begin{split}\bra*{l}\hat{O}^{\text{A},\beta\gamma}\ket*{m}&=\bra*{l}\big{(}\partial_{\gamma}\hat{H}\partial_{\beta}\hat{P}-\partial_{\beta}\hat{H}\partial_{\gamma}\hat{P}-\partial_{\gamma}\hat{P}\partial_{\beta}\hat{H}+\partial_{\beta}\hat{P}\partial_{\gamma}\hat{H}\big{)}\ket*{m}\\\
&=\underbrace{-\bra*{l}[\hat{H},r_{\gamma}][\hat{P},r_{\beta}]\ket*{m}+\bra*{l}[\hat{H},r_{\beta}][\hat{P},r_{\gamma}]\ket*{m}}_{A_{1}}+\underbrace{\bra*{l}[\hat{P},r_{\gamma}][\hat{H},r_{\beta}]\ket*{m}-\bra*{l}[\hat{P},r_{\beta}][\hat{H},r_{\gamma}]\ket*{m}}_{A_{2}}.\end{split}$
(46)

Regarding the $A_{1}$ term, we have

$\begin{split}A_{1}=&-\bra*{l}\big{(}\hat{H}r_{\gamma}-r_{\gamma}\hat{H}\big{)}\big{(}\hat{P}r_{\beta}-r_{\beta}\hat{P}\big{)}\ket*{m}+\bra*{l}\big{(}\hat{H}r_{\beta}-r_{\beta}\hat{H}\big{)}\big{(}\hat{P}r_{\gamma}-r_{\gamma}\hat{P}\big{)}\ket*{m}\\\
=&-\bra*{l}\hat{H}r_{\gamma}\hat{P}r_{\beta}\ket*{m}+\cancel{\bra*{l}\hat{H}r_{\gamma}r_{\beta}\ket*{m}}+\bra*{l}r_{\gamma}\hat{H}\hat{P}r_{\beta}\ket*{m}-\bra*{l}r_{\gamma}\hat{H}r_{\beta}\ket*{m}+\bra*{l}\hat{H}r_{\beta}\hat{P}r_{\gamma}\ket*{m}\\\
&-\cancel{\bra*{l}\hat{H}r_{\beta}r_{\gamma}\ket*{m}}-\bra*{l}r_{\beta}\hat{H}\hat{P}r_{\gamma}\ket*{m}+\bra*{l}r_{\beta}\hat{H}r_{\gamma}\ket*{m}\\\
=&-\epsilon_{l}\bra*{l}r_{\gamma}\hat{P}r_{\beta}\ket*{m}+\bra*{l}r_{\gamma}\hat{H}\hat{P}r_{\beta}\ket*{m}-\bra*{l}r_{\gamma}\hat{H}r_{\beta}\ket*{m}+\epsilon_{m}\bra*{l}r_{\beta}\hat{P}r_{\gamma}\ket*{m}-\bra*{l}r_{\beta}\hat{H}\hat{P}r_{\gamma}\ket*{m}\\\
&+\bra*{l}r_{\beta}\hat{H}r_{\gamma}\ket*{m}.\end{split}$ (47)

For $A_{2}$ some terms identically vanish, as we have a conduction-band state
on the left,

$\begin{split}A_{2}=&\bra*{l}[\hat{P},r_{\gamma}][\hat{H},r_{\beta}]\ket*{m}-\bra*{l}[\hat{P},r_{\beta}][\hat{H},r_{\gamma}]\ket*{m}\\\
=&\bra*{l}\big{(}\hat{P}r_{\gamma}-r_{\gamma}\hat{P}\big{)}\big{(}\hat{H}r_{\beta}-r_{\beta}\hat{H}\big{)}\ket*{m}-\bra*{l}\big{(}\hat{P}r_{\beta}-r_{\beta}\hat{P}\big{)}\big{(}\hat{H}r_{\gamma}-r_{\gamma}\hat{H}\big{)}\ket*{m}\\\
=&-\bra*{l}r_{\gamma}\hat{P}\hat{H}r_{\beta}\ket*{m}+\bra*{l}r_{\gamma}\hat{P}r_{\beta}\hat{H}\ket*{m}+\bra*{l}r_{\beta}\hat{P}\hat{H}r_{\gamma}\ket*{m}-\bra*{l}r_{\beta}\hat{P}r_{\gamma}\hat{H}\ket*{m}\\\
=&-\bra*{l}r_{\gamma}\hat{P}\hat{H}r_{\beta}\ket*{m}+\epsilon_{m}\bra*{l}r_{\gamma}\hat{P}r_{\beta}\ket*{m}+\bra*{l}r_{\beta}\hat{P}\hat{H}r_{\gamma}\ket*{m}-\epsilon_{m}\bra*{l}r_{\beta}\hat{P}r_{\gamma}\ket*{m}.\end{split}$
(48)

Everything together,

$\bra*{l}\hat{O}^{\text{A},\beta\gamma}\ket*{m}=\big{(}\epsilon_{l}-\epsilon_{m}\big{)}\bra*{l}\big{(}r_{\beta}\hat{P}r_{\gamma}-r_{\gamma}\hat{P}r_{\beta}\big{)}\ket*{m}+\bra*{l}\big{(}r_{\beta}\hat{H}r_{\gamma}-r_{\gamma}\hat{H}r_{\beta}\big{)}\ket*{m}.$
(49)

For convenience, we shall write the second term as

$\begin{split}\bra*{l}\big{(}r_{\beta}\hat{H}r_{\gamma}-r_{\gamma}\hat{H}r_{\beta}\big{)}\ket*{m}=&\bra*{l}\big{(}r_{\beta}[\hat{H},r_{\gamma}]-r_{\gamma}[\hat{H},r_{\beta}]\big{)}\ket*{m}\\\
=&-i\bra*{l}\big{(}r_{\beta}v_{\gamma}-r_{\gamma}v_{\beta}\big{)}\ket*{m},\end{split}$
(50)

where we have used $v_{\alpha}=i[\hat{H},r_{\alpha}]$. Therefore, the
antisymmetric contribution reads as

$\bra*{l}\hat{O}^{\text{A},\beta\gamma}\ket*{m}=-i\bra*{l}\big{(}r_{\beta}v_{\gamma}-r_{\gamma}v_{\beta}\big{)}\ket*{m}+\big{(}\epsilon_{l}-\epsilon_{m}\big{)}\bra*{l}\big{(}r_{\beta}\hat{P}r_{\gamma}-r_{\gamma}\hat{P}r_{\beta}\big{)}\ket*{m}.$
(51)

####  VII.2.3 Symmetric + Antisymmetric contribution

We can now compute the $\tilde{\mathcal{W}}_{\mathbf{k}}^{\alpha,\beta\gamma}$
term given by Eq. (37),

$\begin{split}\tilde{\mathcal{W}}_{\mathbf{k}}^{\alpha,\beta\gamma}=&\frac{i}{2}\sum_{m,l}^{o,e}\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}\big{(}\hat{O}^{\text{S},\beta\gamma}+\hat{O}^{\text{A},\beta\gamma}\big{)}\ket*{m}}{(\epsilon_{m}-\epsilon_{l})^{2}}\\\
=&\frac{i}{2}\sum_{m,l}^{o,e}\frac{\bra*{m}r_{\alpha}\ket*{l}}{(\epsilon_{m}-\epsilon_{l})^{2}}\Big{[}\big{(}\epsilon_{l}-\epsilon_{m}\big{)}\bra*{l}\big{(}r_{\beta}r_{\gamma}-\cancel{r_{\beta}\hat{P}r_{\gamma}}-r_{\gamma}\hat{P}r_{\beta}\big{)}\ket*{m}\\\
&+\big{(}\epsilon_{l}-\epsilon_{m}\big{)}\bra*{l}\big{(}\cancel{r_{\beta}\hat{P}r_{\gamma}}-r_{\gamma}\hat{P}r_{\beta}\big{)}\ket*{m}-i\bra*{l}\big{(}r_{\beta}v_{\gamma}-r_{\gamma}v_{\beta}\big{)}\ket*{m}\Big{]}\\\
=&\frac{i}{2}\sum_{m,l}^{o,e}\frac{\bra*{m}r_{\alpha}\ket*{l}}{(\epsilon_{l}-\epsilon_{m})}\bra*{l}\big{(}r_{\beta}r_{\gamma}-2r_{\gamma}\hat{P}r_{\beta}\big{)}\ket*{m}+\frac{1}{2}\sum_{m,l}^{o,e}\frac{\bra*{m}r_{\alpha}\ket*{l}}{(\epsilon_{l}-\epsilon_{m})^{2}}\bra*{l}\big{(}r_{\beta}v_{\gamma}-r_{\gamma}v_{\beta}\big{)}\ket*{m}.\end{split}$
(52)

###  VII.3 All terms together: computing $\tilde{\eta}_{\alpha\beta\gamma}$

$\begin{split}\tilde{\eta}_{\alpha\beta\gamma}=&-\frac{4\pi}{\Omega}2s\text{
Im}\Big{(}\tilde{\mathcal{X}}_{\mathbf{k}}^{\mathcal{E}_{\alpha}k_{\gamma}\mathcal{E}_{\beta}}+\tilde{\mathcal{W}}_{\mathbf{k}}^{\alpha,\beta\gamma}+\big{[}\tilde{\mathcal{W}}_{\mathbf{k}}^{\beta,\alpha\gamma}\big{]}^{*}\Big{)}\\\
=&-\frac{4\pi}{\Omega}2s\text{
Im}\sum_{m,l}^{o,e}\Bigg{\\{}i\frac{\bra*{m}r_{\alpha}\hat{Q}r_{\gamma}\ket*{l}\bra*{l}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})}-i\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}r_{\gamma}\hat{Q}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})}\\\
&\hskip
56.9055pt+\frac{i}{2}\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}\big{(}r_{\beta}r_{\gamma}-2r_{\gamma}\hat{P}r_{\beta}\big{)}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})}+\frac{1}{2}\frac{\bra*{m}r_{\alpha}\ket*{l}\bra*{l}\big{(}r_{\beta}v_{\gamma}-r_{\gamma}v_{\beta}\big{)}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})^{2}}\\\
&\hskip
56.9055pt-\frac{i}{2}\frac{\bra*{m}\big{(}r_{\alpha}r_{\gamma}-2r_{\alpha}\hat{P}r_{\gamma}\big{)}\ket*{l}\bra*{l}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})}+\frac{1}{2}\frac{\bra*{m}\big{(}r_{\alpha}v_{\gamma}-r_{\gamma}v_{\alpha}\big{)}\ket*{l}\bra*{l}r_{\beta}\ket*{m}}{(\epsilon_{l}-\epsilon_{m})^{2}}\Bigg{\\}}.\end{split}$
(53)

We now use $\hat{P}+\hat{Q}=1$ to simplify some terms,

$\begin{split}\tilde{\eta}_{\alpha\beta\gamma}=&-\frac{4\pi}{\Omega}2s\text{
Im}\sum_{m,l}^{o,e}\Bigg{\\{}\frac{1}{(\epsilon_{l}-\epsilon_{m})^{2}}\Big{[}\bra*{m}r_{\alpha}\ket*{l}\bra*{l}\big{(}r_{\beta}v_{\gamma}-r_{\gamma}v_{\beta}\big{)}\ket*{m}+\bra*{m}\big{(}r_{\alpha}v_{\gamma}-r_{\gamma}v_{\alpha}\big{)}\ket*{l}\bra*{l}r_{\beta}\ket*{m}\Big{]}\\\
&\hskip
73.97733pt+\frac{i}{(\epsilon_{l}-\epsilon_{m})}\Big{[}\bra*{m}r_{\alpha}r_{\gamma}\ket*{l}\bra*{l}r_{\beta}\ket*{m}-\bra*{m}r_{\alpha}\ket*{l}\bra*{l}r_{\gamma}r_{\beta}\ket*{m}\Big{]}\Bigg{\\}}.\end{split}$
(54)

As a final step, we use the fact that $\text{Im}\big{[}a-ib\big{]}=-b$ and
$\text{Im}\big{[}i(a-ib)\big{]}=a$, in order to obtain the following
expression for the natural optical activity tensor for finite samples,

$\begin{split}\tilde{\eta}_{\alpha\beta\gamma}&=-\frac{4\pi}{\Omega}s\text{
Im}\sum_{m,l}^{o,e}\Bigg{[}\frac{1}{(\epsilon_{l}-\epsilon_{m})^{2}}\Big{(}\bra*{m}r_{\alpha}\ket*{l}\bra*{l}(r_{\beta}v_{\gamma}-r_{\gamma}v_{\beta})\ket*{m}-\bra*{m}r_{\beta}\ket*{l}\bra*{l}(r_{\alpha}v_{\gamma}-r_{\gamma}v_{\alpha})\ket*{m}\Big{)}\\\
&\hskip
71.13188pt+\frac{i}{(\epsilon_{l}-\epsilon_{m})}\Big{(}\bra*{m}r_{\beta}\ket*{l}\bra*{l}r_{\alpha}r_{\gamma}\ket*{m}-\bra*{m}r_{\alpha}\ket*{l}\bra*{l}r_{\beta}r_{\gamma}\ket*{m}\Big{)}\Bigg{]},\end{split}$
(55)

which is the same expression that can be obtained from Ref. [10].

##  VIII VII. The two-rank gyration tensor and the optical rotatory parameter
for molecules

The optical rotatory parameter for finite systems like molecules, $\beta$, is
usually expressed as the sum of the diagonal components of the two-rank
gyration tensor [11, 12], $g_{\alpha\beta}$,

$\beta=\frac{\Omega}{4\pi}\frac{1}{2}\sum_{a}\frac{1}{3}g_{aa}.$ (56)

In order to go from our three-rank optical activity tensor,
$\tilde{\eta}_{\alpha\beta\gamma}$, to the two-rank gyration tensor,
$g_{\alpha\beta}$, all we need to do is to contract one index with the Levi-
Civita symbol (see Eq. (2) of the main text),

$\begin{split}g_{\alpha\beta}=&\frac{1}{2}\epsilon_{ab\alpha}\tilde{\eta}_{ab\beta}\\\
=&-\frac{1}{2}s\text{
Im}\sum_{m,l}^{o,e}\epsilon_{ab\alpha}\Bigg{[}\frac{1}{(\epsilon_{l}-\epsilon_{m})^{2}}\bra*{m}r_{a}\ket*{l}\bra*{l}\big{(}r_{b}v_{\beta}-r_{\beta}v_{b}\big{)}\ket*{m}+\frac{i}{(\epsilon_{l}-\epsilon_{m})}\bra*{m}r_{a}\ket*{l}\bra*{l}r_{b}v_{\beta}\ket*{m}\Bigg{]}.\end{split}$
(57)

One can then easily obtain the rotatory parameter given by Eq. (56),

$\begin{split}\beta=&-\frac{1}{3}s\text{
Im}\sum_{m,l}^{o,e}\frac{1}{(\epsilon_{l}-\epsilon_{m})^{2}}\Bigg{[}\bra*{m}r_{1}\ket*{l}\bra*{l}\big{(}r_{2}v_{3}-r_{3}v_{2}\big{)}\ket*{m}+\bra*{m}r_{2}\ket*{l}\bra*{l}\big{(}r_{3}v_{1}-r_{1}v_{3}\big{)}\ket*{m}\\\
&\hskip
128.0374pt+\bra*{m}r_{3}\ket*{l}\bra*{l}\big{(}r_{1}v_{2}-r_{2}v_{1}\big{)}\ket*{m}\Bigg{]}.\end{split}$
(58)

## References

  * Perdew and Wang [1992] J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
  * VAN [2018] Comput. Phys. Comm. 226, 39 (2018).
  * Hamann [2013] D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  * Asendorf [1957] R. H. Asendorf, J. Chem. Phys. 27, 11 (1957).
  * Pagola  _et al._ [2014] G. I. Pagola, M. B. Ferraro, P. F. Provasi, S. Pelloni, and P. Lazzeretti, J. Chem. Phys. 141, 094305 (2014).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Royo and Stengel [2019] M. Royo and M. Stengel, Phys. Rev. X 9, 021050 (2019).
  * Zabalo  _et al._ [2022] A. Zabalo, C. E. Dreyer, and M. Stengel, Phys. Rev. B 105, 094305 (2022).
  * Essin  _et al._ [2010] A. M. Essin, A. M. Turner, J. E. Moore, and D. Vanderbilt, Phys. Rev. B 81, 205104 (2010).
  * Malashevich and Souza [2010] A. Malashevich and I. Souza, Phys. Rev. B 82, 245118 (2010).
  * Krykunov and Autschbach [2006] M. Krykunov and J. Autschbach, J. Chem. Phys. 125, 034102 (2006).
  * Buckingham and Dunn [1971] A. D. Buckingham and M. B. Dunn, J. Chem. Soc. A , 1988 (1971).