#  Theory of optical axion electrodynamics and application to the Kerr effect
in topological antiferromagnets

Junyeong Ahn<EMAIL_ADDRESS>Department of Physics, Harvard
University, Cambridge, MA 02138, USA  Su-Yang Xu<EMAIL_ADDRESS>Department of Chemistry and Chemical Biology, Harvard University, Cambridge,
MA 02138, USA  Ashvin Vishwanath<EMAIL_ADDRESS>Department of
Physics, Harvard University, Cambridge, MA 02138, USA

###### Abstract

Emergent axion electrodynamics in magneto-electric media is expected to
provide novel ways to detect and control material properties with
electromagnetic fields. However, despite being studied intensively for over a
decade, its theoretical understanding remains mostly confined to the static
limit. Here, we introduce a theory of axion electrodynamics at general
frequencies. We define a proper optical axion magneto-electric coupling
through its relation to optical surface Hall conductivity and provide ways to
calculate it in lattice systems. By employing our formulas, we show that axion
electrodynamics can lead to a significant Kerr effect in thin-film
antiferromagnets at wavelengths that are seemingly too long to resolve the
spatial modulation of magnetism. We identify the wavelength scale above which
the Kerr effect is suppressed. Our theory is particularly relevant to
materials like MnBi2Te4, a topological antiferromagnet whose magneto-electric
response is shown here to be dominated by the axion contribution even at
optical frequencies.

##  I Introduction

In a medium where spatial inversion and time reversal symmetries are both
broken, magnetic and electric fields couple in a way that is fundamentally
different from the electromagnetic induction described by Maxwell’s equations
in vacuum, thus leading to exotic electrodynamics. A topic of particular
interest in magneto-electric coupling phenomena recently is axion
electrodynamics. Axion electrodynamics is theoretically proposed to be
realized in a class of topological materials called axion insulators [1, 2, 3,
4, 5, 6] for sufficiently slowly oscillating electromagnetic fields. In
particular, the static axion magneto-electric coupling is quantized by a
multiple of the fundamental constant $e^{2}/2h$ [7, 8, 9], which originates
from the half-quantized Hall conductance of the topological surface states.
While such a quantized magneto-electric coupling has not yet been directly
observed, experimental progress has been made to observe its consequences in
the static limit [10, 11, 12, 13, 14, 15, 16].

In contrast to the static limit, the emergent axion electrodynamics at generic
optical frequencies remains largely unexplored theoretically. Formulating a
theory of optical axion electrodynamics has to overcome the challenge that the
axion coupling is ill-defined in periodic three-dimensional systems, making it
hard to calculate and understand. The static axion angle can be calculated by
the Chern-Simons integral of the non-abelian Berry connection [8, 17, 18] at
the cost of being gauge dependent. Despite its gauge dependence, the Chern-
Simons integral is well-defined as an angular variable taking a value between
$0$ and $2\pi$, because a gauge transformation changes its value only by a
multiple of $2\pi$. This is a manifestation that the static surface Hall
conductance changes by a multiple of $e^{2}/h$ under deformations. However, at
optical frequencies, a similar approach does not seem feasible. This is
because the optical surface Hall conductivity is frequency dependent, and thus
a generic surface deformation changes the surface Hall response by a non-
quantized amount. Owing to this difficulty, theoretical understanding of the
optical axion electrodynamics has remained elusive to date [1].

In this paper, we make two important steps toward the complete formulation of
the optical axion electrodynamics. First, we show that a proper definition of
the optical magneto-electric coupling allows us to calculate the optical axion
angle in a fully gauge-independent way in a system with finite thickness.
Although the optical axion electrodynamics is a part of the linear-response
optical magneto-electric effects, it is distinguished from the other
contributions. Non-axionic magneto-electric effects are described within the
well-established theory of gyrotropic birefringence and natural optical
activity, which are based on the bulk current response [20, 21, 1]. On the
other hand, the optical axion electrodynamics is a surface phenomenon and thus
is not captured by those theories [1]. We therefore define the well-define
optical axion angle by analyzing the surface response. Using the well-defined
axion angle, we can understand the optical axion electrodynamics of thin films
and also that of bulk crystals by increasing the thickness. Second, in systems
with periodic boundary conditions along all three directions, we find the
optical axion angle at high optical frequencies can be estimated by the
optical layer Hall conductivity that we define in the main text. These
findings present advantages in numerical calculations as well as conceptual
advances.

Our development of a theory of optical axion electrodynamics opens the door to
understanding novel axion-induced optical phenomena. Here we provide a
concrete example. Magneto-optic effects make electromagnetic waves powerful
probes of the magnetic structure in materials. Conventionally, both Kerr and
Faraday effects are attributed to the optical Hall effect and thus commonly
perceived as probes of net magnetization [22, 23]. Although not well known,
previous theoretical studies proposed that the magneto-electric effect can
also lead to the Kerr effect, providing a novel way to probe fully compensated
antiferromagnets whose symmetries strictly prohibit the Hall effect and net
magnetization [20, 24, 25, 26, 27, 28]. Still, however, there are two aspects
that need further investigation. First, axion electrodynamic contribution to
the Kerr effect has not been well understood. Second, since previous studies
have focused on three-dimensional bulk systems, magneto-optic Kerr effects in
quasi two-dimensional antiferromagnets remain elusive. We present theoretical
analysis revealing the precise conditions for realizing the optical axion
electrodynamic Kerr effect as well as quantitative numerical analysis allowed
by our gauge-invariant formulas.

Our results apply broadly to both topological and non-topological media
because optical magneto-electric effects, as non-quantized phenomena, are not
sensitive to the topological nature of the ground state. Therefore, our work
provides a theoretical basis for the detection and manipulation of
antiferromagnetism in a large class of materials, thus having potential for
wide applications to antiferromagnetic spintronics and the study of magnetic
structure in quantum materials

Meanwhile, topological antiferromagnets need special attention as they are
ideal platforms for optical axion electrodynamics. To understand this, we note
two aspects. First, it is desired to have antiferromagnets having bulk
symmetry that reverses an odd number of spacetime coordinates and is broken on
the surfaces (e.g. spatial inversion symmetry), because such a symmetry
suppresses bulk magneto-electric effects but not the axion magneto-electric
effect. This condition is similar to the requirement of the quantization of
the static axion angle but additionally requires that the quantizing-symmetry
is broken at the surface to allow for a nonzero value. Second, spatially
spreading of electronic quantum states is needed, in order to show a response
distinguished from that of decoupled layered or Mott antiferromagnets. This
condition is again satisfied in topological antiferromagnets. We thus apply
our theory to a model of MnBi2Te4, which is the only experimentally realized
axion topological antiferromagnet to date, and show that the Kerr effect in
this material is significant.

Figure 1:  Optical magneto-electric effect. Magneto-electric (ME) coupling
leads to different electrodynamics within the bulk and on the surface. Inside
the magneto-electric medium, two linearly polarized light propagate with
different speeds because the wave equation is modified by origin-independent
magneto-electric coupling $T_{ij}$ and electric quadrupole susceptibility
$S_{ijk}$ (this effect is called gyrotropic birefringence [20]). On the
surface, axion magneto-electric coupling $\theta^{(z)}$ comes into play
additionally, contributing to the surface Hall conductivity. This is most
readily seen by writing the action for optical axion electrodynamics allowing
for a spatially varying $\theta$ parameter, $S_{\rm OA}=(e^{2}/2\pi
h)\int\theta(x)\bf E\cdot B$. Integrating by parts in the presence of an
interface along the $z$ direction where the axion angle jumps at the interface
$\delta\theta$ and identifying the coefficient of $A_{i}$ with the surface
current density we find ${\bf J}=-(e^{2}/2\pi h)\delta\theta\hat{z}\times{\bf
E}$. Thus an electric field parallel to the surface sets up a current in the
orthogonal direction along the surface, indicative of a surface Hall effect.
$S_{\rm OA}$ does not modify the propagation of electromagnetic fields within
the bulk medium where $\theta$ does not vary spatially.

##  II Results

###  II.1 Optical magneto-electric coupling

Motivated by the equivalence between the surface Hall conductivity and the
axion angle in the static limit, we study the surface current response to
define the optical axion angle. We consider currents generated by
electromagnetic multipole moments. Electromagnetic responses from multipole
moments are smaller for higher order moments [21]: electric dipole $P_{i}$
$\gg$ electric quadrupole $Q_{ij}$ and magnetic dipole $M_{i}$ $\gg$ higher
orders. Here, we consider only up to the electric quadrupole-magnetic dipole
order, giving the leading-order magneto-electric effect. The bulk current
density is related to the multipole moments by:
$J_{i}=\dot{P}_{i}-\frac{1}{2}\partial_{j}\dot{Q}_{ij}+\epsilon_{ijk}\partial_{j}M_{k}$.
While electric quadrupole and magnetic dipole moments do not generate
macroscopic currents in macroscopically homogeneous lattice systems, they
generate currents on the system boundary where the material property changes
spatially. The induced multipole moments have the following form in the
frequency domain [21]:

$\displaystyle P_{i}$
$\displaystyle=\sum_{j}(\chi_{ij}-i\chi_{ij}^{\prime})E_{j}+\frac{1}{2}\sum_{jk}(a_{ijk}-ia^{\prime}_{ijk})\nabla_{k}E_{j}$
$\displaystyle\quad+\sum_{j}(G_{ij}-iG^{\prime}_{ij})B_{j},$ $\displaystyle
Q_{ij}$ $\displaystyle=\sum_{k}(a_{kij}+ia^{\prime}_{kij})E_{k},$
$\displaystyle M_{i}$ $\displaystyle=\sum_{j}(G_{ji}+iG^{\prime}_{ji})E_{j},$
(1)

where $\chi_{ij}=\chi_{ji}$ and $\chi^{\prime}_{ij}=-\chi^{\prime}_{ji}$ are
the electric susceptibility tensors, $G_{ij}$ and $G^{\prime}_{ij}$ are
magneto-electric coupling, and $a^{\prime}_{ijk}=a^{\prime}_{ikj}$ and
$a_{ijk}=a_{ikj}$ are electric quadrupolar susceptibility tensors. Here, we
are interested in the magnetic magneto-electric effects described by $G_{ij}$
and $a^{\prime}_{ijk}$, which occur only when time reversal symmetry is broken
while being compatible with spacetime inversion $PT$ symmetry. Therefore, we
assume $PT$ symmetry to neglect the complications arising from the bulk Hall
conductivity and natural optical activity, both of which are excluded in this
symmetry setting, i.e., $\chi^{\prime}_{ij}=0$ and $G^{\prime}_{ij}=a_{ijk}=0$
[Table. 1]. This assumption does not affect our key results (see Appendix A
for discussions without $PT$ symmetry). As we show below, magneto-electric
effects occur in combination with electric quadrupole responses at optical
frequencies, requiring the consideration of the combination of $G_{ij}$ and
$a^{\prime}_{ijk}$.

Tensor | $P$ | $T$ | $PT$ | Phenomena | Kerr | Faraday
---|---|---|---|---|---|---
$\chi_{ij}$ | $+$ | $+$ | $+$ | Refraction and absorption | No | No
$\chi^{\prime}_{ij}$ | $+$ | $-$ | $-$ | Hall effect | Yes | Yes
$G_{ij},a^{\prime}_{ijk}$ | $-$ | $-$ | $+$ | Optical magneto-electric effect | Yes | No
$G^{\prime}_{ij},a_{ijk}$ | $-$ | $+$ | $-$ | Natural optical activity | No | Yes
Table 1:  Symmetry properties of electromagnetic linear response functions.
Response functions are defined in Eq. (II.1). The sign in the second column
shows that parity of the response functions under spatial inversion $P$, time
reversal $T$, and spacetime inversion $PT$. Kerr and Faraday in the last
column indicate the optical rotation of the light polarization plane in
reflection and transmission, respectively (i.e., the Kerr and Faraday
effects). Kerr and Faraday effects are allowed when $T$ and $PT$ symmetry are
broken, respectively [29, 30, 31].

Let us consider the surface at $z=0$ of a three-dimensional homogeneous
material with the outward normal direction $\hat{z}$ [Fig. 1]. The surface
current density ${\bf j}^{s}=\int^{d_{s}/2}_{-d_{s}/2}dz{\bf J}$, where
$d_{s}$ is the characteristic thickness of the interface where response
functions change rapidly as functions of $z$, is related to the multipole
moments through

$\displaystyle j^{s}_{i}(\omega)$ $\displaystyle=\sum_{j,k}\left[\delta
G_{jk}(\omega)+\sum_{l}\frac{\omega}{2}\epsilon_{kzl}\delta
a^{\prime}_{jlz}(\omega)\right]\epsilon_{kzi}E_{j}(\omega)$
$\displaystyle\quad+\sum_{j}\left[\bar{\sigma}_{ij}(\omega)-i\bar{\sigma}_{izj}(\omega)\right]\delta
E_{j}(\omega),$ (2)

where $\delta f=f(-d_{s}/2)-f(d_{s}/2)$ and
$\bar{f}=[f(-d_{s}/2)+f(d_{s}/2)]/2$ are the difference and average of the
material property $f(z)$ across the interface, and $\sigma_{ij}$ and
$\sigma_{ijk}=i(\epsilon_{jkl}T_{il}+\epsilon_{ikl}T_{jl})+\omega S_{ijk}$ are
the bulk conductivity tensors defined by
$J_{i}=\sigma_{ij}E_{j}+\sum_{j,k}\sigma_{ijk}q_{j}E_{k}+O(q^{2})$ for light
wave vector $q$, where
$T_{ij}=G_{ij}-\frac{1}{3}\delta_{ij}\sum_{k=1}^{3}G_{kk}-\frac{i}{6}\omega\sum_{k,l=1}^{3}\epsilon_{jkl}a^{\prime}_{kli}$,
and $S_{ijk}=(a^{\prime}_{ijk}+a^{\prime}_{jki}+a^{\prime}_{kij})/3$ [21, 1].

We neglect the $\delta E$ term because it is a bulk response whose
contribution to the interface vanishes as $d_{s}\rightarrow 0$. The form of
the surface current density in Eq. (II.1) suggests defining surface-sensitive
magneto-electric coupling by

$\displaystyle G^{(z)}_{ix}(\omega)$
$\displaystyle=G_{ix}(\omega)-\frac{1}{2}\omega a^{\prime}_{iyz}(\omega),$
$\displaystyle G^{(z)}_{iy}(\omega)$
$\displaystyle=G_{iy}(\omega)+\frac{1}{2}\omega a^{\prime}_{ixz}(\omega),$ (3)

where the superscript $(z)$ explicitly shows that these quantities depend on
the choice of the surface normal direction [$(z)$ and $(-z)$ are the same,
though]. The $zz$ component of this surface-sensitive magneto-electric
coupling is defined from the response of the surface charge. Using
$\rho=-\nabla\cdot{\bf P}+\frac{1}{2}\partial_{i}\partial_{j}Q_{ij}+\ldots$,
we obtain
$\rho^{s}=G_{zi}B_{i}+\frac{\omega}{2}(a^{\prime}_{yzx}-a^{\prime}_{xzy})B_{z}+\frac{\omega}{2}a^{\prime}_{zzy}B_{x}-\frac{\omega}{2}a^{\prime}_{zzx}B_{y}+\ldots$
where the ellipsis includes the electric dipole term $\chi_{zi}E_{i}$, the
symmetric $(\partial_{j}E_{k}+\partial_{k}E_{j})$ terms, and higher-order
multipole terms. From this, we define

$\displaystyle G^{(z)}_{zz}(\omega)$
$\displaystyle=G_{zz}(\omega)+\frac{1}{2}\omega\left[a^{\prime}_{yxz}(\omega)-a^{\prime}_{xyz}(\omega)\right]$
(4)

as well as $G^{(z)}_{zx}=G_{zx}+\omega a^{\prime}_{zzy}/2$ and
$G^{(z)}_{zy}=G_{zy}-\omega a^{\prime}_{zzx}/2$.

While various components of the magneto-electric coupling are related to the
surface optical conductivity, there is a unique component that manifests
itself only at the surface: the axion angle. We define the optical axion angle
by the trace part of $G^{(z)}_{ij}$.

$\displaystyle\theta^{(z)}(\omega)\equiv\pi\frac{2h}{e^{2}}\frac{1}{3}\sum_{i=1}^{3}G^{(z)}_{ii}(\omega).$
(5)

To see that this only manifests itself at the surface rather than the bulk, we
note that the bulk current response by the magneto-electric and electric
quadrupole susceptibilities is determined only through $T_{ij}$ and $S_{ijk}$.
Because of this, previous studies focusing on bulk magneto-electric effect did
not capture optical axion electrodynamics [25, 1]. See Supplementary Note 1
for more details on the bulk response. The surface-sensitive magneto-electric
coupling is fully determined by these bulk-response quantities and the axion
angle:

$\displaystyle G^{(z)}_{xx}(\omega)$ $\displaystyle=\frac{e^{2}}{2\pi
h}\theta^{(z)}(\omega)+T_{xx}(\omega)-\frac{\omega}{2}S_{xyz}(\omega),$
$\displaystyle G^{(z)}_{yy}(\omega)$ $\displaystyle=\frac{e^{2}}{2\pi
h}\theta^{(z)}(\omega)+T_{yy}(\omega)+\frac{\omega}{2}S_{xyz}(\omega),$
$\displaystyle G^{(z)}_{zz}(\omega)$ $\displaystyle=\frac{e^{2}}{2\pi
h}\theta^{(z)}(\omega)+T_{zz}(\omega),$ $\displaystyle G^{(z)}_{ij}(\omega)$
$\displaystyle=T_{ij}(\omega)-\frac{\omega}{2}S_{iik}(\omega)\epsilon_{jik},\text{
for }i\neq j\text{ and }j\neq 3.$ (6)

Since $T_{ij}$ is traceless, it does not contribute to the trace of
$G^{(z)}_{ij}$. Note that $G^{(z)}_{ij}(\omega)$ transforms as a tensor under
magnetic layer group actions but not under the full three-dimensional magnetic
space group actions.

###  II.2 Magneto-electric coupling with open boundaries in one direction

Defining the combination $G^{(z)}_{ij}$ of $G_{ij}$ and $a^{\prime}_{ijk}$ has
an advantage in practical calculations as well as in the formulation. As
$G^{(z)}_{ij}$ characterizes the surface current response which is measurable,
it admits a gauge-invariant form in three-dimensional lattice systems with
open boundaries along one direction, or in quasi-two-dimensional systems. This
nice property is absent in the diagonal magneto-electric coupling $G_{ii}$,
making it hard to calculate.

The bare magneto-electric coupling $G_{ij}=\partial P_{i}/\partial B_{j}$ and
$a^{\prime}_{ijk}=\partial Q_{jk}/\partial E_{i}$ at zero temperature have the
following form according to linear response theory (see Appendix B).

$\displaystyle G_{ij}(\omega)$
$\displaystyle=\frac{V}{\hbar}\sum_{n,m}\frac{f_{nm}}{\omega_{mn}-\omega}{\rm
Re}\braket{n}{\hat{P}_{i}}{m}\braket{m}{\hat{M}_{j}}{n},$ $\displaystyle
a^{\prime}_{ijk}(\omega)$
$\displaystyle=-\frac{V}{\hbar}\sum_{n,m}\frac{f_{nm}}{\omega_{mn}-\omega}\frac{\omega_{mn}}{\omega}{\rm
Im}\braket{n}{\hat{P}_{i}}{m}\braket{m}{\hat{Q}_{jk}}{n},$ (7)

where $V$ is the volume of the system, $n,m$ are indices for energy
eigenstates with eigenvalue $\hbar\omega_{n}$, $f_{nm}=f_{n}-f_{m}$ is the
difference between the Fermi-Dirac distribution of the $n$ and $m$ states.
$\hat{P}_{i}=-e\hat{r}^{i}/V$,
$\hat{M}_{j}=-(e\epsilon_{jkl}\hat{r}^{k}\hat{v}^{l}/2+\hat{m}^{s}_{j})/V$,
and $\hat{Q}_{jk}=-e\hat{r}^{j}\hat{r}^{k}/V$ are electric dipole, magnetic
dipole, and electric quadrupole density operators, respectively, where
$\hat{\bf r}$ and $\hat{\bf v}$ are the position and velocity operators of
electrons, and $\hat{m}^{s}$ is the spin magnetic moment operator. Equation
(II.2) can be calculated for molecular systems [21, 32], meaning finite
systems with open boundary conditions along all directions, by using the real-
space representation of $\hat{\bf r}$ and the relation $\hat{\bf
v}=-i\hbar^{-1}[\hat{\bf r},\hat{H}]$, where $\hat{H}$ is the Hamiltonian of
the system.

In periodic systems, however, $G_{ij}$ and $a^{\prime}_{ijk}$ are not well
defined separately because the position operator is not well-defined because
the position is not uniquely defined. This manifests through the momentum-
space representation of the position operator $\braket{\psi_{m\bf
k^{\prime}}}{\hat{r}}{\psi_{n\bf k}}=-\delta_{mn}i\partial_{\bf k}\delta_{\bf
k^{\prime}k}+\delta_{\bf k^{\prime}k}\braket{u_{m\bf k}}{i\partial_{\bf
k}}{u_{n\bf k}}$, whose diagonal matrix element is not well defined because of
$\partial_{\bf k}\delta_{\bf k^{\prime}k}$. On the other hand, the diagonal
matrix elements of the position operator do not appear in the response
functions that have a well defined physical meaning in periodic systems.
$T_{ij}$ and $S_{ijk}$ are such examples that characterizes the bulk current
response [1]. Since $G^{(z)}$ characterizes the surface response, one can
expect that it is well defined in quasi-two-dimensional periodic systems.

After doing some algebra that we relegate to Appendix B, we can write
$G^{(z)}_{ii}$s in two-dimensional momentum space as

$\displaystyle G^{(z)}_{xx}(\omega)$ $\displaystyle=\frac{e^{2}}{\hbar
V}\sum_{n\neq m,k_{x},k_{y}}\frac{f_{nm}}{\omega_{mn}-\omega}{\rm
Re}\left[r^{x}_{nm}\braket{\psi_{m(k_{x},k_{y})}}{-\frac{1}{2}\left(\hat{v}^{y}\hat{r}^{z}+\hat{r}^{z}\hat{v}^{y}\right)+\hat{m}^{\rm
s}_{x}}{\psi_{n(k_{x},k_{y})}}\right],$ $\displaystyle G^{(z)}_{yy}(\omega)$
$\displaystyle=\frac{e^{2}}{\hbar V}\sum_{n\neq
m,k_{x},k_{y}}\frac{f_{nm}}{\omega_{mn}-\omega}{\rm
Re}\left[r^{y}_{nm}\braket{\psi_{m(k_{x},k_{y})}}{\frac{1}{2}\left(\hat{v}^{x}\hat{r}^{z}+\hat{r}^{z}\hat{v}^{x}\right)+\hat{m}^{\rm
s}_{y}}{\psi_{n(k_{x},k_{y})}}\right],$ $\displaystyle G^{(z)}_{zz}(\omega)$
$\displaystyle=\frac{e^{2}}{\hbar V}\sum_{n\neq
m,k_{x},k_{y}}\frac{f_{nm}}{\omega_{mn}-\omega}{\rm
Re}\left[\frac{1}{2}\sum_{p;E_{p}\neq
E_{m}}\left(r^{z}_{nm}r^{x}_{mp}v^{y}_{pn}-r^{z}_{np}r^{x}_{pm}v^{y}_{mn}-(x\leftrightarrow
y)\right)+r^{z}_{nm}(m^{\rm s}_{z})_{mn}\right],$ (8)

where $r^{i}_{mn}=\braket{\psi_{m}}{\hat{r}^{i}}{\psi_{n}}$ and
$v^{i}_{mn}=\braket{\psi_{m}}{\hat{v}^{i}}{\psi_{n}}$ are matrix elements of
the position and velocity operators, $|\psi_{n(k_{x},k_{y})}\rangle$ is the
Bloch state, the subscripts $n$ ,$m$ and $p$ are band indices. While the
position operator matrix element is not well defined in momentum space [33],
the diagonal matrix elements of $\hat{r}^{x}$ and $\hat{r}^{y}$ do not appear
in Eq. (II.2), so that the surface-sensitive magneto-electric coupling is well
defined in two-dimensional momentum space. Combining equations in Eqs. (5) and
(II.2), we can obtain the optical axion angle. In simple tight-binding models
where $\hat{r}^{z}$ commutes with $\hat{v}^{y}$, $G^{(z)}_{xx}(0)$ reduces to
the expression of $G_{xx}(0)$ derived by Liu and Wang [34]. To our knowledge,
our formulas represent the first expressions for calculating the diagonal
components of the optical magneto-electric coupling in crystalline (periodic)
systems.

The form of $G^{(z)}_{xx}$ and $G^{(z)}_{yy}$ suggests that their orbital
magneto-electric part may be interpreted as a real-space dipole of the Berry
curvature. This idea works exactly for two-band systems in the limit of
decoupled layer systems, where the matrix element part can be written as the
product of the Berry curvature $F_{xy}=-2{\rm Im}[r^{x}_{12}r^{y}_{21}]$ and
the $z$ component of the position eigenvalues $r^{z}_{11}$ or $r^{z}_{22}$. We
apply this idea below to the case where the $z$ direction is also periodic.

Figure 2:  Intra-cell and unit-cell-position contributions to optical magneto-
electric coupling a, Preferred choices of the unit cell with open boundaries.
Arrows represent the average moment of a layer, where we bipartition the
degrees of freedom within the unit cell into the upper and lower parts, each
of which defines what we call a layer. The left-hand-side and middle even-
layer systems are related by spatial inversion symmetry and also by the
translation of the blue-arrow layers by a lattice constant if the deformation
of the electronic states at the surface is neglected. The right-hand-side odd-
layer system is inversion symmetric. The inversion centers are shown as yellow
stars. b, Hall conductivity of the unit cell in the presence of surface
states. The vertical displacement $R_{z}$ of a unit cell times its net Hall
conductivity contributes to diagonal magneto-electric coupling $G^{(z)}_{xx}$,
$G^{(z)}_{yy}$, and $G^{(z)}_{zz}$.

###  II.3 Intra-cell magneto-electric coupling and optical layer Hall
conductivity

In fully periodic lattice systems where the $z$ direction is also periodic, it
is hard to calculate the full orbital part of $G^{(z)}_{ii}(\omega)$ because
$\hat{r}^{z}$ is not well defined in momentum space. Nevertheless, we can
still define and calculate the magneto-electric coupling of the unit cell,
which we call the intra-cell magneto-electric coupling. Note, this treatment
will be necessarily approximate, in contrast to our previous discussion of the
slab geometry, but provides an physical understanding for the results. For
example, the use of the intra-cell magneto-electric coupling makes the
relation concrete between the axion angle and the antiferromagnetism in
inversion symmetric systems. Ultimately we must compare the results between
this approach and the previous slab calculation as we do in another section
below.

Let us begin by giving a physical intuition that a nonzero axion angle is
natural in a fully compensated antiferromagnet [35]. To see this, recall the
analogy between electric polarization in 1D and the axion theta angle $\theta$
in 3D. The former can be defined in a system free of net charge, by stacking
alternate positive and negative charges. Similarly, if we stack alternate
planes with Hall conductance $\pm g^{H}_{xy}$ such that the net Hall
conductance vanishes, the axion magneto-electric coupling becomes well defined
and is the analog of electric polarization of the alternating planes. In fact,
with an applied magnetic field perpendicular to the planes, the induced
polarization from having alternating charges in the $\pm g^{H}_{xy}$ layers,
leads to a finite electric polarization which is readily calculated as
$g^{H}_{xy}\delta a_{z}/a_{z}$, where $\delta a_{z}$ is spacing between
alternating antiferromagnetic planes versus the vertical size of the unit cell
$a_{z}$. The charge in each flux quantum area is $g_{H}(h/e)$. For an
antiferromagnet where spacing between planes $=1/2a_{z}$, this is just
$\theta=2\pi g^{H}_{xy}/2$.

To present a more detailed analysis, let us decompose the position operator
into the intra-cell polarization and unit cell position parts by
$r^{z}_{mn}={\mathbb{A}}^{z}_{mn}+R^{z}_{mn}$. Namely,

$\displaystyle\braket{\psi_{m\bf k^{\prime}}}{\hat{r}^{z}}{\psi_{n\bf k}}$
$\displaystyle=\delta_{\bf k,k^{\prime}}\sum_{\beta,\alpha}\braket{\psi_{m\bf
k}}{\psi_{\beta{\bf k}}}{\mathbb{A}}^{z}_{\beta\alpha}({\bf
k})\braket{\psi_{\alpha{\bf k}}}{\psi_{n\bf k}}$
$\displaystyle\quad+\sum_{\alpha,{\bf R}}\braket{\psi_{m\bf
k^{\prime}}}{w_{\alpha{\bf R}}}R^{z}\braket{w_{\alpha{\bf R}}}{\psi_{n\bf k}}$
(9)

where $|w_{\alpha{\bf R}}\rangle$ is the Wannier state with the collective
index $\alpha$ for spin and orbital (cf. $n$ and $m$ are band indices),
$|\psi_{\alpha{\bf k}}\rangle=N^{-1/2}\sum_{\bf R}e^{i{\bf k}\cdot{\bf
R}}|w_{\alpha{\bf R}}\rangle$, and ${\mathbb{A}}^{z}_{\beta\alpha}({\bf
k})=\sum_{\bf R}\braket{w_{\beta{\bf 0}}}{\hat{r}^{z}}{w_{\alpha{\bf
R}}}e^{i{\bf k}\cdot{\bf R}}$. The second term in Eq. (II.3) vanishes for
$n\neq m$ as one can see by writing it as
$-i\delta_{mn}\partial_{k^{z}}\delta_{\bf k,k^{\prime}}$. This decomposition
is independent of choosing the basis $\alpha$ within the unit cell, where each
unit cell is labelled by a given lattice vector ${\bf R}$. However, the
decomposition depends on the choice of the unit cell, which we discuss below.

The intra-cell polarization term, the first term in Eq. (II.3), defines the
magneto-electric coupling of the unit cell. This intra-cell magneto-electric
coupling depends on the choice of the unit cell. In our case, choosing a unit
cell corresponds to fixing the value of the Wannier centers
$\braket{w_{\alpha{\bf 0}}}{\hat{r}^{z}}{w_{\alpha{\bf 0}}}$, which is
ambiguous by respective lattice translations of the Wannier states
$|w_{\alpha{\bf 0}}\rangle$. A physical interpretation of this multi-
valuedness is similar to that of electric polarization [35]: The magneto-
electric coupling depends on how we open the boundary, and there exists a
preferred choice of the unit cell for each boundary condition [Fig. 2(a)].

The unit-cell position term has the form
$\sum_{R_{z}}R^{z}\sigma^{H}_{xy}(R_{z})$, where
$\sigma^{H}_{xy}\equiv(\sigma_{xy}-\sigma_{yx})/2$ is the Hall conductivity,
and $\sigma^{H}_{xy}(R^{z})=-e^{2}\hbar^{-1}V^{-1}\sum_{n\neq
m,k_{x},k_{y}}f_{nm}\omega_{mn}/(\omega_{mn}-\omega)^{-1}{\rm
Im}\left[r^{x}_{nm}r^{y}_{mn}P^{R^{z}}_{nn}\right]$, and
$\hat{P}^{R^{z}}=\sum_{\alpha,R_{x},R_{y}}|w_{\alpha{\bf R}}\rangle\langle
w_{\alpha{\bf R}}|$ is the projection to $R^{z}$. This term is also multi-
valued in periodic systems because the unit cell position ${\bf R}$ is not
uniquely defined in periodic systems. This part, however, does not contribute
to the magneto-electric coupling when the Hall conductivity of the unit layer
vanishes, which is the case in $PT$-symmetric systems. When the boundary is
introduced, the Hall conductivity of the surface unit layer can be nonzero
even when the bulk Hall conductivity vanishes [Fig. 2(b)]. This change of the
$R^{z}$ term is the main source of the difference in the magneto-electric
couplings between finite-size systems and periodic lattice systems. This
effect is significant especially in the static limit of axion insulators,
where the emergent Dirac surface states modify the low-frequency surface Hall
conductivity in the order of $e^{2}/2h$.

As we derive below, in inversion-symmetric even-layer antiferromagnets, the
intra-cell contribution to the magneto-electric coupling is simplified to the
optical layer Hall conductivity. Namely, the optical axion angle is

$\displaystyle\frac{e^{2}}{2\pi
h}\theta^{(z)}(\omega)\approx\frac{1}{2}a_{z}\sigma^{\rm
LH}_{xy}(\omega)+\sum_{R_{z}}R^{z}\sigma^{H}_{xy}(R_{z})$ (10)

for inversion-symmetric AFMs, where $a_{z}$ is the vertical lattice constant,
and $\sigma_{xy}$ is a three-dimensional conductivity taking the unit of
conductance (unit of $e^{2}/h$) divided by the length. As for the definition
of layers, we note that there are two inversion-invariant values of the
vertical displacement $z$ within $a_{z}$. We refer to the two quasi-two-
dimensional bipartite regions centered at those two invariant $z$ values as
layers [Fig. 2(a)]. The layer Hall conductivity is then defined as
$\sigma^{\rm LH}_{xy}=(\sigma^{H,u}_{xy}-\sigma^{H,d}_{xy})/2$ from the bulk
optical Hall conductivity projected to the single layer $l=u,d$:

$\displaystyle\sigma^{H,l}_{xy}(\omega)$ $\displaystyle=-\frac{e^{2}}{\hbar
V}\sum_{n\neq m,k_{x},k_{y}}\frac{f_{nm}\omega_{mn}}{\omega_{mn}-\omega}{\rm
Im}\left[r^{x}_{nm}r^{y}_{mn^{\prime}}P^{l}_{n^{\prime}n}\right],$ (11)

where $P^{l}_{n^{\prime}n}({\bf k})=\sum_{\alpha,{\bf R};r_{\alpha,{\bf
R}}^{z}\in l\text{ layers}}\braket{\psi_{n\bf k}}{w_{\alpha{\bf
R}}}\braket{w_{\alpha{\bf R}}}{\psi_{n^{\prime}\bf k}}$, $|w_{\alpha{\bf
R}}\rangle$s are inversion-symmetric Wannier states, and $r^{z}_{\alpha,{\bf
R}}=\braket{w_{\alpha{\bf R}}}{\hat{r}^{z}}{w_{\alpha{\bf R}}}$ is the Wannier
center.

To obtain Eq. (10), note that the unit cell for even-layer systems is
symmetric under the combination of inversion and a lattice translation of
either the $u$ layer by $-a_{z}$ (or the $l$ layer by $+a_{z}$) [Fig. 2(a)].
If we focus on the intra-cell part (the ${\mathbb{A}}$ part), the combined
symmetry gives a constraint
$\theta^{(z)}_{\mathbb{A}}(\omega)=-\theta_{\mathbb{A}}^{(z)}(\omega)-a_{z}\sigma^{H,d}_{xy}(\omega)$,
where the last term is due to the transformation property of the axion angle
$\delta\theta^{(z)}(\omega)=d_{z}\sigma^{H}_{xy}(\omega)$ under the
translation by a vector ${\bf d}$. As we consider antiferromagnets with zero
net Hall conductivity, such that
$\sigma^{H,d}_{xy}=-\sigma^{H,u}_{xy}=-\sigma^{\rm LH}_{xy}$, we arrive at Eq.
(10). Another useful way of understanding this result is to think of the
electric polarization generated by applying a uniform magnetic field, and then
the displacement of layers of one sign of the Hall effect obviously
contributes to change in electric polarization, as we explained at the
beginning of this section.

Figure 3:  Reflections by $PT$-symmetric magneto-electric media. We assume
$PT$ symmetry to exclude natural optical activity, for simplicity. a, Single
interface. The light propagation within each medium is determined by the
refractive index $n$ and origin-independent magneto-electric coupling $T_{ij}$
and electric quadrupole susceptibility $S_{ijk}$ tensors. We consider the case
where no birefringence occurs such that the light propagation is described by
a unique refractive index in each medium. b, Double interfaces. Infinite
reflections occur between the top and bottom surfaces of medium 2. Light
obtains the complex phase $\phi=kd=n_{2}\omega d/c$ while propagating the
distance $d$ in medium 2.

###  II.4 Magneto-optic effects in fully compensated antiferromagnets

The optical axion angle manifests directly through the magneto-optic Kerr
effect. To gain some intuition for this, recall that the change in axion angle
at the surface leads to a surface Hall effect
$\sigma^{s}_{xy}=(\delta\theta/2\pi)e^{2}/h$. Clearly, such a surface Hall
conductance will lead to a Kerr effect. We present a more systematic analysis
below.

Let us consider the reflection at the single interface between two media $1$
and $2$ [Fig. 3]. For simplicity, we assume that both media have $M_{x}T$,
$C_{3z}$, and $PT$ symmetries, which are shared by magneto-electric materials
Cr2O3 and MnBi2Te4. We also assume normal incidence (i.e., light is incident
along $-\hat{z}$). Then there is no birefringence because of the symmetry we
require, so the propagation of light within each medium is then characterized
by a single complex-valued refractive index $n_{\mu}$, where $\mu=1,2$ labels
the two media [See Eq. (II.1)].

The electric field in medium $1$ consists of incident and reflected fields
${\bf E}_{i}$ and ${\bf E}_{r}$ while that in medium $2$ is the transmitted
field ${\bf E}_{t}$:

$\displaystyle{\bf E}_{1}$ $\displaystyle={\bf E}^{i}+{\bf
E}^{r}\equiv(1+r){\bf E}^{i},$ $\displaystyle{\bf E}_{2}$ $\displaystyle={\bf
E}^{t}\equiv t{\bf E}^{i},$ (12)

where $r$ and $t=1+r$ are $2\times 2$ Jones matrices for reflection and
transmission, respectively. ${\bf E}_{1}={\bf E}_{2}$ at the interface because
electric fields are continuous by Faraday’s law, because we consider normal
incidence such that electric fields are parallel to the interface. The
contribution from the surface conductivity is encoded in the boundary
condition of the magnetic field at the interface.

$\displaystyle{\bf B}^{t}={\bf B}^{i}+{\bf B}^{r}+\mu_{0}\hat{z}\times{\bf
j}^{s}$ (13)

Here, we consider only the surface currents induced from the bulk, i.e.,
$\hat{z}\times{\bf j}^{s}=\sigma^{s}_{xy}{\bf E}$, where
$\sigma^{s}_{xy}=G^{(z),\mu=2}_{xx}-G^{(z),\mu=1}_{xx}$. By solving the
boundary condition equations, we obtain

$\displaystyle r_{xx}$
$\displaystyle=\frac{(n_{1}-n_{2})(n_{1}+n_{2})-(\mu_{0}c\sigma^{s}_{xy})^{2}}{(\mu_{0}c\sigma^{s}_{xy})^{2}+(n_{1}+n_{2})^{2}},$
$\displaystyle r_{xy}$
$\displaystyle=-\frac{2n_{1}(\mu_{0}c\sigma^{s}_{xy})}{(\mu_{0}c\sigma^{s}_{xy})^{2}+(n_{1}+n_{2})^{2}}.$
(14)

Magneto-electric coupling appears in the reflective Jones matrix through the
surface conductivity, which is a manifestation that the Kerr effects here are
surface phenomena.

The complex Kerr angle is defined by
$\phi_{K}=\varphi_{K}+i\eta_{K}=\tan^{-1}(r_{xy}/r_{xx})$. Its real part
measures the rotation of the light polarization plane while its imaginary part
measures the circular dichroism, i.e., the intensity imbalance between the
reflected left and right circularly polarized light. The Kerr angle is
$O(\mu_{0}c\sigma^{s}_{xy})$ in general, but it can be much enhanced when
$n_{1}\approx n_{2}$ because of the suppression of $r_{xx}$.

When the sample thickness is much larger than the wavelength of light, it is
enough to suppose that reflection occurs at a single interface. However, when
thickness $d$ is comparable to or less than wavelength $\lambda$, which is the
case particularly relevant for thin films, we need to consider the response
from the whole sample including top and bottom surfaces (cf. When the photon
energy is smaller than the band gap, double interfaces can be relevant even
for $d\gg\lambda$ for an insulating medium [14] because then light incident on
the top reaches the bottom without attenuation.).

To study this case, we consider three media with refractive index $n_{\mu}$,
where $\mu=1,2,3$, as shown in Fig. 3b. We again assume $M_{x}T$, $C_{3z}$,
and $PT$ symmetries for each media. The Jones reflection matrix of the sample
for light incident from medium 1 is then the infinite sum of multiple
reflections.

$\displaystyle r$
$\displaystyle=r_{T}+e^{2i\phi}(1+r^{\prime}_{T})r_{B}\left(1-e^{2i\phi}r^{\prime}_{T}r_{B}\right)^{-1}(1+r_{T}),$
(15)

where we use $t_{T,B}=1+r_{T,B}$ and $t^{\prime}_{T}=1+r^{\prime}_{T}$, and
$\phi=n_{2}\omega d/c$ is the complex-valued phase obtained by the one-way
propagation through the sample thickness $d$. Here, the subscripts $T$ and $B$
indicate the top and bottom of the sample, and the prime indicates the process
where the light is incident to the top surface from below (we follow the
notation in Ref. [13]). See Appendix D for the expression of Jones matrices
$r_{T}$, $r^{\prime}_{T}$, and $r_{B}$. Similarly, for transmission,

$\displaystyle
t=(1+r_{B})\left(1-e^{2i\phi}r^{\prime}_{T}r_{B}\right)^{-1}e^{i\phi}(1+r_{T}).$
(16)

Note that $t\neq 1+r$ when $\phi\neq 0$.

Considering the case where only medium 2 is magneto-electric, we set
$G^{(z),1}_{xx}=G^{(z),3}_{xx}=0$ and $G^{(z),2}_{xx}\equiv G^{(z)}_{xx}\neq
0$. We first consider $|\delta n|\ll|\phi|$, where $\delta n\equiv
n_{3}-n_{1}$. The Kerr and Faraday rotation angles are then given by

$\displaystyle\tan\phi_{K}$
$\displaystyle=\frac{2n_{1}\mu_{0}cG^{(z)}_{xx}}{n_{2}^{2}-n_{1}^{2}+(\mu_{0}cG^{(z)}_{xx})^{2}},$
$\displaystyle\tan\phi_{F}$
$\displaystyle=-\frac{\mu_{0}cG^{(z)}_{xx}\sin\phi}{\left(n_{1}^{2}+n_{2}^{2}+(\mu_{0}cG^{(z)}_{xx})^{2}\right)\sin\phi+2in_{1}n_{2}\cos\phi}\delta
n$ (17)

in the leading order of $\delta n$, where $\tan\phi_{F}=t_{xy}/t_{xx}$. A
nonzero Faraday rotation requires $\delta n\neq 0$ because $PT$ symmetry needs
to be broken [31]. On the other hand, the Kerr rotation can be nonzero with
$\delta n=0$ because it is compatible with $PT$ symmetry [Table 1]. The Kerr
rotation is independent of $\phi$ in the leading order of $\delta n$ when
$\delta n$ is sufficiently small as if the response comes from the top surface
only, while it actually comes from multiple reflections between the top and
bottom surfaces. In this limit, the way the Kerr effect goes away is highly
nontrivial as the thickness goes to zero (i.e., $|\phi|\rightarrow 0$ with
$|\phi|\gg|\delta n|$). The Kerr angle stays constant, but the amplitudes of
the reflected signals ultimately vanish. Therefore, a finite Kerr effect can
be observed from a thin film when optical equipment is highly sensitive.

However, a nontrivial $\phi$ dependence appears when $\phi$ is the smallest
parameter ($|\delta n|\gg|\phi|$), where we have

$\displaystyle\tan\phi_{K}$
$\displaystyle=\frac{4in_{1}n_{3}\mu_{0}cG^{(z)}_{xx}}{n_{2}(n_{1}^{2}-n_{3}^{2})}\phi+O\left(\phi^{2}\right),$
$\displaystyle\tan\phi_{F}$
$\displaystyle=-\frac{i(n_{1}-n_{3})\mu_{0}cG^{(z)}_{xx}}{n_{2}(n_{1}+n_{3})}\phi+O\left(\phi^{2}\right),$
(18)

which show that both Kerr and Faraday rotation vanish for $\phi=0$. Therefore,
a nonzero $\phi$ is necessary for the Kerr effect as well as the Faraday
effect in fully compensated antiferromagnets. Note that, while the Kerr angle
in Eq. (II.4) is independent of $\phi$, it applies only when $|\phi|\gg|\delta
n|$.

The crossover between two regimes respectively described by Eqs. (II.4) and
(II.4) occurs when $|\phi|\sim|\delta n|$. The corresponding wavelength scale

$\displaystyle\lambda\sim\lambda^{*}\equiv 2\pi\left|\frac{n_{2}}{\delta
n}\right|d$ (19)

is much larger than the sample thickness $d$ when $|\delta n|\ll 1$. This is
remarkable because one might naively expect that, because the response from
the spatially separated top and bottom layers are not resolved when
$\lambda\gg d$, the Kerr effect is vanishingly small in that regime and scales
linearly with $d/\lambda$. However, our analysis shows that a much stricter
$\lambda\gg\lambda^{*}$ is required for the suppression of the Kerr angle. To
explain how this works, we note again that the Kerr angle and the amplitude of
the Kerr rotated signal can behave differently. While the amplitude of the
Kerr rotated signal ($\propto r_{xy}$) is indeed suppressed for $\lambda\gg
d$, the amplitude of the non-rotated signal ($\propto r_{xx}$) is also
suppressed in the same limit. Their suppression at large wavelengths
compensate each other to keep the ratio $\phi_{K}=r_{xy}/r_{xx}$ as long as
$\lambda$ is smaller than a larger length scale $\lambda^{*}$, above which the
Kerr angle ultimately gets suppressed.

In thin film axion insulators, equation (II.4) is typically more relevant at
photon energies below the surface gap. The surface gap of experimentally
realized axion insulators are about 50 meV [36, 37, 16]. For example, if we
take $n_{2}=5$ and $d=1$ nm, we obtain a very small value $|\phi|\leq
1.27\times 10^{-3}$ at $\hbar\omega\leq 50$ meV, which is typically smaller
than $|\delta n|$. On the other hand, in the infrared and visible regime where
the photon energy is in the order of eV, equation (II.4) can become relevant.
We demonstrate this in the following section.

###  II.5 Model calculations

Figure 4:  Optical magneto-electric coupling in the low-energy tight-binding
model of MnBi2Te4. a and b, Real and imaginary part of the orbital magneto-
electric axion angle. Bulk LH indicates the optical layer Hall conductivity
calculated with periodic boundary conditions. Bulk LH approaches the exact
axion magneto-electric coupling as photon energy increases. c,d, Orbital
magneto-electric $T_{xx}$ and the spin magneto-electric coupling. Spin parts
contribute to axion angle and $T_{xx}$ through
$\theta^{(z),s}=(2G^{s}_{xx}+G^{s}_{zz})/3$ and
$T_{xx}^{s}=(G^{s}_{xx}-G^{s}_{zz})/3$. c, Spectra for two layers. These are
three orders of magnitudes smaller than the orbital axion magneto-electric
coupling. d, Layer-number dependence at photon energy $\hbar\omega=2$ eV.
Solid and dashed lines are curves fitted with $a/N_{l}+b$ form, where $a$ and
$b$ are fitting parameters. All optical response functions are calculated with
$\gamma=10$ meV to broaden the resonance through
$\omega\rightarrow\omega+i\gamma$.  Figure 5:  Kerr and Faraday rotation
angles in the low-energy tight-binding model of MnBi2Te4. a, Complex Kerr
rotation angle $\phi_{K}=\tan^{-1}(r_{xy}/r_{xx})$. b, Complex Faraday
rotation angle $\phi_{F}=\tan^{-1}(t_{xy}/t_{xx})$. The real and imaginary
parts of the complex angles are shown as solid and dashed lines, respectively.
We consider the effect of both top and bottom surfaces by using Eqs. (15) and
(16) with the refractive indices $n_{1}=2.2$ and $n_{3}=2.4$ for the
capsulating media. The complex refractive index of the model itself is
obtained from the electric dipole susceptibility through
$n_{2}(\omega)=\sqrt{1+\chi_{xx}(\omega)/\epsilon_{0}}$. For LH, we calculate
the axion angle from the layer Hall conductivity with periodic boundary
conditions and use the thickness of 50 layers for $\phi=n_{2}\omega d/c$. All
optical response functions are calculated with $\gamma=10$ meV to broaden the
resonance through $\omega\rightarrow\omega+i\gamma$.

Let us apply our theory to study the optical axion electrodynamics in
MnBi2Te4. MnBi2Te4 is the only stoichiometric compound that experimentally
realizes the intrinsic antiferromagnetic axion insulator [36, 37], which has
now become an attractive platform for studying axion magneto-electric effects
[38, 10, 11, 12, 39, 40, 41]. As it is a layered antiferromagnet, its few-
layer behavior and layer number dependence is also of interest [34, 42].

Here we calculate its magneto-optic properties based on the low-energy model
in Refs. [38, 43]. The goal of our calculations here is to cement the validity
of our new theory by providing a concrete model example as well as to
understand qualitative features (e.g., the dominance of the axion contribution
and the significant Kerr and negligible Faraday effects) of the magneto-optic
response in MnBi2Te4. Our model is expected to quantitatively capture the low-
energy properties of the material. On the other hand, at photon energies much
larger than the band gap, a precise quantitative calculation of the magneto-
optical spectrum will require a model, like the full ab-initio model, that
captures all the significant optical transitions involving those states
neglected in our model. With this in mind, we consider a nearest-neighbor
tight-binding Hamiltonian on a layer-stacked triangle lattice

$\displaystyle\hat{H}=\sum_{i,\alpha,\beta}\hat{c}^{\dagger}_{i\alpha}(h_{0})_{\alpha\beta}\hat{c}_{i\beta}-\sum_{\braket{i,j},\alpha,\beta}\hat{c}^{\dagger}_{i\alpha}t^{ij}_{\alpha\beta}\hat{c}_{j\beta},$
(20)

where $i,j$ are the site indices, $\braket{i,j}$ means that the summation is
over nearest neighbors, and $\alpha,\beta=1,\ldots,4$ run over two spin and
two orbital degrees of freedom at each site.

As the non-magnetic state has space group $R\bar{3}m$ (No. 166), we impose
time reversal $T=is_{y}K$, inversion $P=\tau_{z}$, and threefold
$C_{3z}=\exp\left(-i\pi s_{z}/3\right)$ and twofold $C_{2x}=-is_{x}$
rotational symmetries, where $s_{i}$ and $\tau_{i}$ are Pauli matrices for
spin and orbital, respectively. The onsite Hamiltonian satisfying all the
symmetries of the nonmagnetic state is $h_{0}=e_{0}+e_{5}\tau_{z}.$ Along the
$z$ direction, the nearest-neighbor hopping matrices are $T_{4}\equiv
t^{j+{\bf a}_{4},j}=t^{z}_{0}+it^{z}_{3}s_{z}\tau_{x}+t^{z}_{5}\tau_{z}$,
where ${\bf a}_{4}=(0,0,a_{z})$, $a_{z}$ is the out-of-plane lattice
parameter. For the in-plane directions, the hopping matrices are $t^{j+{\bf
a}_{1},j}=t_{0}+it_{1}s_{x}\tau_{x}+it_{4}\tau_{y}+t_{5}\tau_{z}=(t^{j,j+{\bf
a}_{1}})^{\dagger}$, $t^{j+{\bf a}_{2},j}=C_{3z}T_{1}C_{3z}^{-1}$, and
$t^{j+{\bf a}_{3},j}=C_{3z}T_{2}C_{3z}^{-1}$, where ${\bf a}_{1}=(a,0,0)$,
${\bf a}_{2}=C_{3z}{\bf a}_{1}$, ${\bf a}_{3}=C_{3z}{\bf a}_{2}$, $a$ is the
in-plane lattice parameter (see Appendix F for further details).

We consider the effect of the layer-alternating (i.e., A-type)
antiferromagnetism by adding $(-1)^{n-1}m\sigma_{z}$ to $h_{0}$, where $n$ is
the layer index. While this term breaks time reversal symmetry, inversion
symmetry remains the symmetry of the lattice system. However, finite even-
layer systems break inversion symmetry and have non-zero axion angle according
to Eq. (10).

Figure 4a shows the orbital part of the axion angle calculated with the tight-
binding parameters of MnBi2Te4 derived in Ref. [38] (We make a momentum-
dependent overall energy shift to obtain an insulating filling at half filling
as we describe in Appendix F. The band structure and electric susceptibility
are shown in Supplementary Figs. 1 and 2.). At high energies above 1 eV, the
axion angle is well approximated by the optical layer Hall conductivity for
any number of layers. However, the axion angle deviates significantly from the
optical layer Hall conductivity as the photon energy gets lower below 1 eV.
The deviation at the low energy increases with the number of layers because
the surface massive Dirac fermion is then more localized at the surfaces and
increases the second term in Eq. (10), making the static axion angle reach the
quantized value $\theta=\pi$.

$T_{xx}$ and spin magneto-electric coupling are much smaller than the axion
angle, as shown in Figs. 4b and 4c, respectively. This is consistent with our
expectation that these origin-independent magneto-electric couplings are
strongly suppressed in systems with local inversion symmetry. Furthermore,
they decrease inversely with the number of layers $N_{l}$ [34], because only
$O(1/N_{l})$ portion of layers near the top and bottom generates a nontrivial
response. This contrasts to the case with a finite $T_{xx}$ for
$N_{l}\rightarrow\infty$, where the response is coming from $O(N_{l})$ layers.

As $T_{xx}$ is relatively small, optical axion electrodynamics dominates the
Kerr and Faraday effects in this system. Figure 5 shows the Kerr and Faraday
rotation angles calculated with the magneto-electric coupling. We consider the
case where the model system (medium 2) is encapsulated by medium 1 and medium
3, having frequency-independent refractive indices $n_{1}=2.2$ and
$n_{3}=2.4$, respectively, corresponding to those of the hexagonal Boron
nitride and diamond at photon energy around $1$ eV [44, 45]. The calculated
Kerr rotation angle $\varphi_{K}$ is about $0.02^{\circ}$ at photon energies
larger than $1$ eV, which is about one order of magnitude smaller than
$\varphi_{K}\lesssim 1^{\circ}$ in typical ferromagnets although our
antiferromagnetic system has zero net magnetic moment. The Kerr angle in real
MnBi2Te4 can even be much enhanced because of the contributions from higher-
energy bands that we do not include here. As we consider $n_{1}\neq n_{3}$,
the Faraday rotation is nonzero because the cancellation between the top and
bottom surfaces is incomplete. The Faraday effect is two orders of magnitude
weaker than the Kerr effect.

##  III Discussion

Our theory of optical axion electrodynamics fills a crucial missing piece in
the macroscopic theory of magneto-optic effects in antiferromagnets developed
mostly by Graham and Raab [25, 26, 27, 21]. They used only origin-independent
$T_{ij}$ and $S_{ijk}$ to ensure physically meaningful results such as the
consistency with the reciprocal relations. However, our approach shows that we
can include the origin-dependent axion magneto-electric coupling in the
theory, and it is precisely the axion angle that controls the Kerr effect in
antiferromagnets with local inversion symmetry. As we show by using the low-
energy tight-binding model of MnBi2Te4, the omission of the axion
electrodynamics can underestimate the Kerr effect by orders of magnitudes. In
general, the same suppression of $T_{ij}$ and $S_{ijk}$ is expected in systems
with bulk symmetries that reverses an odd number of spacetime coordinates. As
long as those symmetries are broken at the surfaces, the axion magneto-
electric coupling is not much affected by the symmetries, such that axion
electrodynamics dominates the response.

In fact, the trace part of the magneto-electric coupling was included in the
study by Hornreich and Shtrikman [20] prior to Refs. [25, 26, 27, 21],
However, their estimate of the effect was four orders of magnitudes smaller
than the value experimentally observed subsequently [46]. This inconsistency
lead to the appearance of theories based on different approaches [24, 25, 26,
27], all of which do not include the axion electrodynamics. Our introduction
of the surface-sensitive magneto-electric coupling and the study of double
interfaces allow for precise quantitative understanding of the Kerr effect
including the axion electrodynamics, especially for thin films.

While we focus on the Kerr effect due to macroscopic magneto-electric coupling
for fully compensated antiferromagnets, there is another microscopic mechanism
based on the spatial modulation (phase change as well as the attenuation) of
the electric field $E(z)=E_{0}e^{ik_{z}z}$ [24], where $k_{z}=n\omega/c$ is
complex valued. However, the microscopic mechanism produces only minor
effects. To see this, let $\phi_{K,0}$ be the Kerr rotation angle by a single
layer with net magnetization in A-type antiferromagnet. Considering the
reflection of each layer as well as the complex phase rotation during the
propagation, one obtains the Kerr angle
$\phi_{K}=\phi_{K,0}(1-e^{2ik_{z}a_{z}}+e^{4ik_{z}a_{z}}-e^{6ik_{z}a_{z}}+\ldots)/(1+e^{2ik_{z}a_{z}}+e^{4ik_{z}a_{z}}+e^{6ik_{z}a_{z}}+\ldots)=-ik_{z}a_{z}\phi_{K,0}+O((k_{z}a_{z})^{2})$
for an even number layers [24], where $a_{z}$ is the layer spacing. This Kerr
angle is smaller than the axion-induced $\phi_{K}\approx\phi_{K,0}$ because
$a_{z}\ll\lambda$ for the wavelength down to the UV regime.

In the static limit, both Faraday and Kerr effects are often considered as
manifestations of the axion electrodynamics [8, 13, 14]. It is because the
systems under consideration have finite net Hall conductivity. The same sign
of the Hall conductivity is induced on the top and bottom surfaces of a
$Z_{2}$ topological insulator by either external magnetic fields [8, 13] or
coupling to ferromagnets [14]. The main focus of those studies is the
manifestation of the half-quantized surface Hall conductivity, rather than the
magneto-electric response of antiferromagnets.

An open question we leave for future studies is to formulate a quantum
geometric theory of optical axion electrodynamics in periodic systems,
generalizing the Chern-Simons integral in the static limit. A drawback of
calculating intra-cell optical magneto-electric coupling through Eq. (II.3)
(or calculating the layer Hall conductivity) is that it does not capture the
topological magneto-electric effect because we drop the second term in Eq.
(10). A unified formula that captures both intra-cell optical magneto-electric
coupling and topological magneto-electric effect is desired, and it is likely
to require extending the quantum geometric formulation for electric dipole
moments [47] to magnetic dipole and electric quadrupole moments and defining a
proper optical Chern-Simons integral. However, this may not be feasible, in
which case we are forced to work with quasi-two-dimensional systems.

Finally, we note that the optical axion angle we define should be
distinguished from the dymanical axion fields [48]. In our optical axion
electrodynamics based on linear response theory, the effective action $S_{\rm
OA}\propto\int d\omega d{\bf x}\theta(\omega,{\bf x}){\bf E}(\omega,{\bf
x})\cdot{\bf B}(\omega,{\bf x})$ for non-absorptive media describes the
propagation of the electromagnetic fields modified by elastic scattering by
the medium. The optical axion angle is a ground-state property, which is non-
dynamical. On the other hand, the dynamical axion field interacts spacetime-
locally with the electromagnetic fields through $S_{\rm dynamic}\propto\int
dtd{\bf x}\theta(t,{\bf x}){\bf E}(t,{\bf x})\cdot{\bf B}(t,{\bf x})$. This
describes Raman scattering where the energy or momentum of the incoming
electromagnetic field is tranfered to the dynamical axion field. The interplay
between the two distinct phenomena is an interesting research direction.

###### Acknowledgements.

We appreciate Jian-Xiang Qiu , Philip Kim, Ari Turner, Ivo Souza, and Allan
MacDonald for helpful discussions. This work was supported by the Center for
Advancement of Topological Semimetals, an Energy Frontier Research Center
funded by the U.S. Department of Energy Office of Science, Office of Basic
Energy Sciences, through the Ames Laboratory under contract No. DE-
AC02-07CH11358 (J.A., S.-Y.X., and A.V.).

##  Appendix A Generalization to include natural optical activity

###  A.1 Electromagnetic multipole moments

Let us consider electric dipole $P_{i}$, electric quadrupole $Q_{ij}$, and
magnetic dipole $M_{i}$ moment densities induced by electric and magnetic
fields [21]:

$\displaystyle P_{i}$
$\displaystyle=P^{0}_{i}+\chi_{ij}E_{j}+\frac{1}{2}a_{ijk}\nabla_{k}E_{j}+G_{ij}B_{j}$
$\displaystyle+\omega^{-1}\left[\chi^{\prime}_{ij}\dot{E}_{j}+\frac{1}{2}a^{\prime}_{ijk}\nabla_{k}\dot{E}_{j}+G^{\prime}_{ij}\dot{B}_{j}\right]\ldots$
$\displaystyle Q_{ij}$
$\displaystyle=Q^{0}_{ij}+{\mathfrak{a}}_{ijk}E_{k}+\omega^{-1}{\mathfrak{a}}^{\prime}_{ijk}\dot{E}_{k}\ldots$
$\displaystyle M_{i}$
$\displaystyle=M^{0}_{i}+{\mathfrak{G}}_{ij}E_{j}+\omega^{-1}{\mathfrak{G}}^{\prime}_{ij}\dot{E}_{j}\ldots$
(21)

for monochromatic electromagnetic fields in time domain, where $P^{0}_{i}$,
$Q^{0}_{ij}$, and $M^{0}_{i}$ are permanent multipole moments,
${\mathfrak{a}}_{ijk}=a_{kij}$,
${\mathfrak{a}}^{\prime}_{ijk}=-a^{\prime}_{kij}$,
${\mathfrak{G}}_{ij}=G_{ji}$, ${\mathfrak{G}}^{\prime}_{ij}=-G^{\prime}_{ji}$,
The ellipsis $``\ldots"$ indicates electric-octupole/magnetic-quadrupole or
higher-order multipole contributions that we neglect here. Here, the primed
susceptibility tensors transform oppositely under time reversal compared to
the non-primed ones. For example, while $\chi_{ij}$ is even under time
reversal, $\chi_{ij}^{\prime}$ is odd under time reversal.

Electromagnetic multipole moment densities are defined by [21]

$\displaystyle\hat{P}_{i}$ $\displaystyle=-e\hat{r}_{i}/V,$
$\displaystyle\hat{Q}_{ij}$ $\displaystyle=-e\hat{r}_{i}\hat{r}_{j}/V,$
$\displaystyle\hat{M}_{i}$
$\displaystyle=\frac{1}{4}\epsilon_{ijk}(\hat{r}_{j}\hat{J}^{\rm
orb}_{k}-\hat{J}^{\rm orb}_{j}\hat{r}_{k})+\hat{M}^{\rm spin}_{i}$
$\displaystyle=\hat{M}^{\rm orb}_{i}+\hat{M}^{\rm spin}_{i},$ (22)

where we split the orbital magnetic and spin parts, which respectively
originates from the minimal coupling $\nabla\rightarrow\nabla+ie{\bf A}$ and
the explicit dependence on ${\bf B}$ independent of the minimal coupling.

$\displaystyle\hat{J}^{\rm orb}_{i}$
$\displaystyle=\frac{1}{V}\frac{\partial\hat{H}}{\partial
A_{i}}{\bigg{|}}_{{\bf B}\text{ fixed}},$ $\displaystyle\hat{M}^{\rm
spin}_{i}$ $\displaystyle=-\frac{1}{V}\frac{\partial\hat{H}}{\partial
B_{i}}{\bigg{|}}_{{\bf A}\text{ fixed}}.$ (23)

###  A.2 Surface-sensitive magneto-electric coupling

By generalizing the procedure in the main text to include both natural optical
activity and gyrotropic birefringence, we define

$\displaystyle\tilde{G}^{(z)}_{ix}(\omega)$
$\displaystyle=\tilde{G}_{ix}(\omega)-\frac{i}{2}\omega\tilde{a}_{iyz}(\omega),$
$\displaystyle\tilde{G}^{(z)}_{iy}(\omega)$
$\displaystyle=\tilde{G}_{iy}(\omega)+\frac{i}{2}\omega\tilde{a}_{ixz}(\omega),$
$\displaystyle\tilde{G}^{(z)}_{zz}(\omega)$
$\displaystyle=\tilde{G}_{zz}(\omega)-\frac{i}{2}\omega\left[\tilde{\mathfrak{a}}_{zxy}(\omega)-\tilde{\mathfrak{a}}_{zyx}(\omega)\right]$
(24)

from the surface response

$\displaystyle j^{s}_{x}$
$\displaystyle=\sum_{i=1}^{3}\left(\tilde{\mathfrak{G}}_{yi}-\frac{i}{2}\omega\tilde{\mathfrak{a}}_{xzi}\right)E_{i}+\ldots,$
$\displaystyle j^{s}_{y}$
$\displaystyle=-\sum_{i=1}^{3}\left(\tilde{\mathfrak{G}}_{xi}+\frac{i}{2}\omega\tilde{\mathfrak{a}}_{yzi}\right)E_{i}+\ldots,$
$\displaystyle\rho^{s}$
$\displaystyle=\tilde{G}_{zi}B_{i}-\frac{i}{2}\omega(\tilde{\mathfrak{a}}_{zxy}-\tilde{\mathfrak{a}}_{zyx})B_{z}-\frac{i}{2}\omega\tilde{\mathfrak{a}}_{yzz}B_{x}+\frac{i}{2}\omega\tilde{\mathfrak{a}}_{xzz}B_{y}+\ldots,$
(25)

where $\tilde{G}_{ij}=G_{ij}-iG^{\prime}_{ij}$,
$\tilde{\mathfrak{G}}_{ij}=G_{ji}+iG^{\prime}_{ji}$,
$\tilde{a}_{ijk}=a_{ijk}-ia^{\prime}_{ijk}$, and
$\tilde{\mathfrak{a}}_{ijk}=a_{kij}+ia^{\prime}_{kij}$. We can also define

$\displaystyle\tilde{G}^{(z,2)}_{zi}(\omega)$
$\displaystyle=\tilde{G}_{zi}(\omega)-\frac{i}{2}\omega\sum_{k}\tilde{\mathfrak{a}}_{kzz}(\omega)\epsilon_{kzi}\text{
for }i=x,y,$ (26)

which is distinguished from $\tilde{G}^{(z)}_{zi}(\omega)$. However, we focus
on $\tilde{G}^{(z)}_{zi}(\omega)$ because we are interested in current
generation rather than the charge fluctuation on the surface. We decompose
$\tilde{G}^{(z)}_{ij}$ into

$\displaystyle\tilde{G}^{(z)}_{xx}(\omega)$ $\displaystyle=\frac{e^{2}}{2\pi
h}\theta^{(z)}(\omega)+\tilde{T}_{xx}(\omega)-\frac{\omega}{2}\left[S_{xyz}(\omega)+i\Gamma_{xyz}(\omega)\right],$
$\displaystyle\tilde{G}^{(z)}_{yy}(\omega)$ $\displaystyle=\frac{e^{2}}{2\pi
h}\theta^{(z)}(\omega)+\tilde{T}_{yy}(\omega)+\frac{\omega}{2}\left[S_{xyz}(\omega)+i\Gamma_{xyz}(\omega)\right],$
$\displaystyle\tilde{G}^{(z)}_{zz}(\omega)$ $\displaystyle=\frac{e^{2}}{2\pi
h}\theta^{(z)}(\omega)+\tilde{T}_{zz}(\omega),$
$\displaystyle\tilde{G}^{(z)}_{i\neq j}(\omega)$
$\displaystyle=\tilde{T}_{ij}(\omega)+\frac{\omega}{2}\epsilon_{zij}\left[S_{iiz}(\omega)+i\Gamma_{iiz}(\omega)\right],\text{
for }i,j\in 1,2,$ $\displaystyle\tilde{G}^{(z)}_{zi}(\omega)$
$\displaystyle=\tilde{T}_{zj}(\omega)+\frac{\omega}{2}\sum_{k=1}^{3}\epsilon_{kzi}\left[S_{zzk}(\omega)-ia_{kzz}\right],\text{
for }i\in 1,2,$ (27)

where

$\displaystyle\Gamma_{ijk}$ $\displaystyle=a_{ijk}+a_{jik}-a_{kij},$
$\displaystyle\tilde{T}_{ij}$
$\displaystyle=G_{ij}-\frac{1}{3}\delta_{ij}\sum_{k=1}^{3}G_{kk}-\frac{i}{6}\omega\sum_{k,l=1}^{3}\epsilon_{jkl}a^{\prime}_{kli}$
$\displaystyle\quad-i\left(G^{\prime}_{ij}-\sum_{k,l=1}^{3}\epsilon_{jkl}\frac{1}{2}\omega
a_{kli}\right),$ (28)

and $\theta^{(z)}(\omega)=2\pi he^{-2}\sum_{i=1}^{3}G^{(z)}_{ii}(\omega)/3$
and
$S_{ijk}(\omega)=\frac{1}{3}\left[a^{\prime}_{ijk}(\omega)+a^{\prime}_{jki}(\omega)+a^{\prime}_{kij}(\omega)\right]$
are the same as in the main text.

$\tilde{G}^{(z)}_{ij}(\omega)$ have nontrivial dependence under the change of
the spatial origin by ${\bf d}=(d_{x},d_{y},d_{z})$ because

$\displaystyle\delta\theta^{(z)}(\omega)$
$\displaystyle=-\sigma^{H}_{xy}(\omega)d_{z},$
$\displaystyle\delta\left(-\frac{i}{2}\omega\Gamma_{ijz}(\omega)\right)$
$\displaystyle=-\sigma^{\rm sym}_{ij}(\omega)d_{z},$
$\displaystyle\delta\left(-\frac{i}{2}\omega a_{kzz}(\omega)\right)$
$\displaystyle=-\sigma^{\rm sym}_{kz}(\omega)d_{z},$ (29)

where $\tilde{\sigma}_{ij}(\omega)=-i\omega\tilde{\chi}_{ij}(\omega)$ is the
complex-valued optical conductivity tensor, and
$\sigma^{H}_{ij}(\omega)=[\tilde{\sigma}_{ij}(\omega)-\tilde{\sigma}_{ji}(\omega)]/2=-\omega\chi^{\prime}_{ij}(\omega)$
is its anti-symmetric part.

The origin dependence shows the ambiguity of defining the surface degrees of
freedom. Let us recall that we define the surface current density
$j^{s}=\int^{d_{s}/2}_{-d_{s}/2}dzj^{i}(z)$ over the region $-d_{s}/2\leq
z\leq d_{s}/2$. While we can define a physically meaningful value of $d_{s}$
based on the surface property of a system, this value is still not completely
uniquely defined (e.g., if $d_{s}$ is the surface thickness, $1.01d_{s}$ also
makes sense as a surface thickness). Because of the ambiguity of $d_{s}$, the
amount of the bulk-conductivity contribution to $j^{s}$ can also vary. For
example, let us consider the $xy$ component of the surface conductivity. It
has a bulk-conductivity term as well as the magnetic dipole and electric
quadrupole contributions. Let us suppose that the surface is defined as the
interface between medium 1 (at $z<0$) and medium 2 at ($z>0$).

$\displaystyle\sigma^{s}_{xy}(\omega)$
$\displaystyle=\tilde{\mathfrak{G}}^{(z)}_{yy}(-d_{s}/2)-\tilde{\mathfrak{G}}^{(z)}_{yy}(d_{s}/2)+\int^{d_{s}/2}_{-d_{s}/2}dz\sigma_{xy}(\omega,z)$
$\displaystyle=\tilde{\mathfrak{G}}^{(z)}_{yy}(-d_{s}/2)-\tilde{\mathfrak{G}}^{(z)}_{yy}(d_{s}/2)$
$\displaystyle\quad+\sigma_{xy,1}(\omega)\frac{d_{s}}{2}+\sigma_{xy,2}(\omega)\frac{d_{s}}{2}$
$\displaystyle=\tilde{\mathfrak{G}}^{(z)}_{yy,1}+\delta\tilde{\mathfrak{G}}^{(z)}_{yy,1}-\tilde{\mathfrak{G}}^{(z)}_{yy,2}-\delta\tilde{\mathfrak{G}}^{(z)}_{yy,2},$
(30)

where

$\displaystyle\delta\tilde{\mathfrak{G}}^{(z)}_{yy}=-\sigma_{xy}(\omega)d_{z}$
(31)

is the change of $\tilde{\mathfrak{G}}^{(z)}_{yy}$ by the shifting of the
origin by ${\bf d}=(0,0,d_{z})$. This shows that, while we can define
$\sigma^{s}_{xy}(\omega)$ as the difference of
$\tilde{\mathfrak{G}}^{(z)}_{yy}$ between two media for any value of $d_{s}$,
we have to shift the spatial origin by $-d_{s}/2$ and $d_{s}/2$ for medium 1
and medium 2, respectively.

##  Appendix B Quantum mechanical expressions of the magneto-electric coupling

###  B.1 Linear response theory

The susceptibility tensor for the linear response of an operator $\hat{A}$ to
the external field $F_{B}$

$\displaystyle A(t)=\braket{\hat{A}(t)}_{F_{B}=0}+\int
dt^{\prime}\tilde{\chi}_{AB}(t-t^{\prime})F_{B}(t^{\prime})$ (32)

is given by

$\displaystyle\tilde{\chi}_{AB}(t-t^{\prime})$
$\displaystyle=-\frac{i}{\hbar}\Theta(t-t^{\prime})\braket{[\hat{A}(t),\hat{B}(t^{\prime})]},$
(33)

where $\hat{B}=\partial\hat{H}(F)/\partial F_{B}|_{F_{B}=0}$ is conjugate to
the external field. Here, we take the Heisenberg picture
$\hat{O}(t)=e^{i\hat{H}t}\hat{O}e^{-i\hat{H}t}$.

In the frequency domain, the susceptibility tensor can be written as

$\displaystyle\tilde{\chi}_{AB}(\omega)$
$\displaystyle=\int^{\infty}_{-\infty}dt^{\prime}e^{i\omega(t-t^{\prime})}\tilde{\chi}_{AB}(t-t^{\prime})$
$\displaystyle=-\frac{1}{\hbar}\sum_{n,m}f_{nm}\frac{\omega_{mn}}{\omega}\frac{\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}}{\omega_{mn}-\omega}$
(34)

plus FS terms, where $|n\rangle$ is the energy eigenstate of the unperturbed
Hamiltonian with energy $E_{n}=\hbar\omega_{n}$,
$\omega_{mn}=\omega_{m}-\omega_{n}$, and $f_{nm}=f_{n}-f_{m}$ is the
difference between the Fermi-Dirac distribution function $f$ of the
$|n\rangle$ and $|m\rangle$ states, and the FS terms originate from the Fermi
surface, i.e., they have momentum-space derivatives acting on the Fermi-Dirac
distribution function ($\partial_{\bf k}f_{n}$) in the momentum space
representation. The derivation goes as follows. At zero temperature, we have

$\displaystyle\tilde{\chi}_{AB}(\omega+i\Gamma)$
$\displaystyle=\int^{\infty}_{-\infty}dt^{\prime}e^{i(\omega+i\Gamma)(t-t^{\prime})}\tilde{\chi}_{AB}(t-t^{\prime})$
$\displaystyle=-\frac{i}{\hbar}\int^{\infty}_{-\infty}dt^{\prime}e^{i\omega(t-t^{\prime})}\Theta(t-t^{\prime})\braket{[\hat{A}(t),\hat{B}(t^{\prime})]}$
$\displaystyle=-\frac{i}{\hbar}\sum_{n,m}\int^{t}_{-\infty}dt^{\prime}e^{i(\omega+i\Gamma)(t-t^{\prime})}f_{n}(\braket{n}{\hat{A}(t)}{m}\braket{m}{\hat{B}(t^{\prime})}{n}-\braket{n}{\hat{B}(t^{\prime})}{m}\braket{m}{\hat{A}(t)}{n})$
$\displaystyle=-\frac{i}{\hbar}\sum_{n,m}\int^{t}_{-\infty}dt^{\prime}e^{i(\omega+i\Gamma)(t-t^{\prime})}f_{n}(e^{-i\omega_{mn}(t-t^{\prime})}\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}-e^{i\omega_{mn}(t-t^{\prime})}\braket{n}{\hat{B}}{m}\braket{m}{\hat{A}}{n})$
$\displaystyle=-\frac{i}{\hbar}\sum_{n,m}f_{n}\left(\frac{1}{-i(\omega+i\Gamma-\omega_{mn})}\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}-\frac{1}{-i(\omega+i\Gamma+\omega_{mn})}\braket{n}{\hat{B}}{m}\braket{m}{\hat{A}}{n}\right)$
$\displaystyle=-\frac{1}{\hbar}\sum_{n,m}\frac{f_{nm}\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}}{\omega_{mn}-(\omega+i\Gamma)}-\frac{1}{\hbar}\sum_{n,m}\left(f_{n}\frac{\braket{n}{\hat{B}}{m}\braket{m}{\hat{A}}{n}}{\omega+i\Gamma+\omega_{mn}}-f_{m}\frac{\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}}{\omega+i\Gamma+\omega_{nm}}\right)$
$\displaystyle=-\frac{1}{\hbar}\sum_{n,m}\frac{f_{nm}\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}}{\omega_{mn}-(\omega+i\Gamma)}-\frac{1}{\hbar}\frac{1}{\omega+i\Gamma}\sum_{n}f_{n}\braket{n}{[\hat{B}_{n},\hat{A}_{n}]}{n}$
$\displaystyle=-\frac{1}{\hbar}\sum_{n,m}f_{nm}\frac{\omega_{mn}}{\omega}\frac{\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}}{\omega_{mn}-(\omega+i\Gamma)}.$
(35)

where we introduce a finite relaxation rate $\Gamma$ for convergence of time
integral, $\hat{O}_{n}=\hat{P}_{n}\hat{O}\hat{P}_{n}$ with
$\hat{P}_{n}=|n\rangle\langle n|$ is the projection of $\hat{O}$ to states
$|n\rangle$, and
$\sum_{n}f_{n}\braket{n}{[\hat{B}_{n},\hat{A}_{n}]}{n}=\sum_{n\neq
m}f_{nm}A_{nm}B_{mn}$.

It is often convenient to separate the real and imaginary parts of
$\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}$ as follows.

$\displaystyle\tilde{\chi}_{AB}(\omega)$
$\displaystyle=-\frac{1}{2\hbar}\sum_{n,m}f_{nm}\frac{1}{\omega}\left[\frac{\omega_{mn}}{\omega_{mn}-\omega}-\frac{\omega_{nm}}{\omega_{nm}-\omega}\right]{\rm
Re}[\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}]-\frac{i}{2\hbar}\sum_{n,m}f_{nm}\frac{1}{\omega}\left[\frac{\omega_{mn}}{\omega_{mn}-\omega}+\frac{\omega_{nm}}{\omega_{nm}-\omega}\right]{\rm
Im}[\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}]$
$\displaystyle=-\frac{1}{\hbar}\sum_{n,m}f_{nm}\left(\frac{\omega_{mn}}{\omega_{mn}^{2}-\omega^{2}}{\rm
Re}[\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}]+i\frac{\omega_{mn}^{2}/\omega}{\omega_{mn}^{2}-\omega^{2}}{\rm
Im}[\braket{n}{\hat{A}}{m}\braket{m}{\hat{B}}{n}]\right)$
$\displaystyle=\chi_{AB}(\omega)-i\chi^{\prime}_{AB}(\omega),$ (36)

where we assume $\hat{A}$ and $\hat{B}$ are Hermitian operators.

Let us take electric susceptibility $\tilde{\chi}_{ij}$ as an example, which
is defined by

$\displaystyle P_{i}(t)=\braket{\hat{P}_{i}}_{E=0}+\int
dt^{\prime}\tilde{\chi}_{ij}(t-t^{\prime})E_{j}(t^{\prime}).$ (37)

Then, $\hat{A}=\hat{P}_{i}$ is the electric polarization density and
$\hat{B}=-\hat{P}_{j}V$ is the polarization, so we have

$\displaystyle\tilde{\chi}_{ij}(\omega)$
$\displaystyle=\frac{V}{\hbar}\sum_{n,m}f_{nm}\frac{\omega_{mn}}{\omega(\omega_{mn}-\omega)}\braket{n}{\hat{P}_{i}}{m}\braket{m}{\hat{P}_{j}}{n}$
$\displaystyle=\frac{2V}{\hbar}\sum_{n,m}f_{n}\frac{\omega_{mn}}{\omega_{mn}^{2}-\omega^{2}}{\rm
Re}[\braket{n}{\hat{P}_{i}}{m}\braket{m}{\hat{P}_{j}}{n}]$
$\displaystyle\quad+i\frac{2V}{\hbar}\sum_{n,m}f_{n}\frac{\omega_{mn}^{2}/\omega}{\omega_{mn}^{2}-\omega^{2}}{\rm
Im}[\braket{n}{\hat{P}_{i}}{m}\braket{m}{\hat{P}_{j}}{n}]$
$\displaystyle=\chi_{ij}(\omega)-i\chi^{\prime}_{ij}(\omega)$ (38)

Similarly, other susceptibility tensors are given by

$\displaystyle\tilde{a}_{ijk}(\omega)$
$\displaystyle=\frac{V}{\hbar}\sum_{n,m}f_{nm}\frac{\omega_{mn}}{\omega(\omega_{mn}-\omega)}\braket{n}{\hat{P}_{i}}{m}\braket{m}{\hat{Q}_{jk}}{n}=a_{ijk}-ia^{\prime}_{ijk},$
$\displaystyle\tilde{\mathfrak{a}}_{ijk}(\omega)$
$\displaystyle=\frac{V}{\hbar}\sum_{n,m}f_{nm}\frac{\omega_{mn}}{\omega(\omega_{mn}-\omega)}\braket{n}{\hat{Q}_{ij}}{m}\braket{m}{\hat{P}_{k}}{n}=a_{kij}+ia^{\prime}_{kij},$
$\displaystyle\tilde{G}_{ij}(\omega)$
$\displaystyle=\frac{V}{\hbar}\sum_{n,m}f_{nm}\frac{\omega_{mn}}{\omega(\omega_{mn}-\omega)}\braket{n}{\hat{P}_{i}}{m}\braket{m}{\hat{M}_{j}}{n}=G_{ij}-iG^{\prime}_{ij},$
$\displaystyle\tilde{\mathfrak{G}}_{ij}(\omega)$
$\displaystyle=\frac{V}{\hbar}\sum_{n,m}f_{nm}\frac{\omega_{mn}}{\omega(\omega_{mn}-\omega)}\braket{n}{\hat{M}_{i}}{m}\braket{m}{\hat{P}_{j}}{n}=G_{ji}+iG^{\prime}_{ji}.$
(39)

###  B.2 Structure of the optical magneto-electric coupling

The bare magneto-electric coupling can be decomposed into three parts as
follows.

$\displaystyle G_{ij}(\omega)$ $\displaystyle=G^{\rm K}_{ij}(\omega)+G^{\rm
C}_{ij}(\omega)+G^{A}_{ij}(\omega)$ $\displaystyle=\frac{e^{2}}{\hbar
V}\epsilon_{jkl}\sum_{n\in{\rm occ},m\in{\rm unocc},{\bf
k}}\frac{\omega_{mn}}{\omega^{2}_{mn}-\omega^{2}}{\rm
Re}\left[\sum_{n^{\prime}}r^{i}_{nm}r^{k}_{mn^{\prime}}v^{l}_{n^{\prime}n}+\sum_{m^{\prime}}r^{i}_{nm}v^{l}_{mm^{\prime}}r^{k}_{m^{\prime}n}\right]$
$\displaystyle+\frac{e^{2}}{\hbar V}\epsilon_{jkl}\sum_{n,n^{\prime}\in{\rm
occ},m\in{\rm unocc},{\bf
k}}\frac{\omega_{mn}^{2}}{\omega_{mn}^{2}-\omega^{2}}{\rm
Im}\left[r^{i}_{nm}r^{k}_{mn^{\prime}}r^{l}_{n^{\prime}n}\right]$
$\displaystyle+\frac{e^{2}}{\hbar V}\epsilon_{jkl}\sum_{n\in{\rm occ},m\in{\rm
unocc},{\bf
k}}\frac{\omega^{2}\omega_{mn}}{(\omega_{mn}^{2}-\omega^{2})^{2}}{\rm
Re}\left[r^{i}_{nm}r^{k}_{mn}\right](v^{l}_{n}+v^{l}_{m}).$ (40)

The $K$-term and $C$-term correspond to the Kubo-like and Chern-Simons terms
in Ref. [1] of the static limit. The last $A$-term is nonzero only when
$\omega\neq 0$. While the $C$-term was called the Chern-Simons term, it has
additional terms, actually. In the static limit,

$\displaystyle G^{\rm C}_{ij}(0)$ $\displaystyle=\delta_{ij}\frac{\theta_{\rm
CS}e^{2}}{2\pi h}+\frac{1}{3}\epsilon_{jkl}\lim_{\omega\rightarrow
0}\left[\omega
a^{\prime}_{kli}(\omega)\right]+\frac{e^{2}}{6\hbar}\epsilon_{jkl}\int_{\bf
k}\partial_{k}g_{il},$ (41)

Here, $\theta_{\rm CS}$ is the axion angle given by the Chern-Simons integral
in the three-dimensional Brillouin zone [17, 18]

$\displaystyle\theta_{\rm CS}$
$\displaystyle=-\frac{4\pi^{2}}{3V}\epsilon_{ijk}{\rm Im}{\rm
Tr}\left[P\hat{r}^{i}P\hat{r}^{j}P\hat{r}^{k}\right]$
$\displaystyle=-\frac{1}{4\pi}\epsilon_{ijk}\int_{\rm
BZ}d^{3}k\bigg{[}\left(A^{i}\frac{1}{2}F^{jk}+\frac{i}{3}A^{i}A^{j}A^{k}\right)-\sum_{n,n^{\prime}\in{\rm
occ}}\partial_{j}{\rm
Re}\left(\braket{i\partial_{i}\psi_{n}}{\psi_{n^{\prime}}}A^{k}_{n^{\prime}n}\right)\bigg{]},$
(42)

where $P$ is the projection to the occupied states, and
$A^{k}_{n^{\prime}n}=\braket{u_{n^{\prime}}}{i\partial_{k}}{u_{n}}$ is the
non-abelian Berry connection for the occupied states. Since the last term
vanishes in insulators with vanishing Chern number, only the Chern-Simons
integral remains in the expression. The quadrupole term takes the form
$\lim_{\omega\rightarrow 0}\left[\omega
a^{\prime}_{kli}(\omega)\right]=-\frac{e^{2}}{\hbar V}\epsilon_{jkl}{\rm
Im}{\rm Tr}\left[P\hat{r}^{k}Q\hat{r}^{l}\hat{r}^{i}\right]$. Lastly,
$g_{ij}=\sum_{n\in{\rm occ},m\in{\rm unocc}}r^{i}_{nm}r^{j}_{mn}$ is the
quantum metric of the occupied states. The quantum metric term in Eq. (41)
cancels the Fermi surface contribution from the quadrupole term (the quantum
metric contribution in electric quadrupole responses was discussed in [49,
50]), such that there is no Fermi surface contribution to the magneto-electric
coupling $G_{ij}$. In comparison, note that its time-reversal-symmetric
counterpart $G^{\prime}_{ij}$ has a Fermi surface contribution, which gives
rise to natural optical activity in metals, termed gyrotropic magnetic effect
[2, 52]. The quadrupole and quantum metric terms were missed in previous
studies [17, 18], but they do not affect the axion angle.

Let us derive Eq. (B.2). A key equation in our derivation is

$\displaystyle\sum_{j,k}\epsilon_{ijk}\braket{i\partial_{j}\psi_{m}}{\psi_{p}}\braket{\psi_{p}}{i\partial_{k}\psi_{n}}=\sum_{j,k}\epsilon_{ijk}\braket{\partial_{j}\psi_{m}}{\partial_{k}\psi_{n}}=0,$
(43)

which follows from the Wannier representation:

$\displaystyle\sum_{j,k}\epsilon_{ijk}\braket{\partial_{j}\psi_{m}}{\partial_{k}\psi_{n}}$
$\displaystyle=\epsilon_{ijk}\frac{1}{N}\sum_{j,k,\bf
R,R^{\prime}}R_{j}R^{\prime}_{k}e^{i({\bf k}\cdot{\bf R}-{\bf
k}^{\prime}\cdot{\bf R}^{\prime})}\braket{w_{m{\bf R}^{\prime}}}{w_{n\bf R}}$
$\displaystyle=\delta_{mn}\frac{1}{N}\sum_{\bf R}e^{i{\bf k}\cdot({\bf R}-{\bf
R}^{\prime})}\sum_{j,k}\epsilon_{ijk}R_{j}R_{k}$ $\displaystyle=0.$ (44)

Using this identity, we can show that

$\displaystyle\sum_{i,j,k}\epsilon_{ijk}{\rm
TrIm}\left[P\hat{r}^{i}P\hat{r}^{j}P\hat{r}^{k}\right]$
$\displaystyle=-\sum_{i,j,k}\epsilon_{ijk}{\rm
TrIm}\left[P\hat{r}^{i}P\hat{r}^{j}Q\hat{r}^{k}\right]$
$\displaystyle=\sum_{i,j,k}\epsilon_{ijk}{\rm
Re}\left[\left(A^{i}_{nn^{\prime}}-\braket{i\partial_{i}\psi_{n}}{\psi_{n^{\prime}}}\right)\frac{1}{2}F^{jk}_{n^{\prime}n}\right]$
$\displaystyle=\sum_{i,j,k}\epsilon_{ijk}{\rm
Re}\left[A^{i}_{nn^{\prime}}\frac{1}{2}F^{jk}_{n^{\prime}n}-\braket{i\partial_{i}\psi_{n}}{\psi_{n^{\prime}}}\left(\partial_{j}A^{k}_{n^{\prime}n}-i(A^{j}A^{k})_{n^{\prime}n}\right)\right]$
$\displaystyle=\sum_{i,j,k}\epsilon_{ijk}\left[A^{i}_{nn^{\prime}}\frac{1}{2}F^{jk}_{n^{\prime}n}+i\braket{i\partial_{i}\psi_{n}}{\psi_{n^{\prime}}}(A^{j}A^{k})_{n^{\prime}n}-\partial_{j}{\rm
Re}\left(\braket{i\partial_{i}\psi_{n}}{\psi_{n^{\prime}}}A^{k}_{n^{\prime}n}\right)\right]$
$\displaystyle=\sum_{i,j,k}\epsilon_{ijk}\left[A^{i}_{nn^{\prime}}\frac{1}{2}F^{jk}_{n^{\prime}n}-i\frac{1}{3}{\rm
Tr}\left[P\hat{r}^{i}P\hat{r}^{j}P\hat{r}^{k}-A^{i}A^{j}A^{k}\right]-\partial_{j}{\rm
Re}\left(\braket{i\partial_{i}\psi_{n}}{\psi_{n^{\prime}}}A^{k}_{n^{\prime}n}\right)\right]$
$\displaystyle=\frac{1}{3}\sum_{i,j,k}\epsilon_{ijk}{\rm
TrIm}\left[P\hat{r}^{i}P\hat{r}^{j}P\hat{r}^{k}\right]+\sum_{i,j,k}\epsilon_{ijk}\left[{\rm
Tr}\left(A^{i}\frac{1}{2}F^{jk}+\frac{i}{3}A^{i}A^{j}A^{k}\right)-\partial_{j}{\rm
Re}\left(\braket{i\partial_{i}\psi_{n}}{\psi_{n^{\prime}}}A^{k}_{n^{\prime}n}\right)\right].$
(45)

It follows that

$\displaystyle\frac{1}{3}{\rm Tr}G$
$\displaystyle=\frac{e^{2}}{\hbar}\frac{1}{3}\epsilon_{ijk}{\rm
TrIm}\left[P\hat{r}^{i}P\hat{r}^{j}P\hat{r}^{k}\right]$
$\displaystyle=\frac{e^{2}}{\hbar}\frac{1}{2}\epsilon_{ijk}\bigg{[}{\rm
Tr}\left(A^{i}\frac{1}{2}F^{jk}+\frac{i}{3}A^{i}A^{j}A^{k}\right)-\partial_{j}{\rm
Re}\left(\braket{i\partial_{i}\psi_{n}}{\psi_{n^{\prime}}}A^{k}_{n^{\prime}n}\right)\bigg{]}.$
(46)

This is equivalent to Eq. (B.2). At finite $\omega$, the trace part of the
magneto-electric coupling is not represented as a Chern-Simons integral by the
same approach.

###  B.3 Gauge-invariant expressions

Since $\tilde{G}^{(z)}_{ij}(\omega)$ and $\tilde{G}^{(z,2)}_{ij}(\omega)$ are
origin dependent only along the $z$ direction, they can be calculated gauge
independently when the $z$ direction has open boundaries. Here we derive such
gauge-invariant expressions. We focus on the orbital magneto-electric coupling
to simplify expressions. It is straightforward to include spin parts as the
spin magnetic moment is well defined in momentum space. The expressions of
$T_{ij}$, $T_{ij}^{\prime}$ and $S_{ijk}$ are also given for completeness.

Let us begin with $\tilde{G}^{(z)}_{ix}(\omega)$.

$\displaystyle\tilde{G}^{(z)}_{ix}(\omega)$
$\displaystyle=\tilde{G}_{ix}(\omega)-\frac{i}{2}\omega\tilde{a}_{iyz}(\omega)$
$\displaystyle=\frac{V}{\hbar}\sum_{n\neq
m}f_{nm}\frac{\braket{n}{\hat{P}_{i}}{m}\braket{m}{\hat{m}_{x}}{n}-\frac{i}{2}\omega_{mn}\braket{n}{\hat{P}_{i}}{m}\braket{m}{\hat{Q}_{yz}}{n}}{\omega_{mn}-\omega}$
$\displaystyle=\frac{e^{2}}{\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}\left(r^{i}_{nm}\braket{m}{\frac{1}{2}\left(\hat{r}^{y}\hat{v}^{z}-\hat{r}^{z}\hat{v}^{y}\right)}{n}-\frac{i}{2}\omega_{mn}r^{i}_{nm}\braket{m}{\hat{r}^{y}\hat{r}^{z}}{n}\right)$
$\displaystyle=\frac{e^{2}}{\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}r^{i}_{nm}\braket{m}{\frac{1}{2}\left(\hat{r}^{y}\hat{v}^{z}-\hat{r}^{z}\hat{v}^{y}-\frac{1}{i\hbar}[\hat{r}^{y}\hat{r}^{z},\hat{H}]\right)}{n}$
$\displaystyle=\frac{e^{2}}{\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}r^{i}_{nm}\braket{m}{\frac{1}{2}\left(-\hat{v}^{y}\hat{r}^{z}-\hat{r}^{z}\hat{v}^{y}\right)}{n}.$
(47)

Similarly,

$\displaystyle\tilde{G}^{(z)}_{iy}(\omega)=\tilde{G}_{iy}(\omega)+\frac{i}{2}\omega\tilde{a}_{ixz}(\omega)$
$\displaystyle=\frac{e^{2}}{\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}r^{i}_{nm}\braket{m}{\frac{1}{2}\left(\hat{v}^{x}\hat{r}^{z}+\hat{r}^{z}\hat{v}^{x}\right)}{n}.$
(48)

For the $zz$ component, we obtain

$\displaystyle\tilde{G}^{(z)}_{zz}(\omega)$
$\displaystyle=\tilde{G}_{zz}-\frac{i}{2}\omega\left[\tilde{\mathfrak{a}}_{zxy}(\omega)-\tilde{\mathfrak{a}}_{zyx}(\omega)\right]$
$\displaystyle=\frac{e^{2}}{2\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}\sum_{p;E_{p}\neq
E_{m}}\left[r^{z}_{nm}r^{x}_{mp}v^{y}_{pn}-r^{z}_{np}r^{x}_{pm}v^{y}_{mn}-(x\leftrightarrow
y)\right]$ (49)

by using

$\displaystyle\tilde{G}_{zz}$ $\displaystyle=\frac{V}{\hbar}\sum_{n\neq
m}f_{nm}\frac{\braket{n}{\hat{P}_{z}}{m}\braket{m}{\hat{M}_{z}}{n}}{\omega_{mn}-\omega}$
$\displaystyle=\frac{e^{2}}{2\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}r^{z}_{nm}\left[\sum_{p\neq
m}(r^{x}_{mp}v^{y}_{pn}-r^{y}_{mp}v^{x}_{pn})+i\omega_{mn}(r^{x}_{mm}r^{y}_{mn}-r^{y}_{mm}r^{x}_{mn})\right]$
$\displaystyle=\frac{e^{2}}{2\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}\left[r^{z}_{nm}\sum_{p\neq
m}(r^{x}_{mp}v^{y}_{pn}-r^{y}_{mp}v^{x}_{pn})+i\omega_{mn}\left(\left[(\hat{r}^{z}\hat{r}^{x})_{nm}-\sum_{p\neq
m}r^{z}_{np}r^{x}_{pm}\right]r^{y}_{mn}-(x\rightarrow y)\right)\right]$
$\displaystyle=\frac{e^{2}}{2\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}\left[r^{z}_{nm}\sum_{p;E_{p}\neq
E_{m}}(r^{x}_{mp}v^{y}_{pn}-r^{y}_{mp}v^{x}_{pn})-i\omega_{mn}\sum_{p;E_{p}\neq
E_{m}}\left[r^{z}_{np}r^{x}_{pm}-(x\leftrightarrow y)\right]r^{y}_{mn}\right]$
$\displaystyle\qquad+\frac{i}{2}\omega\left[\tilde{\mathfrak{a}}_{zxy}(\omega)-\tilde{\mathfrak{a}}_{zyx}(\omega)\right].$
(50)

For the $zx$ component, we can follow the strategy for the $zz$ component

$\displaystyle\tilde{G}_{zx}$ $\displaystyle=\frac{V}{\hbar}\sum_{n\neq
m}f_{nm}\frac{\braket{n}{\hat{P}_{z}}{m}\braket{m}{\hat{M}_{x}}{n}}{\omega_{mn}-\omega}$
$\displaystyle=\frac{e^{2}}{2\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}r^{z}_{nm}\left[\sum_{p\neq
m}r^{y}_{mp}v^{z}_{pn}-\sum_{p}r^{z}_{mp}v^{y}_{pn}+i\omega_{mn}r^{y}_{mm}r^{z}_{mn}\right]$
$\displaystyle=\frac{e^{2}}{2\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}\left[r^{z}_{nm}\sum_{p\neq
m}r^{y}_{mp}v^{z}_{pn}-r^{z}_{nm}\sum_{p}r^{z}_{mp}v^{y}_{pn}+i\omega_{mn}\left(\left[(\hat{r}^{z}\hat{r}^{y})_{nm}-\sum_{p\neq
m}r^{z}_{np}r^{y}_{pm}\right]r^{z}_{mn}\right)\right]$
$\displaystyle=\frac{e^{2}}{2\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}\left[r^{z}_{nm}\sum_{p;E_{p}\neq
E_{m}}r^{y}_{mp}v^{z}_{pn}-r^{z}_{nm}\sum_{p}r^{z}_{mp}v^{y}_{pn}-i\omega_{mn}\sum_{p;E_{p}\neq
E_{m}}r^{z}_{np}r^{y}_{pm}r^{z}_{mn}\right]$
$\displaystyle\qquad+\frac{i}{2}\omega\tilde{\mathfrak{a}}_{zyz}(\omega)$ (51)

to define

$\displaystyle\tilde{G}^{(z,2)}_{zx}$
$\displaystyle=\tilde{G}_{zx}-\frac{i}{2}\omega\tilde{\mathfrak{a}}_{yzz}(\omega)$
$\displaystyle=\frac{e^{2}}{2\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}\left[r^{z}_{nm}\left(\sum_{p;E_{p}\neq
E_{m}}r^{y}_{mp}v^{z}_{pn}-\sum_{p}r^{z}_{mp}v^{y}_{pn}\right)-\sum_{p;E_{p}\neq
E_{m}}r^{z}_{np}r^{y}_{pm}v^{z}_{mn}\right].$ (52)

Similarly, we define

$\displaystyle\tilde{G}^{(z,2)}_{zy}$
$\displaystyle=\tilde{G}_{zy}+\frac{i}{2}\omega\tilde{\mathfrak{a}}_{xzz}(\omega)$
$\displaystyle=-\frac{e^{2}}{2\hbar V}\sum_{n\neq
m}\frac{f_{nm}}{\omega_{mn}-\omega}\left[r^{z}_{nm}\left(\sum_{p;E_{p}\neq
E_{m}}r^{x}_{mp}v^{z}_{pn}-\sum_{p}r^{z}_{mp}v^{x}_{pn}\right)-\sum_{p;E_{p}\neq
E_{m}}r^{z}_{np}r^{x}_{pm}v^{z}_{mn}\right].$ (53)

Note that $G^{(z)}_{zi}=G^{(z,2)}_{zi}$ while $G^{\prime(z)}_{zi}\neq
G^{\prime(z,2)}_{zi}$ for $i=x,y$ .

Fully origin-independent bulk response functions $T_{ij}$, $T_{ij}^{\prime}$
and $S_{ijk}$ can be calculated from the bulk conductivity tensor defined by
$J_{i}=\sum_{j,k}\sigma_{ijk}q_{j}E_{k}$ [1]:

$\displaystyle T_{ij}(\omega)$
$\displaystyle=\frac{1}{3i}\epsilon_{jkl}\sigma_{ikl},$ $\displaystyle
T^{\prime}_{ij}(\omega)$
$\displaystyle=\frac{1}{8i}\epsilon_{jkl}\left[2(\sigma_{ikl}-\sigma_{kil})-(\sigma_{kli}-\sigma_{lki})\right],$
$\displaystyle S_{ijk}(\omega)$
$\displaystyle=-\frac{1}{6i}\left(\sigma_{ijk}+\sigma_{jki}+\sigma_{kij}+\sigma_{jik}+\sigma_{ikj}+\sigma_{kji}\right),$
(54)

where

$\displaystyle\sigma_{ikl}$
$\displaystyle=\frac{ie^{2}}{\hbar\omega}\sum_{n,m}\int_{\bf
k}f_{nm}\bigg{[}\frac{\omega_{mn}}{\omega_{mn}-\omega}\left((r^{j}_{nm}B^{ik}_{mn})^{*}+r^{i}_{nm}B^{jk}_{mn}\right)$
$\displaystyle\qquad\quad+\frac{\omega_{mn}^{2}}{2(\omega_{mn}-\omega)^{2}}r^{i}_{nm}r^{j}_{mn}(v^{k}_{mm}+v^{k}_{nn})\bigg{]}$
$\displaystyle\qquad-\frac{ie^{2}}{\hbar\omega^{2}}\sum_{n,m}\int_{\bf
k}\partial_{k}f_{n}v^{i}_{nn}v^{j}_{nn}$ (55)

and

$\displaystyle B^{kl}_{mn}=\frac{1}{2}\left(\sum_{p\neq
m}r^{k}_{mp}v^{l}_{pn}+\sum_{p\neq n}v^{l}_{mp}r^{k}_{pn}\right).$ (56)

##  Appendix C Decomposition of the position operator

Here we derive Eq. (II.3). Let us consider the matrix element of the position
operator in the Bloch state basis. We transform it to the Wannier basis
$|w_{\alpha{\bf R}}\rangle$ by

$\displaystyle\braket{\psi_{m\bf k^{\prime}}}{\hat{r}^{z}}{\psi_{n\bf k}}$
$\displaystyle=\sum_{\beta,{\bf R}^{\prime};\alpha,{\bf R}}\braket{\psi_{m\bf
k^{\prime}}}{w_{\beta{\bf R^{\prime}}}}\braket{w_{\beta{\bf
R}^{\prime}}}{\hat{r}^{z}}{w_{\alpha{\bf R}}}\braket{w_{\alpha{\bf
R}}}{\psi_{n\bf k}}$ $\displaystyle=\sum_{\beta,{\bf R}^{\prime};\alpha,{\bf
R}}\braket{\psi_{m\bf k^{\prime}}}{w_{\beta{\bf
R^{\prime}}}}\braket{w_{\beta{\bf 0}}}{\hat{T}_{\bf
R^{\prime}}^{\dagger}\hat{r}^{z}\hat{T}_{\bf R^{\prime}}}{w_{\alpha{\bf
R-R^{\prime}}}}\braket{w_{\alpha{\bf R}}}{\psi_{n\bf k}}$
$\displaystyle=\sum_{\beta,{\bf R}^{\prime};\alpha,{\bf R}}\braket{\psi_{m\bf
k^{\prime}}}{w_{\beta{\bf R^{\prime}}}}\braket{w_{\beta{\bf
0}}}{\hat{r}^{z}}{w_{\alpha({\bf R}-{\bf R}^{\prime})}}\braket{w_{\alpha{\bf
R}}}{\psi_{n\bf k}}$ $\displaystyle+\quad\sum_{\beta,{\bf
R}^{\prime};\alpha,{\bf R}}\braket{\psi_{m\bf k^{\prime}}}{w_{\beta{\bf
R^{\prime}}}}R^{z}\delta_{\alpha\beta}\delta_{{\bf R},{\bf
R}^{\prime}}\braket{w_{\alpha{\bf R}}}{\psi_{n\bf k}}$
$\displaystyle=\sum_{\beta,{\bf R}^{\prime};\alpha,{\bf R}}\braket{\psi_{m\bf
k^{\prime}}}{w_{\beta{\bf R^{\prime}}}}\braket{w_{\beta{\bf
0}}}{\hat{r}^{z}}{w_{\alpha({\bf R}-{\bf R}^{\prime})}}\braket{w_{\alpha{\bf
R}}}{\psi_{n\bf k}}$ $\displaystyle\quad+\sum_{\alpha,{\bf
R}}\braket{\psi_{m\bf k^{\prime}}}{w_{\alpha{\bf
R}}}R^{z}\braket{w_{\alpha{\bf R}}}{\psi_{n\bf k}},$ (57)

where we use $|w_{\alpha{\bf R}}\rangle=\hat{T}_{\bf R}|w_{\alpha{\bf
0}}\rangle$ in the second line, where $\hat{T}_{\bf R}$ is the translation
operator by ${\bf R}$, and use $\hat{T}_{\bf
R^{\prime}}^{\dagger}\hat{r}^{z}\hat{T}_{\bf R^{\prime}}=\hat{r}^{z}+{\bf
R^{\prime}}^{z}$ and the orthonormality of the Wannier states
$\braket{w_{\beta{\bf R}^{\prime}}}{w_{\alpha{\bf
R}}}=\delta_{\alpha\beta}\delta_{{\bf R},{\bf R}^{\prime}}$ in the third line.
To get Eq. (II.3) from Eq. (C), we define ${\mathbb{A}}^{z}_{\beta\alpha}({\bf
k})=\sum_{\bf R}\braket{w_{\beta{\bf 0}}}{\hat{r}^{z}}{w_{\alpha{\bf
R}}}e^{i{\bf k}\cdot{\bf R}}$, such that

$\displaystyle\braket{w_{\beta{\bf 0}}}{\hat{r}^{z}}{w_{\alpha{\bf
R}}}=\frac{1}{N}\sum_{\bf k}{\mathbb{A}}^{z}_{\beta\alpha}({\bf k})e^{-i{\bf
k}\cdot{\bf R}},$ (58)

where $N$ is the number of ${\bf k}$ points. Then, the first term in the last
line of Eq. (C) becomes

$\displaystyle\sum_{\beta,{\bf R}^{\prime};\alpha,{\bf R}}\braket{\psi_{m\bf
k^{\prime}}}{w_{\beta{\bf R^{\prime}}}}\braket{w_{\beta{\bf
0}}}{\hat{r}^{z}}{w_{\alpha({\bf R}-{\bf R}^{\prime})}}\braket{w_{\alpha{\bf
R}}}{\psi_{n\bf k}}$ $\displaystyle=\sum_{\beta,{\bf R}^{\prime};\alpha,{\bf
R}}\braket{\psi_{m\bf k^{\prime}}}{w_{\beta{\bf 0}}}\frac{1}{N}\sum_{\bf
k^{\prime\prime}}{\mathbb{A}}^{z}_{\beta\alpha}({\bf
k^{\prime\prime}})e^{-i{\bf k^{\prime\prime}}\cdot{\bf(R-R^{\prime})}}$
$\displaystyle\quad\times e^{i({\bf k\cdot R-k^{\prime}\cdot
R^{\prime}})}\braket{w_{\alpha{\bf 0}}}{\psi_{n\bf k}}$
$\displaystyle=\sum_{\alpha,\beta,{\bf R}^{\prime}}\braket{\psi_{m\bf
k^{\prime}}}{w_{\beta{\bf 0}}}{\mathbb{A}}^{z}_{\beta\alpha}({\bf k})e^{i{\bf
k}\cdot{\bf R^{\prime}}}e^{-i{\bf k^{\prime}\cdot
R^{\prime}}}\braket{w_{\alpha{\bf 0}}}{\psi_{n\bf k}}$
$\displaystyle=\delta_{\bf k,k^{\prime}}N\sum_{\alpha,\beta}\braket{\psi_{m\bf
k}}{w_{\beta{\bf 0}}}{\mathbb{A}}^{z}_{\beta\alpha}({\bf
k})\braket{w_{\alpha{\bf 0}}}{\psi_{n\bf k}}$ $\displaystyle=\delta_{\bf
k,k^{\prime}}\sum_{\alpha,\beta}\braket{\psi_{m\bf k}}{\psi_{\beta{\bf
k}}}{\mathbb{A}}^{z}_{\beta\alpha}({\bf k})\braket{\psi_{\alpha{\bf
k}}}{\psi_{n\bf k}},$ (59)

and we arrive at Eq. (II.3).

$\displaystyle\braket{\psi_{m\bf k^{\prime}}}{\hat{r}^{z}}{\psi_{n\bf k}}$
$\displaystyle=\delta_{\bf k,k^{\prime}}\sum_{\beta,\alpha}\braket{\psi_{m\bf
k}}{\psi_{\beta{\bf k}}}{\mathbb{A}}^{z}_{\beta\alpha}({\bf
k})\braket{\psi_{\alpha{\bf k}}}{\psi_{n\bf k}}$
$\displaystyle\quad+\sum_{\alpha,{\bf R}}\braket{\psi_{m\bf
k^{\prime}}}{w_{\alpha{\bf R}}}R^{z}\braket{w_{\alpha{\bf R}}}{\psi_{n\bf
k}}.$ (60)

The second term is nonzero only when $n=m$ because

$\displaystyle\sum_{\alpha,{\bf R}}\braket{\psi_{m\bf
k^{\prime}}}{w_{\alpha{\bf R}}}R^{z}\braket{w_{\alpha{\bf R}}}{\psi_{n\bf k}}$
$\displaystyle=\frac{1}{N}\sum_{\alpha,{\bf R},{\bf k}_{1},{\bf
k}_{2}}e^{i{({\bf k}_{2}-{\bf k}_{1})\cdot{\bf R}}}\braket{\psi_{m\bf
k^{\prime}}}{\psi_{\alpha{\bf k}_{1}}}R^{z}\braket{\psi_{\alpha{\bf
k}_{2}}}{\psi_{n\bf k}}$ $\displaystyle=\frac{1}{N}\sum_{{\bf R}}e^{i{({\bf
k}-{\bf k}^{\prime})\cdot{\bf R}}}R^{z}\sum_{\alpha}\braket{\psi_{m\bf
k^{\prime}}}{\psi_{\alpha{\bf k}^{\prime}}}\braket{\psi_{\alpha{\bf
k}}}{\psi_{n\bf k}}$ $\displaystyle=-i\partial_{k^{z}}\delta_{\bf
k,k^{\prime}}\sum_{\alpha}\braket{\psi_{m\bf k^{\prime}}}{\psi_{\alpha{\bf
k}^{\prime}}}\braket{\psi_{\alpha{\bf k}}}{\psi_{n\bf k}}$
$\displaystyle=-i\delta_{mn}\partial_{k^{z}}\delta_{\bf k,k^{\prime}},$ (61)

where we use that $\partial_{k^{z}}\delta_{\bf k,k^{\prime}}\neq 0$ requires
${\bf k}^{\prime}\rightarrow{\bf k}$, and $\sum_{\alpha}\braket{\psi_{m\bf
k}}{\psi_{\alpha{\bf k}}}\braket{\psi_{\alpha{\bf k}}}{\psi_{n\bf
k}}=\delta_{mn}$.

##  Appendix D Macroscopic electrodynamics in the medium

Electric displacement ${\bf D}$ and magnetic field ${\bf H}$ satisfying
Maxwell’s equations

$\displaystyle\nabla\cdot{\bf D}$ $\displaystyle=\rho^{f},$
$\displaystyle\nabla\times{\bf H}$ $\displaystyle={\bf J}^{f}+\dot{\bf D}$
(62)

are defined by

$\displaystyle\begin{pmatrix}{\bf D}\\\ {\bf
H}\end{pmatrix}=\begin{pmatrix}\tilde{A}&\tilde{T}\\\
\tilde{U}&\tilde{X}\end{pmatrix}\begin{pmatrix}{\bf E}\\\ {\bf
B}\end{pmatrix},$ (63)

where $\tilde{F}=F-iF^{\prime}$ for $F=A,T,U,X$, and

$\displaystyle A_{ij}$
$\displaystyle=\epsilon_{0}\delta_{ij}+\chi_{ij}+\frac{1}{3}(a^{\prime}_{ijk}+a^{\prime}_{jki}+a^{\prime}_{kij})k_{k},$
$\displaystyle A^{\prime}_{ij}$ $\displaystyle=\chi_{ij}^{\prime},$
$\displaystyle T_{ij}$
$\displaystyle=G_{ij}-\frac{1}{3}\delta_{ij}G_{kk}-\frac{1}{6}\omega\epsilon_{jkl}a^{\prime}_{kli},$
$\displaystyle T^{\prime}_{ij}$
$\displaystyle=G^{\prime}_{ij}-\frac{1}{2}\omega\epsilon_{jkl}a_{kli},$
$\displaystyle U_{ij}$ $\displaystyle=-U_{ji},$ $\displaystyle
U^{\prime}_{ij}$ $\displaystyle=U_{ji},$ $\displaystyle X_{ij}$
$\displaystyle=\mu_{0}^{-1}\delta_{ij},$ $\displaystyle X^{\prime}_{ij}$
$\displaystyle=0$ (64)

up to electric quadrupole-magnetic dipole order. While Maxwell’s equations do
not uniquely specify the form of ${\bf D}$ and ${\bf H}$, additional
requirements from the reciprocity relations and spatial-origin independence
gives Eq. (D) as shown in [25]. Here, the free charge and current $\rho^{f}$
and ${\bf J}^{f}$ are boundary charge and current appearing due to the change
of material properties across the interface that is not described by
$\tilde{T}_{ij}$, $\tilde{U}_{ij}$, and
$S_{ijk}=(a^{\prime}_{ijk}+a^{\prime}_{jki}+a^{\prime}_{kij})/3$.

In $PT$-symmetric systems where $G^{\prime}_{ij}=a_{ijk}=0$,

$\displaystyle\rho^{f}$
$\displaystyle=-B_{j}\partial_{i}\left(G_{ij}-T_{ij}-\frac{1}{2}\epsilon_{jkl}a^{\prime}_{kli}\right)$
$\displaystyle\quad+\frac{i}{2}(\partial_{j}{E}_{k}+\partial_{k}{E}_{j})\partial_{i}(a^{\prime}_{kij}-S_{kij})$
$\displaystyle\quad+\frac{i}{2}{E}_{k}\partial_{i}\partial_{j}(a^{\prime}_{kij}-S_{kij}),$
$\displaystyle J^{f}_{i}$
$\displaystyle=E_{j}\partial_{k}\left[\epsilon_{ikl}(G_{jl}-T_{ji})-\frac{\omega}{2}(a^{\prime}_{jik}-S_{jik})\right],$
(65)

###  D.1 Wave equation

The wave equation up to electric quadrupole/magnetic dipole takes the
following form [21]:

$\displaystyle\left[\delta_{ij}+\epsilon_{0}^{-1}\tilde{\chi}_{ij}+i\mu_{0}nc\sum_{k}\kappa_{k}\tilde{\sigma}_{ijk}+n^{2}(\kappa_{i}\kappa_{j}-\delta_{ij})\right]E_{j}=0$
(66)

where $\tilde{\chi}_{ij}=\chi_{ij}-i\chi^{\prime}_{ij}$ satisfying
$\chi_{ij}=\chi_{ji}$ and $\chi^{\prime}_{ij}=-\chi^{\prime}_{ji}$,
$\kappa_{i}=k_{i}/|{\bf k}|$ is the propagation direction of light, $n$ is the
refractive index, $\tilde{\sigma}_{ijk}=\sigma_{ijk}-i\sigma^{\prime}_{ijk}$
is the complex bulk conductivity coefficient defined by
$\tilde{\sigma}_{ij}({\bf
q})=\tilde{\sigma}_{ij}(0)+\tilde{\sigma}_{ijk}q_{k}+\ldots$, and

$\displaystyle\sigma_{ijk}$
$\displaystyle=i\left[\epsilon_{ikl}G_{jl}+\epsilon_{jkl}G_{il}-\frac{1}{2}\omega(a_{ijk}^{\prime}+a_{jik}^{\prime})\right]=\sigma_{jik},$
$\displaystyle\sigma^{\prime}_{ijk}$
$\displaystyle=i\left[-\epsilon_{ikl}G^{\prime}_{jl}+\epsilon_{jkl}G^{\prime}_{il}+\frac{1}{2}\omega(a_{ijk}-a_{jik})\right]=-\sigma^{\prime}_{jik}.$
(67)

Let us assume $C_{3z}$ and $C_{2x}$ symmetries for simplicity. We further
impose that the bulk Hall response is zero, i.e., $\chi^{\prime}_{ij}=0$, in
order to focus on the magneto-electric and electric-quadrupole effects. For
${\bm{\kappa}}=\pm\hat{z}$, the wave equation is

$\displaystyle\begin{pmatrix}1+\epsilon_{0}^{-1}\chi_{xx}-n^{2}&n\kappa_{z}\mu_{0}c\sigma^{\prime}_{xyz}\\\
-n\kappa_{z}\mu_{0}c\sigma^{\prime}_{xyz}&1+\epsilon_{0}^{-1}\chi_{xx}-n^{2}\end{pmatrix}\begin{pmatrix}E_{x}\\\
E_{y}\end{pmatrix}=0.$ (68)

The refractive index satisfying the wave equation is given by

$\displaystyle
n_{\pm}^{\kappa_{z}}=\sqrt{1+\epsilon_{0}^{-1}\chi_{xx}+(i\mu_{0}c\sigma^{\prime}_{xyz}/2)^{2}}\mp
i\kappa_{z}\mu_{0}c\sigma^{\prime}_{xyz}/2,$ (69)

for circular polarization $\hat{\pm}=\hat{x}+i\hat{y}$.

###  D.2 Reflection and transmission from a single interface

We consider the interface of medium 1 ($z>0$) and medium 2 ($z<0$) with the
surface normal $\hat{z}$. For normal incidence, in the circularly polarized
basis,

$\displaystyle\begin{pmatrix}H_{+}\\\
H_{-}\end{pmatrix}=\frac{1}{\mu_{0}c}\begin{pmatrix}g_{+}&0\\\
0&g_{-}\end{pmatrix}\begin{pmatrix}E_{+}\\\ E_{-}\end{pmatrix}$ (70)

within the media $\mu=1$ or $2$, where
$g_{\pm}=-\mu_{0}c\tilde{T}^{\mu}_{xx}\mp in^{\kappa_{z}}_{\mu\pm}\kappa_{z}$,
$n_{\pm}$ depends on the sign $\kappa_{z}$, and $\kappa_{z}=-1$ for incident
and transmitted light, while $\kappa_{z}=1$ for reflected light. Here,

$\displaystyle\tilde{T}_{xx}$
$\displaystyle=\frac{1}{3}(G_{xx}-G_{zz})-\frac{1}{6}\omega(a^{\prime}_{yzx}-a^{\prime}_{zyx})$
$\displaystyle\quad-i\left[G^{\prime}_{xx}-\frac{1}{2}\omega(a_{yzx}-a_{zyx})\right]$
$\displaystyle=\frac{i}{3}\sigma_{zxy}-\frac{1}{2}\sigma^{\prime}_{xyz}.$ (71)

As we consider light incident from medium $1$ to medium $2$, the electric
field in medium $1$ consists of incident and reflected fields while that in
medium $2$ is the transmitted field.

$\displaystyle{\bf E}_{1}$ $\displaystyle={\bf E}^{i}+{\bf
E}^{r}\equiv(1+r){\bf E}^{i},$ $\displaystyle{\bf E}_{2}$ $\displaystyle={\bf
E}^{t}\equiv t{\bf E}^{i},$ (72)

where

$\displaystyle r=\begin{pmatrix}r_{++}&r_{++}\\\
r_{-+}&r_{--}\end{pmatrix}=\begin{pmatrix}r_{++}&0\\\
0&r_{--}\end{pmatrix},\quad t=1+r$ (73)

by $C_{3z}$ symmetry and the continuity of ${\bf E}$ at the interface.

The ${\bf H}$ field satisfies the boundary condition

$\displaystyle{\bf H}^{t}={\bf H}^{i}+{\bf H}^{r}+\hat{z}\times{\bf j}_{f},$
(74)

where $\hat{z}\times{\bf j}_{f}=(e^{2}/2\pi
h)(\theta^{(z)}_{2}-\theta^{(z)}_{1}){\bf E}$ is the surface current due to
the axion magneto-electric coupling. In terms of ${\bf B}$ fields, the
boundary condition has the form of Eq. (13):

$\displaystyle{\bf B}^{t}$ $\displaystyle={\bf B}^{i}+{\bf
B}^{r}+\mu_{0}\hat{z}\times{\bf j}_{s},$ (75)

where ${\bf j}_{s}={\bf j}_{f}+(\tilde{T}^{2}_{xx}-\tilde{T}^{1}_{xx}){\bf
E}\times\hat{z}$ is the total two-dimensional surface current density. By
solving the boundary condition, we obtain

$\displaystyle r_{++}$
$\displaystyle=\frac{n_{1L}-n_{2L}-i\mu_{0}c\sigma^{s}_{xy}}{n_{1R}+n_{2L}+i\mu_{0}c\sigma^{s}_{xy}},$
$\displaystyle r_{--}$
$\displaystyle=\frac{n_{1R}-n_{2R}+i\mu_{0}c\sigma^{s}_{xy}}{n_{1L}+n_{2R}-i\mu_{0}c\sigma^{s}_{xy}},$
$\displaystyle r_{+-}$ $\displaystyle=r_{-+}=0,$ (76)

where $n_{\mu L}=n_{\mu+}^{-}=n_{\mu-}^{+}$ is the refractive index for left
circularly polarization, and $n_{\mu R}=n_{\mu-}^{-}=n_{\mu+}^{+}$ is the
refractive index for the right circularly polarization, and

$\displaystyle\sigma^{s}_{xy}=\tilde{G}^{(z),2}_{xx}-\tilde{G}^{(z),1}_{xx}=\tilde{T}^{\mu=2}_{xx}-\tilde{T}^{\mu=1}_{xx}+\sigma^{f}_{xy}$
(77)

is the two-dimensional surface conductivity, and $\sigma^{f}_{xy}=(e^{2}/2\pi
h)(\theta^{(z)}_{2}-\theta^{(z)}_{1})$. From the expressions of $r_{++}$ and
$r_{--}$ and Eq. (69), we obtain the Kerr angle

$\displaystyle\varphi_{K}$
$\displaystyle=\tan^{-1}\frac{r_{xy}}{r_{xx}}=\tan^{-1}\left[\frac{-i(r_{++}-r_{--})}{r_{++}+r_{--}}\right].$
(78)

In non-magnetic systems where $G_{ij}=a^{\prime}_{ijk}=0$, the Kerr angle
vanishes systems because

$\displaystyle r_{++}-r_{--}$ $\displaystyle\propto
i\mu_{0}c(\sigma^{\prime\mu=2}_{xyz}-\sigma^{\prime\mu=1}_{xyz})$
$\displaystyle\quad+n_{2R}-n_{2L}-(n_{1R}-n_{1L})$ $\displaystyle=0,$ (79)

where we use that $\tilde{T}^{\mu}_{xx}=-\sigma^{\prime\mu}_{xyz}/2$ when
$\sigma_{ijk}=0$. One may think of this cancellation as a compensation between
bulk and surface responses. The refractive indices are responsible for
circular birefringence in the bulk, while $\delta\sigma^{\prime}_{xyz}$ is
responsible for the surface current that leads to the jump of ${\bf B}$ field
at the surface. Their effects cancel such that there is no net polar Kerr
rotation (i.e., no Kerr rotation at normal incidence), compatible with the
reciprocal relation imposed by time reversal symmetry [53, 54, 29, 30, 31].

In contrast, a nonzero polar Kerr rotation $PT$-symmetric antiferromagnets
occurs because no compensation occurs due to the absence of circular
birefringence in the bulk. By imposing $PT$ symmetry and breaking $T$
symmetry, we get a nonzero Kerr angle from

$\displaystyle r_{xx}$
$\displaystyle=\frac{1}{2}(r_{++}+r_{--})=\frac{(n_{1}-n_{2})(n_{1}+n_{2})-(\mu_{0}c\sigma^{s}_{xy})^{2}}{(\mu_{0}c\sigma^{s}_{xy})^{2}+(n_{1}+n_{2})^{2}},$
$\displaystyle r_{xy}$
$\displaystyle=\frac{1}{2i}(r_{++}-r_{--})=-\frac{2n_{1}(\mu_{0}c\sigma^{s}_{xy})}{(\mu_{0}c\sigma^{s}_{xy})^{2}+(n_{1}+n_{2})^{2}}.$
(80)

Note that our derivation using the ${\bf H}$ field make it manifest that the
magneto-optic reflection has two different origins: bulk-propagation and
surface effects, respectively contained in $T_{xx}$ and $\theta^{(z)}$.

###  D.3 Reflection and transmission from two interfaces

Now we consider three media with $(n_{\mu},m_{\mu})$, where $\mu=1,2,3$, as
shown in Fig. 3b. We consider the limit where the wavelength of light
$\lambda$ is much larger than the sample thickness $d$ and neglect the
variation of the electric field between the top and bottom within the sample.

The Jones reflection matrix of the sample for light incident from medium 1 is
then

$\displaystyle r$
$\displaystyle=r_{T}+t^{\prime}_{T}e^{i\phi}r_{B}e^{i\phi}t_{T}+t^{\prime}_{T}e^{i\phi}r_{B}(e^{i\phi}r^{\prime}_{T}e^{i\phi}r_{B})e^{i\phi}t_{T}$
$\displaystyle\quad+t^{\prime}_{T}e^{i\phi}r_{B}(e^{i\phi}r^{\prime}_{T}e^{i\phi}r_{B})^{2}e^{i\phi}t_{T}+\ldots$
$\displaystyle=r_{T}+e^{2i\phi}t^{\prime}_{T}r_{B}\left(1-e^{2i\phi}r^{\prime}_{T}r_{B}\right)^{-1}t_{T}$
$\displaystyle=r_{T}+e^{2i\phi}(1+r^{\prime}_{T})r_{B}\left(1-e^{2i\phi}r^{\prime}_{T}r_{B}\right)^{-1}(1+r_{T}),$
(81)

where we used $t_{T,B}=1+r_{T,B}$ and $t^{\prime}_{T}=1+r^{\prime}_{T}$, and

$\displaystyle\phi=\frac{n_{2}\omega d}{c}$ (82)

is the complex-valued phase obtained by the propagation across the sample.
Here, the subscript for $r$ and $t$ indicates the top and bottom of the
sample, and the prime indicates the process where the light propagation is
reversed (we follow the notation in Ref. [13]).

For transmission, the Jones matrix is

$\displaystyle t$
$\displaystyle=t_{B}e^{i\phi}t_{T}+t_{B}(e^{i\phi}r^{\prime}_{T}e^{i\phi}r_{B})e^{i\phi}t_{T}$
$\displaystyle\quad+t_{B}(e^{i\phi}r^{\prime}_{T}e^{i\phi}r_{B})(e^{i\phi}r^{\prime}_{T}e^{i\phi}r_{B})t_{T}+\ldots$
$\displaystyle=t_{B}\left(1-e^{2i\phi}r^{\prime}_{T}r_{B}\right)^{-1}e^{i\phi}t_{T}$
$\displaystyle=(1+r_{B})\left(1-e^{2i\phi}r^{\prime}_{T}r_{B}\right)^{-1}e^{i\phi}(1+r_{T}).$
(83)

The expressions above show that both $r$ and $t$ can be obtained from the
reflective Jones matrices at the top and bottom interfaces. We suppose that
each media has $PT$, $C_{3z}$-, and $M_{x}T$ symmetries, which is the setup in
main text. Then we have

$\displaystyle(r_{T})_{xx}$
$\displaystyle=\frac{(n_{1}-n_{2})(n_{1}+n_{2})-(\mu_{0}c\sigma^{T}_{xy})^{2}}{(\mu_{0}c\sigma^{T}_{xy})^{2}+(n_{1}+n_{2})^{2}},$
$\displaystyle(r_{T})_{xy}$
$\displaystyle=-\frac{2n_{1}(\mu_{0}c\sigma^{T}_{xy})}{(\mu_{0}c\sigma^{s}_{xy})^{2}+(n_{1}+n_{2})^{2}},$
$\displaystyle(r^{\prime}_{T})_{xx}$
$\displaystyle=\frac{-(n_{1}-n_{2})(n_{1}+n_{2})-(\mu_{0}c\sigma^{T}_{xy})^{2}}{(\mu_{0}c\sigma^{T}_{xy})^{2}+(n_{1}+n_{2})^{2}},$
$\displaystyle(r^{\prime}_{T})_{xy}$
$\displaystyle=-\frac{2n_{2}(\mu_{0}c\sigma^{T}_{xy})}{(\mu_{0}c\sigma^{T}_{xy})^{2}+(n_{1}+n_{2})^{2}},$
$\displaystyle(r_{B})_{xx}$
$\displaystyle=\frac{(n_{2}-n_{3})(n_{2}+n_{3})-(\mu_{0}c\sigma^{B}_{xy})^{2}}{(\mu_{0}c\sigma^{B}_{xy})^{2}+(n_{2}+n_{3})^{2}},$
$\displaystyle(r_{B})_{xy}$
$\displaystyle=\frac{2n_{2}(\mu_{0}c\sigma^{B}_{xy})}{(\mu_{0}c\sigma^{B}_{xy})^{2}+(n_{2}+n_{3})^{2}},$
(84)

where the superscript $T$ or $B$ for the surface Hall conductivity
$\sigma_{xy}$ indicates the top and bottom surfaces.

##  Appendix E Kerr angle and Stokes parameters

The reflective circular dichroism and Kerr rotation angle can be defined with
Stokes parameters $s_{i=0,1,2,3}$ for reflected light by

$\displaystyle{\rm RCD}$ $\displaystyle\equiv\frac{s_{3}}{s_{0}},$
$\displaystyle\vartheta_{K}$
$\displaystyle\equiv\frac{1}{2}\tan^{-1}\frac{s_{2}}{s_{1}},$ (85)

where Stokes parameters for reflected light are

$\displaystyle s_{0}$ $\displaystyle=I^{r}_{+}+I^{r}_{-},$ $\displaystyle
s_{1}$ $\displaystyle=I^{r}_{x}-I^{r}_{y},$ $\displaystyle s_{2}$
$\displaystyle=I^{r}_{x+y}-I^{r}_{x-y},$ $\displaystyle s_{3}$
$\displaystyle=I^{r}_{+}-I^{r}_{-},$ (86)

where $\hat{\pm}=\hat{x}\pm i\hat{y}$. While this definition of the complex
Kerr angle is different from the more popular
$\phi_{K}=r_{xy}/r_{xx}=\varphi_{K}+i\eta_{K}$ we use in the main text, it has
the advantage that it can be obtained simply by measuring the intensity of the
linearly polarized light. In our case where $C_{3z}$ symmetry is present, the
reflective coefficients in circularly polarized basis are

$\displaystyle r_{++}$
$\displaystyle=\frac{1}{2}\left[r_{xx}+r_{yy}+i(r_{xy}-r_{yx})\right]=r_{xx}+ir_{xy},$
$\displaystyle r_{--}$
$\displaystyle=\frac{1}{2}\left[r_{xx}+r_{yy}-i(r_{xy}-r_{yx})\right]=r_{xx}-ir_{xy},$
$\displaystyle r_{+-}$
$\displaystyle=\frac{1}{2}\left[r_{xx}-r_{yy}-i(r_{xy}+r_{yx})\right]=0,$
$\displaystyle r_{-+}$
$\displaystyle=\frac{1}{2}\left[r_{xx}-r_{yy}+i(r_{xy}+r_{yx})\right]=0.$ (87)

It follows that ${\rm RCD}=2{\rm
Im}\left[r_{xx}^{*}r_{xy}\right]/(|r_{xx}|^{2}+|r_{xy}|^{2})\simeq 2{\rm
Im}[r_{xy}/r_{xx}]=2\eta_{K}$ and $\vartheta_{K}=\frac{1}{2}\tan^{-1}[2{\rm
Re}\left(r_{xx}r_{yx}^{*}\right)/(|r_{xx}|^{2}-|r_{xy}|^{2})]\simeq{\rm
Re}[r_{xy}/r_{xx}]=\varphi_{K}$ for small $|r_{xy}/r_{xx}|$.

##  Appendix F Tight-binding model of MnBi2Te4

We begin with the lattice version of the three-dimensional low-energy model of
MnBi2Te4 at $\Gamma=(0,0,0)$ presented in Ref. [38]. We consider the four
basis states.

$\displaystyle|P1_{z}^{+},\uparrow\rangle,\quad|P2_{z}^{-},\uparrow\rangle,\quad|P1_{z}^{+},\downarrow\rangle,\quad|P2_{z}^{-},\downarrow\rangle,$
(88)

where $P1$ and $P2$ states originates from the $p$ orbitals in Bi and Te,
respectively, the sign $\pm$ indicates the inversion parities, and the arrows
indicate the spin-$z$ direction. The symmetry operators in the nonmagnetic
state are time reversal $T$, inversion $P$, threefold rotation $C_{3z}$, and
twofold rotation $C_{2x}$.

$\displaystyle T=is_{y}K,\quad P=\tau_{z},\quad
C_{3z}=\exp\left[-i\frac{\pi}{3}s_{z}\right],\quad C_{2x}=-is_{x}.$ (89)

Because of $PT=is_{y}\tau_{z}K$ symmetry, only the following five Gamma
matrices are allowed in the Hamiltonian in addition to the overall energy
shift proportional to the identity matrix $\Gamma_{0}$.

$\displaystyle\Gamma_{1}$ $\displaystyle=s_{x}\tau_{x}\quad(+,-),$
$\displaystyle\Gamma_{2}$ $\displaystyle=s_{y}\tau_{x}\quad(-,-),$
$\displaystyle\Gamma_{3}$ $\displaystyle=s_{z}\tau_{x}\quad(-,-),$
$\displaystyle\Gamma_{4}$ $\displaystyle=s_{0}\tau_{y}\quad(+,-),$
$\displaystyle\Gamma_{5}$ $\displaystyle=s_{0}\tau_{z}\quad(+,+).$ (90)

Here, the pair signs show the commutation ($+$) and anticommutation ($-$)
relations with $C_{2x}$ and $P$ in order, i.e.,
$(\epsilon_{C_{2x}},\epsilon_{P})$ where $M\Gamma_{i}=\epsilon_{M}\Gamma_{i}M$
for $i=1,\ldots,5$.

The low-energy effective Hamiltonian up to second order in ${\bf k}$ has the
form

$\displaystyle h_{\rm eff}$ $\displaystyle=\epsilon({\bf
k})s_{0}\tau_{0}+A_{1}k_{z}s_{z}\tau_{x}+A_{2}(k_{x}s_{x}+k_{y}s_{y})\tau_{x}+M({\bf
k})\tau_{z}$ $\displaystyle=\epsilon({\bf
k})\Gamma_{0}+A_{1}k_{z}\Gamma_{3}+A_{2}(k_{x}\Gamma_{1}+k_{y}\Gamma_{2})+M({\bf
k})\Gamma_{5},$ (91)

where

$\displaystyle\epsilon({\bf k})$
$\displaystyle=C_{0}+C_{1}k_{z}^{2}+C_{2}(k_{x}^{2}+k_{y}^{2}),$
$\displaystyle M({\bf k})$
$\displaystyle=M_{0}+M_{1}k_{z}^{2}+M_{2}(k_{x}^{2}+k_{y}^{2}).$ (92)

Let us find a tight-binding Hamiltonian on the two-dimensional triangular
lattice that leads to the above low-energy Hamiltonian. We suppose that each
lattice site has four degrees of freedom given by Eq. (88).

$\displaystyle\hat{H}_{\rm
TB}=\sum_{i,\alpha,\beta}\hat{c}^{\dagger}_{i\alpha}(h_{0})_{\alpha\beta}\hat{c}_{i\beta}-\sum_{\braket{i,j},\alpha,\beta}\hat{c}^{\dagger}_{i\alpha}t^{ij}_{\alpha\beta}\hat{c}_{j\beta}.$
(93)

Here, the onsite Hamiltonian satisfying all the symmetries of the nonmagnetic
state is

$\displaystyle h_{0}=e_{0}\Gamma_{0}+e_{5}\Gamma_{5}.$ (94)

Along the $z$ direction, the nearest-neighbor hopping matrices are

$\displaystyle T_{4}\equiv t^{j+{\bf
a}_{4},j}=t^{z}_{0}\Gamma_{0}+it^{z}_{3}\Gamma_{3}+t^{z}_{5}\Gamma_{5},$ (95)

where ${\bf a}_{4}=(0,0,a_{z})$, $a_{z}$ is the inter-layer lattice parameter.
the form of $T_{4}$ is constrained by the following symmetry conditions

$\displaystyle TT_{4}T^{-1}$ $\displaystyle=T_{4},$ $\displaystyle
C_{3z}T_{4}C_{3z}^{-1}$ $\displaystyle=T_{4},$ $\displaystyle
C_{2x}T_{4}C_{2x}^{-1}$ $\displaystyle=T_{4}^{\dagger},$ $\displaystyle
PT_{4}P^{-1}$ $\displaystyle=T_{4}^{\dagger},$ (96)

For the in-plane directions, the hopping matrices are

$\displaystyle T_{1}$ $\displaystyle\equiv t^{j+{\bf
a}_{1},j}=t_{0}\Gamma_{0}+it_{1}\Gamma_{1}+it_{4}\Gamma_{4}+t_{5}\Gamma_{5}=(t^{j,j+{\bf
a}_{1}})^{\dagger},$ $\displaystyle T_{2}$ $\displaystyle\equiv t^{j+{\bf
a}_{2},j}=C_{3z}T_{1}C_{3z}^{-1}$
$\displaystyle=t_{0}\Gamma_{0}+it_{1}\frac{-\Gamma_{1}+\sqrt{3}\Gamma_{2}}{2}+it_{4}\Gamma_{4}+t_{5}\Gamma_{5}=(t^{j,j+{\bf
a}_{2}})^{\dagger},$ $\displaystyle T_{3}$ $\displaystyle\equiv t^{j+{\bf
a}_{3},j}=C_{3z}T_{2}C_{3z}^{-1}$
$\displaystyle=t_{0}\Gamma_{0}+it_{1}\frac{-\Gamma_{1}-\sqrt{3}\Gamma_{2}}{2}+it_{4}\Gamma_{4}+t_{5}\Gamma_{5}=(t^{j,j+{\bf
a}_{3}})^{\dagger}$ (97)

along the in-plane directions, where ${\bf a}_{1}=(a,0,0)$, ${\bf
a}_{2}=C_{3z}{\bf a}_{1}$, ${\bf a}_{3}=C_{3z}{\bf a}_{2}$, $a$ is the in-
plane lattice parameter, and the form of $T_{1}$ is constrained by the
following symmetry conditions

$\displaystyle TT_{1}T^{-1}$ $\displaystyle=T_{1},$ $\displaystyle
C_{2x}T_{1}C_{2x}^{-1}$ $\displaystyle=T_{1},$ $\displaystyle PT_{1}P^{-1}$
$\displaystyle=T_{1}^{\dagger},$ (98)

and $T_{i}$s are related by $C_{3z}$ because of $C_{3z}$ symmetry imposing,
for example,

$\displaystyle\hat{c}^{\dagger}_{j+{\bf
a}_{2}\alpha}(T_{2})_{\alpha\beta}\hat{c}_{j\beta}$
$\displaystyle=\hat{C}_{3z}\hat{c}^{\dagger}_{j+{\bf
a}_{1}\alpha}(T_{1})_{\alpha\beta}\hat{c}_{j\beta}\hat{C}_{3z}^{-1}$
$\displaystyle=\left(\hat{C}_{3z}\hat{c}^{\dagger}_{j+{\bf
a}_{1}\alpha}\hat{C}_{3z}^{-1}\right)(T_{1})_{\alpha\beta}\left(\hat{C}_{3z}\hat{c}_{j\beta}\hat{C}_{3z}^{-1}\right)$
$\displaystyle=\hat{c}^{\dagger}_{j+{\bf
a}_{2}\gamma}(C_{3z})_{\gamma\alpha}(T_{1})_{\alpha\beta}\hat{c}_{j\delta}(C_{3z}^{*})_{\delta\beta}$
$\displaystyle=\hat{c}^{\dagger}_{j+{\bf
a}_{2}\alpha}(C_{3z}T_{1}C_{3z}^{-1})_{\alpha\beta}\hat{c}_{j\beta}.$ (99)

In momentum space, the tight-binding Hamiltonian becomes

$\displaystyle h_{\rm TB}$
$\displaystyle=h_{0}-(T_{1}e^{-ik_{1}a}+T_{2}e^{-ik_{2}a}+T_{3}e^{-ik_{3}a}+T_{4}e^{-ik_{4}a_{z}}+h.c.)$
$\displaystyle=\left[e_{0}-2t_{0}(\cos k_{1}a+\cos k_{2}a+\cos
k_{3}a)-2t_{0}^{z}\cos k_{4}a_{z}\right]\Gamma_{0}$ $\displaystyle-t_{1}(2\sin
k_{1}a-\sin k_{2}a-\sin k_{3}a)\Gamma_{1}$ $\displaystyle-\sqrt{3}t_{1}(\sin
k_{2}a-\sin k_{3}a)\Gamma_{2}-2t_{3}^{z}\sin k_{4}a_{z}\Gamma_{3}$
$\displaystyle-2t_{4}(\sin k_{1}a+\sin k_{2}a+\sin k_{3}a)\Gamma_{4}$
$\displaystyle+\left[e_{5}-2t_{5}(\cos k_{1}a+\cos k_{2}a+\cos
k_{3}a)-2t_{5}^{z}\cos k_{4}a_{z}\right]\Gamma_{5},$ (100)

where

$\displaystyle k_{1}$ $\displaystyle=k_{x},$ $\displaystyle k_{2}$
$\displaystyle=\frac{1}{2}\left(-k_{x}+\sqrt{3}k_{y}\right),$ $\displaystyle
k_{3}$ $\displaystyle=\frac{1}{2}\left(-k_{x}-\sqrt{3}k_{y}\right),$
$\displaystyle k_{4}$ $\displaystyle=k_{z}.$ (101)

By expanding the tight-binding Hamiltonian up to second order in ${\bf k}$, we
obtain

$\displaystyle h_{\rm TB}$
$\displaystyle=\left[e_{0}-6t_{0}-2t_{0}^{z}+\frac{3}{2}t_{0}a^{2}(k_{x}^{2}+k_{y}^{2})+t_{0}^{z}a_{z}^{2}k_{z}^{2}\right]\Gamma_{0}$
$\displaystyle\quad-3t_{1}ak_{x}\Gamma_{1}-3t_{1}ak_{y}\Gamma_{2}-2t_{3}^{z}ak_{z}\Gamma_{2}$
$\displaystyle+\left[e_{5}-6t_{5}-2t_{5}^{z}+\frac{3}{2}t_{5}a^{2}(k_{x}^{2}+k_{y}^{2})+t_{5}^{z}a_{z}^{2}k_{z}^{2}\right]\Gamma_{5}$
$\displaystyle\quad+O(k^{3}).$ (102)

Comparing this with $h_{\rm eff}$, we find

$\displaystyle e_{0}$ $\displaystyle=C_{0}+2C_{1}/a_{z}^{2}+4C_{2}/a^{2},$
$\displaystyle e_{5}$ $\displaystyle=M_{0}+2M_{1}/a_{z}^{2}+4M_{2}/a^{2},$
$\displaystyle t_{0}$ $\displaystyle=\frac{2C_{2}}{3a^{2}},$ $\displaystyle
t_{0}^{z}$ $\displaystyle=\frac{C_{1}}{a_{z}^{2}},$ $\displaystyle t_{1}$
$\displaystyle=-\frac{A_{2}}{3a},$ $\displaystyle t_{3}^{z}$
$\displaystyle=-\frac{A_{1}}{2a_{z}},$ $\displaystyle t_{5}$
$\displaystyle=\frac{2M_{2}}{3a^{2}},$ $\displaystyle t_{5}^{z}$
$\displaystyle=\frac{M_{1}}{a_{z}^{2}}.$ (103)

We use the parameters derived in Ref. [38] with a modification, $C_{2}=0$, we
take in order to obtain an insulating filling:

$\displaystyle C_{0}$ $\displaystyle=-0.0048\;{\rm eV},$ $\displaystyle C_{1}$
$\displaystyle=2.7232\;{\rm eV\AA^{2}},$ $\displaystyle C_{2}$
$\displaystyle=0\;{\rm eV\AA^{2}},$ $\displaystyle M_{0}$
$\displaystyle=-0.1165\;{\rm eV},$ $\displaystyle M_{1}$
$\displaystyle=11.9048\;{\rm eV\AA^{2}},$ $\displaystyle M_{2}$
$\displaystyle=9.4048\;{\rm eV\AA^{2}},$ $\displaystyle A_{1}$
$\displaystyle=2.7023\;{\rm eV\AA},$ $\displaystyle A_{2}$
$\displaystyle=3.1964\;{\rm eV\AA},$ $\displaystyle a$
$\displaystyle=4.334\;{\rm\AA},$ $\displaystyle a_{z}$
$\displaystyle=\frac{1}{3}c=\frac{1}{3}40.91\;{\rm\AA}=13.64\;{\rm\AA}.$ (104)

The spin part is described by

$\displaystyle h_{\rm spin}$ $\displaystyle=-\mu_{B}{\bf B}\cdot{\bf
s}\tau_{0},$ (105)

where the Bohr magneton is

$\displaystyle\mu_{B}$ $\displaystyle=3.8099\frac{e}{\hbar}\;{\rm eV\AA^{2}}.$
(106)

We consider the antiferromagnetic state of a few-layer MnBi2Te4 with layer-
alternating out-of-plane moments. The antiferromagnetic moment is described by

$\displaystyle h_{\rm AFM}$ $\displaystyle=ml_{z}s_{z}\tau_{0},$ (107)

where $l_{z}$ is a Pauli matrix in the sublattice (i.e., layer) space. We use
$m=0.03$ eV for our calculations. The full $8\times 8$ tight-binding
Hamiltonian in momentum space is then

$\displaystyle h_{\rm TB}$ $\displaystyle=h_{\rm para}+h_{\rm AFM}$
$\displaystyle=\left[e_{0}-2t_{0}(\cos k_{1}a+\cos k_{2}a+\cos
k_{3}a)\right]l_{0}\Gamma_{0}$ $\displaystyle\quad-2t_{0}^{z}\cos
k_{4}a_{z}l_{x}\Gamma_{0}$ $\displaystyle\quad-t_{1}(2\sin k_{1}a-\sin
k_{2}a-\sin k_{3}a)l_{0}\Gamma_{1}$ $\displaystyle\quad-\sqrt{3}t_{1}(\sin
k_{2}a-\sin k_{3}a)l_{0}\Gamma_{2}+2t_{3}^{z}\sin k_{4}a_{z}l_{y}\Gamma_{3}$
$\displaystyle\quad-2t_{4}(\sin k_{1}a+\sin k_{2}a+\sin
k_{3}a)l_{0}\Gamma_{4}$ $\displaystyle\quad+\left[e_{5}-2t_{5}(\cos
k_{1}a+\cos k_{2}a+\cos k_{3}a)\right]\Gamma_{5}$
$\displaystyle\quad-2t_{5}^{z}\cos k_{4}a_{z}l_{x}\Gamma_{5}$
$\displaystyle\quad+ml_{z}\Gamma_{12},$ (108)

where $\Gamma_{12}=\Gamma_{1}\Gamma_{2}/2i=s_{z}\tau_{0}$.

## References

  * [1] Wan, X., Turner, A. M., Vishwanath, A. & Savrasov, S. Y.  Topological semimetal and Fermi-arc surface states in the electronic structure of pyrochlore iridates.  _Phys. Rev. B_ 83, 205101 (2011).
  * [2] Hughes, T. L., Prodan, E. & Bernevig, B. A.  Inversion-symmetric topological insulators.  _Phys. Rev. B_ 83, 245132 (2011).
  * [3] Turner, A. M., Zhang, Y., Mong, R. S. K. & Vishwanath, A.  Quantized response and topology of magnetic insulators with inversion symmetry.  _Phys. Rev. B_ 85, 165120 (2012).
  * [4] Teo, J. C. Y., Fu, L. & Kane, C. L.  Surface states and topological invariants in three-dimensional topological insulators: Application to Bi1-xSbx.  _Phys. Rev. B_ 78, 045426 (2008).
  * [5] Fang, C. & Fu, L.  New classes of three-dimensional topological crystalline insulators: Nonsymmorphic and magnetic.  _Phys. Rev. B_ 91, 161105(R) (2015).
  * [6] Varjas, D., de Juan, F. & Lu, Y.-M.  Bulk invariants and topological response in insulators and superconductors with nonsymmorphic symmetries.  _Phys. Rev. B_ 92, 195116 (2015).
  * [7] Fu, L., Kane, C. L. & Mele, E. J.  Topological insulators in three dimensions.  _Phys. Rev. Lett._ 98, 106803 (2007).
  * [8] Qi, X.-L., Hughes, T. L. & Zhang, S.-C.  Topological field theory of time-reversal invariant insulators.  _Phys. Rev. B_ 78, 195424 (2008).
  * [9] Essin, A. M., Moore, J. E. & Vanderbilt, D.  Magnetoelectric polarizability and axion electrodynamics in crystalline insulators.  _Phys. Rev. Lett._ 102, 146805 (2009).
  * [10] Liu, C. _et al._ Robust axion insulator and Chern insulator phases in a two-dimensional antiferromagnetic topological insulator.  _Nat. Mater._ 19, 522–527 (2020).
  * [11] Deng, Y. _et al._ Quantum anomalous Hall effect in intrinsic magnetic topological insulator MnBi2Te4.  _Science_ 367, 895–900 (2020).
  * [12] Gao, A. _et al._ Layer Hall effect in a 2d topological axion antiferromagnet.  _Nature_ 595, 521–525 (2021).
  * [13] Tse, W.-K. & MacDonald, A. H.  Giant magneto-optical Kerr effect and universal Faraday effect in thin-film topological insulators.  _Phys. Rev. Lett._ 105, 057401 (2010).
  * [14] Maciejko, J., Qi, X.-L., Drew, H. D. & Zhang, S.-C.  Topological quantization in units of the fine structure constant.  _Phys. Rev. Lett._ 105, 166803 (2010).
  * [15] Wu, L. _et al._ Quantized Faraday and Kerr rotation and axion electrodynamics of a 3D topological insulator.  _Science_ 354, 1124–1127 (2016).
  * [16] Okada, K. N. _et al._ Terahertz spectroscopy on Faraday and Kerr rotations in a quantum anomalous Hall state.  _Nat. Commun._ 7, 1–6 (2016).
  * [17] Malashevich, A., Souza, I., Coh, S. & Vanderbilt, D.  Theory of orbital magnetoelectric response.  _New J. Phys._ 12, 053032 (2010).
  * [18] Essin, A. M., Turner, A. M., Moore, J. E. & Vanderbilt, D.  Orbital magnetoelectric coupling in band insulators.  _Phys. Rev. B_ 81, 205104 (2010).
  * [19] Malashevich, A. & Souza, I.  Band theory of spatial dispersion in magnetoelectrics.  _Phys. Rev. B_ 82, 245118 (2010).
  * [20] Hornreich, R. & Shtrikman, S.  Theory of gyrotropic birefringence.  _Phys. Rev._ 171, 1065 (1968).
  * [21] Raab, R. E. & De Lange, O. L.  _Multipole theory in electromagnetism: classical, quantum, and symmetry aspects, with applications_ , vol. 128 (OUP Oxford, 2004).
  * [22] Qiu, Z. & Bader, S. D.  Surface magneto-optic Kerr effect.  _Rev. Sci. Instrum_ 71, 1243–1255 (2000).
  * [23] Sivadas, N., Okamoto, S. & Xiao, D.  Gate-controllable magneto-optic Kerr effect in layered collinear antiferromagnets.  _Phys. Rev. Lett._ 117, 267203 (2016).
  * [24] Dzyaloshinskii, I. & Papamichail, E.  Nonreciprocal optical rotation in antiferromagnets.  _Phys. Rev. Lett._ 75, 3004 (1995).
  * [25] Graham, E. & Raab, R.  Macroscopic theory of reflection from antiferromagnetic Cr2O3.  _J. Phys. Condens. Matter_ 9, 1863 (1997).
  * [26] Graham, E. & Raab, R.  Surface magneto-optical effects in cubic antiferromagnets.  _Phys. Rev. B_ 62, 9561 (2000).
  * [27] Graham, E. & Raab, R.  The role of the macroscopic surface density of bound electric dipole moment in reflection.  _Proc. R. Soc. A_ 457, 471–484 (2001).
  * [28] Orenstein, J.  Optical nonreciprocity in magnetic structures related to high-$T_{c}$ superconductors.  _Phys. Rev. Lett._ 107, 067002 (2011).
  * [29] Halperin, B.  The hunt for anyon superconductivity.  In _The Physics and Chemistry of Oxide Superconductors_ , 439–450 (Springer, 1992).
  * [30] Shelankov, A. & Pikus, G.  Reciprocity in reflection and transmission of light.  _Phys. Rev. B_ 46, 3326 (1992).
  * [31] Armitage, N. P.  Constraints on Jones transmission matrices from time-reversal invariance and discrete spatial symmetries.  _Phys. Rev. B_ 90, 035135 (2014).
  * [32] Barron, L. D.  _Molecular light scattering and optical activity_ (Cambridge University Press, 2009).
  * [33] Marzari, N., Mostofi, A. A., Yates, J. R., Souza, I. & Vanderbilt, D.  Maximally localized Wannier functions: Theory and applications.  _Rev. Mod. Phys._ 84, 1419 (2012).
  * [34] Liu, Z. & Wang, J.  Anisotropic topological magnetoelectric effect in axion insulators.  _Phys. Rev. B_ 101, 205130 (2020).
  * [35] Spaldin, N. A., Fechner, M., Bousquet, E., Balatsky, A. & Nordström, L.  Monopole-based formalism for the diagonal magnetoelectric response.  _Phys. Rev. B_ 88, 094429 (2013).
  * [36] Li, J. _et al._ Intrinsic magnetic topological insulators in van der waals layered MnBi2Te4-family materials.  _Sci. Adv._ 5, eaaw5685 (2019).
  * [37] Otrokov, M. M. _et al._ Prediction and observation of an antiferromagnetic topological insulator.  _Nature_ 576, 416–422 (2019).
  * [38] Zhang, D. _et al._ Topological axion states in the magnetic insulator MnBi2Te4 with the quantized magnetoelectric effect.  _Phys. Rev. Lett._ 122, 206401 (2019).
  * [39] Sekine, A. & Nomura, K.  Axion electrodynamics in topological materials.  _J. Appl. Phys._ 129, 141101 (2021).
  * [40] He, Q. L., Hughes, T. L., Armitage, N. P., Tokura, Y. & Wang, K. L.  Topological spintronics and magnetoelectronics.  _Nat. Mater._ 21, 15–23 (2022).
  * [41] Bernevig, B. A., Felser, C. & Beidenkopf, H.  Progress and prospects in magnetic topological materials.  _Nature_ 603, 41–51 (2022).
  * [42] Yang, S. _et al._ Odd-even layer-number effect and layer-dependent magnetic phase diagrams in mnbi 2 te 4.  _Phys. Rev. X_ 11, 011003 (2021).
  * [43] Zhang, R.-X., Wu, F. & Sarma, S. D.  Möbius insulator and higher-order topology in MnBi2nTe3n+1.  _Phys. Rev. Lett._ 124, 136407 (2020).
  * [44] Lee, S.-Y., Jeong, T.-Y., Jung, S. & Yee, K.-J.  Refractive index dispersion of hexagonal boron nitride in the visible and near-infrared.  _Phys. Status Solidi B_ 256, 1800417 (2019).
  * [45] Zaitsev, A. M.  _Optical properties of diamond: a data handbook_ (Springer Science & Business Media, 2013).
  * [46] Krichevtsov, B., Pavlov, V., Pisarev, R. & Gridnev, V.  Spontaneous non-reciprocal reflection of light from antiferromagnetic Cr2O3.  _J. Phys. Condens. Matter_ 5, 8233 (1993).
  * [47] Ahn, J., Guo, G.-Y., Nagaosa, N. & Vishwanath, A.  Riemannian geometry of resonant optical responses.  _Nat. Phys._ 18, 290–295 (2022).
  * [48] Li, R., Wang, J., Qi, X.-L. & Zhang, S.-C.  Dynamical axion field in topological magnetic insulators.  _Nat. Phys._ 6, 284–288 (2010).
  * [49] Gao, Y. & Xiao, D.  Nonreciprocal directional dichroism induced by the quantum metric dipole.  _Phys. Rev. Lett._ 122, 227402 (2019).
  * [50] Lapa, M. F. & Hughes, T. L.  Semiclassical wave packet dynamics in nonuniform electric fields.  _Phys. Rev. B_ 99, 121111 (2019).
  * [51] Zhong, S., Moore, J. E. & Souza, I.  Gyrotropic magnetic effect and the magnetic moment on the Fermi surface.  _Phys. Rev. Lett._ 116, 077201 (2016).
  * [52] Ma, J. & Pesin, D.  Chiral magnetic effect and natural optical activity in metals with or without Weyl points.  _Phys. Rev. B_ 92, 235205 (2015).
  * [53] Agranovich, V. & Ydson, V.  Transition layer effects in gyrotropic and nongyrotropic media.  _Opt. Commun._ 5, 422–424 (1972).
  * [54] Agranovich, V. & Yudson, V.  On phenomenological electrodynamics of gyrotropic media.  _Opt. Commun._ 9, 58–60 (1973).

###  Supplementary Note 1 Gyrotropic birefringence and natural optical
activity

Here we review the theoretical analysis of gyrotropic birefringence and
natural optical activity. This is an analysis of the bulk current response as
opposed to the surface current response that we present in the main text.
Therefore, axion electrodyamics is not described by the approach in this
section [1]. Below we mainly follow the work of Malashevich and Souza [1].
Along the way, we compare our notations with the notations in that paper [see
Eqs. (Supplementary Note 1) and (Supplementary Note 1)].

We consider a three-dimensional homogeneous system without boundaries. By the
multipole expansion up to electric-quadrupole/magnetic-dipole order, the bulk
current density is written as

$\displaystyle J_{i}$
$\displaystyle=\dot{P}_{i}-\frac{1}{2}\sum_{j}\partial_{j}\dot{Q}_{ij}+\sum_{jk}\epsilon_{ijk}\partial_{j}m_{k}$
$\displaystyle=-i\omega\tilde{\alpha}_{ij}E_{j}-i\bigg{[}i\sum_{l}(\tilde{G}_{il}\epsilon_{ljk}+\epsilon_{ikl}\tilde{\mathfrak{G}}_{lj})+\frac{\omega}{2}(\tilde{a}_{ijk}-\tilde{\mathfrak{a}}_{ikj})\bigg{]}\partial_{k}E_{j}$
$\displaystyle=\sum_{j}\sigma_{ij}E_{j}-i\sum_{j,k}\sigma_{ijk}\partial_{k}E_{j},$
(S1)

where we use that induced multipole moments are

$\displaystyle P_{i}$
$\displaystyle=\tilde{\chi}_{ij}E_{j}+\frac{1}{2}\tilde{a}_{ijk}\nabla_{k}E_{j}+\tilde{G}_{ij}B_{j}$
$\displaystyle Q_{ij}$ $\displaystyle=\tilde{\mathfrak{a}}_{ijk}E_{k}$
$\displaystyle M_{i}$ $\displaystyle=\tilde{\mathfrak{G}}_{ij}E_{j}$ (S2)

and that response functions do not change spatially. We are interested in the
first order in momentum.

$\displaystyle\sigma_{ijk}=i\sum_{l}(\tilde{G}_{il}\epsilon_{ljk}+\epsilon_{ikl}\tilde{\mathfrak{G}}_{lj})+\frac{\omega}{2}(\tilde{a}_{ijk}-\tilde{\mathfrak{a}}_{ikj})$
(S3)

After symmetrizing the first two indices, we have

$\displaystyle\sigma^{S}_{ijk}$
$\displaystyle=\frac{1}{2}\left(\sigma_{ijk}+\sigma_{jik}\right)$
$\displaystyle=i\sum_{l}(\alpha_{il}\epsilon_{ljk}+\epsilon_{ikl}\alpha_{jl})+\omega\gamma_{ijk},$
(S4)

and, similarly, we have the anti-symmetrized part

$\displaystyle\sigma^{A}_{ijk}$
$\displaystyle=\frac{1}{2}\left(\sigma_{ijk}-\sigma_{jik}\right)$
$\displaystyle=i\sum_{l}(\beta_{il}\epsilon_{ljk}-\epsilon_{ikl}\beta_{jl})+\omega\xi_{ijk},$
(S5)

where we define

$\displaystyle\alpha_{ij}$ $\displaystyle=G_{ij},$ $\displaystyle\beta_{ij}$
$\displaystyle=-iG^{\prime}_{ij},$ $\displaystyle\gamma_{ijk}$
$\displaystyle=-\frac{i}{2}(a^{\prime}_{ijk}+a^{\prime}_{jik}),$
$\displaystyle\xi_{ijk}$ $\displaystyle=\frac{1}{2}(a_{ijk}-a_{jik}),$ (S6)

and

$\displaystyle a_{ijk}$
$\displaystyle=\frac{\tilde{a}_{ijk}+\tilde{\mathfrak{a}}_{jki}}{2},$
$\displaystyle a^{\prime}_{ijk}$
$\displaystyle=i\frac{\tilde{a}_{ijk}-\tilde{\mathfrak{a}}_{jki}}{2},$
$\displaystyle G_{ij}$
$\displaystyle=\frac{\tilde{G}_{ij}+\tilde{\mathfrak{G}}_{ji}}{2},$
$\displaystyle G^{\prime}_{ij}$
$\displaystyle=i\frac{\tilde{G}_{ij}-\tilde{\mathfrak{G}}_{ji}}{2}.$ (S7)

The conductivity tensor is related to the origin-independent combinations of
the magneto-electric coupling and electric quadrupole susceptibility by

$\displaystyle\sigma^{S}_{ijk}$
$\displaystyle=i\sum_{l}(\tilde{\alpha}_{il}\epsilon_{ljk}+\epsilon_{ikl}\tilde{\alpha}_{jl})+\omega\tilde{\gamma}_{ijk},$
$\displaystyle\sigma^{A}_{ijk}$
$\displaystyle=i\sum_{l}(\tilde{\beta}_{il}\epsilon_{ljk}-\epsilon_{ikl}\tilde{\beta}_{jl}),$
(S8)

where

$\displaystyle\tilde{\alpha}_{ij}$
$\displaystyle=\frac{1}{3i}\sum_{l}\epsilon_{jkl}\sigma_{ikl}^{S}=\alpha_{ij}-\frac{1}{3}\delta_{ij}\sum_{k}\alpha_{kk}+\frac{\omega}{3i}\sum_{k,l}\epsilon_{jkl}\gamma_{ikl}$
$\displaystyle=T_{ij},$ $\displaystyle\tilde{\gamma}_{ijk}$
$\displaystyle=\frac{1}{3\omega}(\sigma_{ijk}^{S}+\sigma_{jki}^{S}+\sigma_{kij}^{S})=\frac{1}{3\omega}(\gamma_{ijk}+\gamma_{jki}+\gamma_{kij})$
$\displaystyle=-iS_{ijk},$ $\displaystyle\tilde{\beta}_{ij}$
$\displaystyle=\frac{1}{4i}\sum_{k,l}\epsilon_{jkl}(2\sigma_{ikl}^{A}-\sigma_{kli}^{A})=\beta_{ij}+\frac{1}{4i}\sum_{k,l}\epsilon_{jkl}(2\xi_{ikl}-\xi_{kli})$
$\displaystyle=-iT^{\prime}_{ij},$ (S9)

and $T_{ij}$, $T^{\prime}_{ij}$ and $S_{ijk}$ are the notations we use in the
main text. Note that only the traceless part of the magneto-electric coupling
appears in the bulk conductivity tensor $\sigma_{ijk}$, indiciating that the
axion magneto-electric effect is missing in this approach.

## References

  * [1] Malashevich, A. & Souza, I.  Band theory of spatial dispersion in magnetoelectrics.  _Phys. Rev. B_ 82, 245118 (2010).
  * [2] Zhong, S., Moore, J. E. & Souza, I.  Gyrotropic magnetic effect and the magnetic moment on the Fermi surface.  _Phys. Rev. Lett._ 116, 077201 (2016).

Figure S1:  Layer dependence of the band structure. Figure S2:  Layer-
dependent spectrum of the electric susceptibility $\chi_{xx}$. Two peaks in
the spectrum correspond to the surface and bulk band gaps.