# Electrodynamics of Thin Sheets of Twisted Material

Dung Xuan Nguyen  Brown Theoretical Physics Center and Department of Physics,
Brown University, 182 Hope Street, Providence, Rhode Island 02912, USA  Dam
Thanh Son  Kadanoff Center for Theoretical Physics, University of Chicago, 933
East 56th Street, Illinois 60637, USA

###### Abstract

We construct a minimal theory describing the optical activity of a thin sheet
of a twisted material, the simplest example of which is twisted bilayer
graphene. We introduce the notion of “twisted electrical conductivity,” which
parametrizes the parity-odd response of a thin film to a perpendicularly
falling electromagnetic waves with wavelength larger than the thickness of the
sheet. We show that the low-frequency Faraday rotation angle has different
behaviors in different phases. For an insulator, the Faraday angle behaves as
$\omega^{2}$ at low frequencies, with the coefficient being determined by the
linear relationship between a component of the electric quadrupole moment and
the external electric field. For superconductors, the Faraday rotation angle
is constant when the frequency of the incoming EM waves is below the
superconducting gap and is determined by the coefficient of the Lifshitz
invariant in the Ginzburg-Landau functional describing the superconducting
state. In the metallic state, we show that the twisted conductivity is
proportional to the “magnetic helicity” (scalar product of the velocity and
the magnetic moment) of the quasiparticle, averaged around the Fermi surface.
The theory is general and is applicable to strongly correlated phases.

_Introduction._ —Recently, twisted quasi-two-dimensional materials have
attracted intense intention. The most well-known example of such materials is
twisted bilayer graphene, where flat bands have been predicted and found and
correlated insulating and superconducting states observed Bistritzer and
MacDonald (2011); Lopes dos Santos  _et al._ (2007); Suárez Morell  _et al._
(2010); Lopes dos Santos  _et al._ (2012); Cao  _et al._ (2018a, b). Flat
bands have also been discovered in the other twisted systems Carr  _et al._ ;
Kariyado and Vishwanath (2019), including twisted double bilayer graphene Shen
_et al._ (2020); Liu  _et al._ (2020); Cao  _et al._ (2020) and twisted
trilayer graphene Tsai  _et al._ .

One probe of the states of the twisted quasi-2D materials is through their
interaction with the external electromagnetic field. Due to the twisted
structure, a circularly polarized electromagnetic wave passing through such a
material can experience different responses depending on the sign of the
polarization—the so-called Faraday rotation and circular dichroism Kim  _et
al._ (2016); Addison  _et al._ (2019). In this paper we study the question:
what does the chiral response of a twisted quasi-2D material to long-
wavelength electromagnetic waves reveal about the structure of the low-energy
excitations in such a materials?

We first construct here a minimal theory that describes the chiral response of
a thin sheet of material with thickness $d$ much smaller than the wavelength
$\lambda$ of the incoming light. We show that in addition to the usual 2D
electrical conductivity, one just needs a new kinetic coefficient which we dub
“twisted conductivity,” which determines the electric current generated by a
gradient of the electric field, or the time dependence of the magnetic field.
The same quantity also determines the “dipole current,” to be defined, in an
external electric field.

In the metallic phase this twisted conductivity is shown to be proportional to
the “magnetic helicity” of the quasiparticle excitations on the Fermi surface,
defined as the product of the velocity and the in-plane magnetic moment of the
Landau quasiparticle.

Previous theoretical work include Refs. Suárez Morell  _et al._ (2017);
Stauber  _et al._ (2018); Ochoa and Asenjo-Garcia (2020), where one considers
a model of two layers, the current on one layer is linearly dependent on the
electric field on that layer and on the other layer. This approach is limited
to two-layer systems within the approximation that the electrons are localized
in a very thin shell around each layer, whose thickness is much less than the
distance between the layers. In Ref. Wang  _et al._ (2020) the Kubo formula
for the chiral response is derived for commensurate twist angles.

_General consideration._ —Consider a thin sheet of material, stretched along
the $(x,y)$ directions. We will have in mind twisted bilayer graphene as a
prototype. We are interested in the chiral response of this materials to
electromagnetic waves which have wavelength much larger than the thickness of
the sheet. For example, we want to understand Faraday’s rotation of the plane
of polarization of an incoming electromagnetic wave.

We will be mostly interested in the case when the incoming electromagnetic
waves fall perpendicularly to the plane of the sheet. In this case the is
electric field parallel to the sheet.

We assume that the sheet is invariant under a discrete group of rotations
around the $z$ axis, sufficiently large for the conductivity tensor to be
isotropic (for example $C_{3z}$ in the case of twisted graphene sheets). The
twisted sheet is assumed to be not invariant under 2D reflection but only
under a 2D reflection (say, $x\to x$, $y\to-y$) combined with $z\to-z$, which
is simply the two-fold rotation around the $x$ axis $C_{2x}$. This symmetry is
present in twisted bilayer graphene Balents (2019). We will assume that this
symmetry is not spontaneously broken. We will also assume time reversal
invariance.

In the standard treatment, the response of a single thin layer (say, of
graphene) to an incident electromagnetic wave is characterized by the 2D
(complex) conductivity $\sigma(\omega)$ Nair  _et al._ (2008); Mak  _et al._
(2008); Li and Heinz (2018). To treat effects like Faraday rotation or
circular dichroism, the surface conductivity is not sufficient. We now develop
a formalism, which deals exclusively with two-dimensional quantities, but
still allows one to captures these effects.

Let us denote the current in a finite-width piece of material by
${\mathbf{J}}(z,{\mathbf{x}})$. If the sheet is thin, we can define the 2D
current as

${\mathbf{j}}({\mathbf{x}})=\int\\!dz\,{\mathbf{J}}(z,{\mathbf{x}}),$ (1)

and the linear electric response is given by the frequency-dependent
electrical conductivity: ${\mathbf{j}}=\sigma{\mathbf{E}}$. The transmission
and reflection amplitudes of an electromagnetic wave falling perpendicularly
onto the sheet then can be expressed through the real and imaginary parts of
$\sigma$ Li and Heinz (2018). To go beyond $\sigma$, we introduce a new
current

$\bm{\zeta}({\mathbf{x}})=\int\\!dz\,z{\mathbf{J}}(z,{\mathbf{x}}).$ (2)

For a system which consists of two layers separated by a distance $d$, if the
thickness of the electron orbitals on each of the layers is much smaller than
$d$, then $\bm{\zeta}=\frac{d}{2}({\mathbf{j}}_{1}-{\mathbf{j}}_{2})$, where
${\mathbf{j}}_{1}$ and ${\mathbf{j}}_{2}$ are the electric currents on the
layers $1$ and $2$. We term this quantity the “dipole current,” because if
there is no tunneling between the two layers, then the electric dipole moment
along the $z$ direction is conserved, and $\bm{\zeta}$ is the current of that
conserved charge.

The source of the current ${\mathbf{j}}$ are the electric field
${\mathbf{E}}$, and the source of the dipole current $\bm{\zeta}$ is the
gradient of the electric field along the perpendicular direction
$\partial_{z}{\mathbf{E}}$. This can be seen by noticing that the Ohmic heat
generated in the system per unit area is

$\int\\!dz\,{\mathbf{J}}(z)\cdot{\mathbf{E}}(z)\approx\\!\int\\!dz\,{\mathbf{J}}(z)\cdot({\mathbf{E}}+z\partial_{z}{\mathbf{E}})={\mathbf{j}}\cdot{\mathbf{E}}+\bm{\zeta}\cdot\partial_{z}{\mathbf{E}}.$
(3)

We now introduce the conductivities $\sigma$, $\sigma_{1}$,
$\tilde{\sigma}_{1}$, and $\sigma_{2}$:

$\displaystyle{\mathbf{j}}$
$\displaystyle=\sigma{\mathbf{E}}-\sigma_{1}\partial_{z}{\mathbf{E}}\times\mathbf{\hat{z}},$
(4) $\displaystyle\bm{\zeta}$
$\displaystyle=\tilde{\sigma}_{1}({\mathbf{E}}\times\mathbf{\hat{z}})+\sigma_{2}\partial_{z}{\mathbf{E}}.$
(5)

In general the conductivities are functions of the frequency. Time reversal
invariance implies the Onsager relation
$\tilde{\sigma}_{1}(\omega)=\sigma_{1}(\omega)$. Positivity of entropy
production implies that both $\sigma$ and $\sigma_{2}$ have positive real
parts, and
$\mathop{\mathrm{Re}}\sigma\mathop{\mathrm{Re}}\sigma_{2}\geq(\mathop{\mathrm{Re}}\sigma_{1})^{2}$.

On dimensional ground, we expect that typically $\sigma_{1}\sim d\sigma$,
$\sigma_{2}\sim d^{2}\sigma$. For electromagnetic waves with wavelength
$\lambda$, the effects of $\sigma_{1}$ and $\sigma_{2}$ will be suppressed by
$d/\lambda$ and $(d/\lambda)^{2}$, respectively. However, since $\sigma_{d}$
is the first transport coefficient that breaks reflection symmetry, we need to
keep it in order to compute, e.g., Faraday rotation. Although the effect
$\sigma_{2}$ is always small, we will however keep it in our formalism for the
sake of symmetry.

Considering an electromagnetic plane wave with frequency $\omega$ falling
perpendicularly onto the plane. We can work in the gauge where
$A_{0}=A_{z}=0$, where the Maxwell equation for the perpendicular components
of the vector potential is (we use the Gauss units)

$-\partial_{z}^{2}\mathbf{A}-\frac{\omega^{2}}{c^{2}}\mathbf{A}=\frac{4\pi}{c}{\mathbf{J}}.$
(6)

(We assume here that the thin sheet is immerse in vacuum. Our calculation can
be trivially modified when the medium on one or both sides of the sheet has a
nontrivial dielectric constant.) To find the boundary conditions at $z=0$ we
can integrate over $z$ from $z=-\epsilon$ to $z=+\epsilon$

$-{\mathbf{A}}^{\prime}|^{\epsilon}_{-\epsilon}=\frac{4\pi}{c}{\mathbf{j}}=\frac{4\pi
i\omega}{c^{2}}(\sigma{\mathbf{A}}-\sigma_{1}{\mathbf{A}}^{\prime}\times\mathbf{\hat{z}}).$
(7)

We can also multiply Eq. (6) by $z$ and integrate over $z$ to find

$\delta{\mathbf{A}}=\frac{4\pi}{c}\bm{\zeta}=\frac{4\pi
i\omega}{c^{2}}(\sigma_{1}{\mathbf{A}}\times\mathbf{\hat{z}}+\sigma_{2}{\mathbf{A}}^{\prime}),$
(8)

where we have denoted the jump of a quantity $A$ across the sheet by $\delta
A\equiv A(+\epsilon)-A(-\epsilon)$. These boundary conditions can be written
for the two circular polarizations $A_{\pm}=A_{x}\pm iA_{y}$ separately,

$\displaystyle-A_{\pm}^{\prime}|^{\epsilon}_{-\epsilon}$
$\displaystyle=\frac{4\pi i\omega}{c^{2}}(\sigma A_{\pm}\pm
i\sigma_{1}A_{\pm}^{\prime}),$ (9) $\displaystyle
A_{\pm}|^{\epsilon}_{-\epsilon}$ $\displaystyle=\frac{4\pi i\omega}{c^{2}}(\mp
i\sigma_{1}A_{\pm}+\sigma_{2}A^{\prime}_{\pm}).$ (10)

We note here some peculiarities of these boundary conditions. In the usual
problem, only $\sigma$ is nonzero, and the boundary condition (9) specifies
that $A_{a}$ is continuous across $z=0$, while its first derivative
$A_{a}^{\prime}$ has a discontinuity proportional to the value of $A_{a}$ at
$z=0$. This leads to a consistent mathematical problem. But with nonzero
$\sigma_{1}$ and $\sigma_{2}$, Eqs. (9) and (10) imply that both $A_{a}$ and
$A^{\prime}$ are discontinuous, and their discontinuity are proportional to
their values at $z=0$. The boundary conditions become ambiguous: at which
value of $z$ should the right-hand sides of Eqs. (9) and (10) be calculated?

To fully resolve the ambiguity, one needs to solve the scattering problem for
a slab of finite width and carefully take the limit of small width. This is
done in the Supplementary Materials. We will use a shortcut that leads to the
correct answer when the jumps of $A$ and $A^{\prime}$ are small compared to
their values. This requires

$\frac{4\pi}{c^{2}}\sigma_{1}\omega\ll 1.$ (11)

This condition is usually satisfied. For example for a sheet a few atomic
layer thick, if $\sigma\sim e^{2}/h$ then $\sigma_{1}\sim\frac{e^{2}}{h}d$.
Then

$\frac{4\pi\omega}{c^{2}}\sigma_{1}\sim\frac{e^{2}}{\hbar
c}\frac{d}{\lambda}\ll 1.$ (12)

When the jump of $A$ is much smaller than $A$, we can replace $A$ on the
right-hand side of the boundary condition by the average values on the two
sides. This prescription is similar to Griffiths’s treatment of the so-called
$\delta^{\prime}(x)$ potential in quantum mechanics Griffiths (1993). The
boundary conditions can now be written as

$\displaystyle-\delta A_{\pm}^{\prime}$ $\displaystyle=\frac{4\pi
i\omega}{c^{2}}(\sigma\bar{A}_{\pm}\pm i\sigma_{1}\bar{A}_{\pm}^{\prime}),$
(13a) $\displaystyle\delta A_{\pm}$ $\displaystyle=\frac{4\pi
i\omega}{c^{2}}(\mp
i\sigma_{1}\bar{A}_{\pm}+\sigma_{2}\bar{A}^{\prime}_{\pm}).$ (13b)

where we have denoted $\bar{A}\equiv\frac{1}{2}[A(+\epsilon)+A(-\epsilon)]$.

Let us now consider a problem of scattering of a plane wave onto the plane. A
plane wave coming from $z=-\infty$ gives rise to a transmitted and a scattered
waves.

$A_{\pm}(z)=\begin{cases}e^{ikz}+R_{\pm}e^{-ikz},&z<0,\\\
T_{\pm}e^{ikz},&z>0.\end{cases}$ (14)

With the boundary conditions, we find

$\displaystyle T_{\pm}$
$\displaystyle=\frac{1\pm\tilde{\sigma}_{1}k+\frac{1}{4}(\tilde{\sigma}_{1}^{2}-\tilde{\sigma}\tilde{\sigma}_{2})k^{2}}{1+\frac{1}{2}\tilde{\sigma}+\frac{1}{4}(2\tilde{\sigma}_{2}-\tilde{\sigma}_{1}^{2}+\tilde{\sigma}\tilde{\sigma}_{2})k^{2}}\,,$
(15a) $\displaystyle R_{\pm}$
$\displaystyle=-\frac{1}{2}\frac{\tilde{\sigma}-\tilde{\sigma}_{2}k^{2}}{1+\frac{1}{2}\tilde{\sigma}+\frac{1}{4}(2\tilde{\sigma}_{2}-\tilde{\sigma}_{1}^{2}+\tilde{\sigma}\tilde{\sigma}_{2})k^{2}}\,,$
(15b)

where

$\tilde{\sigma}=\frac{4\pi}{c}\sigma,\quad\tilde{\sigma}_{1}=\frac{4\pi}{c}\sigma_{1},\quad\tilde{\sigma}_{2}=\frac{4\pi}{c}\sigma_{2}.$
(16)

As a consistency check of the symmetrized prescription (13), one can verify
that if
$\mathop{\mathrm{Re}}\sigma=\mathop{\mathrm{Re}}\sigma_{1}=\mathop{\mathrm{Re}}\sigma_{2}=0$,
implying no dissipation, then $|T|^{2}+|R|^{2}=1$.

Ignoring $\sigma_{2}$ and assuming $\tilde{\sigma}_{1}k\ll 1$, we have

$T_{\pm}=\frac{{1\pm\frac{4\pi}{c}k\sigma_{1}}}{1+\frac{2\pi\sigma}{c}}\,,\quad
R_{\pm}=-\frac{\frac{2\pi}{c}\sigma}{1+\frac{2\pi}{c}\sigma}\,.$ (17)

The Faraday rotation angle is then

$\theta_{\rm
F}=\frac{1}{2}\arg\frac{T_{+}}{T_{-}}=\frac{4\pi\omega}{c^{2}}\mathop{\mathrm{Im}}\sigma_{1}.$
(18)

while the ratio of the absolute value of the transmission amplitude of the two
polarizations (which is different from 1 if there is circular dichroism) is

$\frac{|T_{+}|^{2}}{|T_{-}|^{2}}=1+\frac{16\pi\omega}{c^{2}}\mathop{\mathrm{Re}}\sigma_{1}.$
(19)

Both Eqs. (18) and (19) do not require $\sigma$ to be small to be valid, but
only (11).

In the Supplementary Materials we reproduce the formulas by solving the
Maxwell equations for a slab of finite thickness in the limit of long
wavelength. Equations (18) and (19), as written, are correct also in the
presence of a dielectric constant on one or or both sides of the sheet.

_Metallic state._ —Consider a state with a Fermi surface. In the regime where
the incoming photon has energy less than the Fermi energy (typically that
means wavelength in the infrared range) we can use the Fermi liquid theory to
treat the problem. In this description, we have quasiparticles states
$|{\mathbf{p}},\alpha\rangle$ where ${\mathbf{p}}$ is a momentum along the
sheet and $\alpha$ is an internal index (for twisted bilayer graphene $\alpha$
would corresponds to the spin, valley, and layer degeneracies). The state has
a magnetic moment $\bm{\mu}^{\alpha}_{\mathbf{p}}$, parallel to the sheet. One
can visualize this magnetic moment as arising from the spiral-like motion of a
wave packet when it moves along the sheet, jumping back and forth between the
two layers. The energy in a magnetic field ${\mathbf{B}}$ has the form (from
now on we suppress the $\alpha$ index)

$E_{\mathbf{p}}=\varepsilon_{\mathbf{p}}-\bm{\mu}_{\mathbf{p}}\cdot{\mathbf{B}}.$
(20)

Note that for fields that are constant on the $(x,y)$ plane,
${\mathbf{B}}=-\partial_{z}{\mathbf{A}}\times\mathbf{\hat{z}}$. Ignoring
Fermi-liquid effects, the kinetic equation for the distribution function
$f_{\mathbf{p}}$ in the relaxation-time approximation is

$\frac{\partial f_{\mathbf{p}}}{\partial t}-e{\mathbf{E}}\cdot\frac{\partial
f_{\mathbf{p}}}{\partial{\mathbf{p}}}=-\frac{f_{\mathbf{p}}-f_{0}(E_{\mathbf{p}})}{\tau}\,,$
(21)

where $f_{0}(\varepsilon)=[e^{\beta(\varepsilon-\mu)}+1]^{-1}$. Linearizing
the equation: $f_{\mathbf{p}}=f_{0}(\varepsilon_{\mathbf{p}})+\delta
f_{\mathbf{p}}$, we find

$\delta
f_{\mathbf{p}}=\frac{1}{1-i\omega\tau}\left(-e\tau{\mathbf{E}}\cdot{\bm{v}}_{\mathbf{p}}+\bm{\mu}_{\mathbf{p}}\cdot{\mathbf{B}}\right)\left(-\frac{\partial
f_{0}}{\partial\varepsilon}\right).$ (22)

To compute the current, we use

$\mathbf{j}=-e\\!\int_{\mathbf{p}}\frac{\partial
E_{\mathbf{p}}}{\partial{\mathbf{p}}}(f_{0}+\delta f_{\mathbf{p}}),$ (23)

where $\int_{\mathbf{p}}\equiv\int d^{2}{\mathbf{p}}/(2\pi)^{2}$, and the
group velocity is computed using the full dispersion (20). Inserting Eq. (22)
we find

$\mathbf{j}=\frac{e\tau}{1-i\omega\tau}\\!\int_{\mathbf{p}}[e({\mathbf{E}}\cdot{\bm{v}}_{\mathbf{p}})-i\omega(\bm{\mu}_{\mathbf{p}}\cdot{\mathbf{B}})]{\bm{v}}_{\mathbf{p}}\left(-\frac{\partial
f_{0}}{\partial\epsilon}\right).$ (24)

By using the Maxwell equation
$\bm{\nabla}\times{\mathbf{E}}=-c^{-1}\dot{\mathbf{B}}$, one gets

$\sigma=\frac{1}{2}\frac{e^{2}\tau\nu(\varepsilon_{\rm
F})}{1-i\omega\tau}\langle{\bm{v}}_{\mathbf{p}}^{2}\rangle,\quad\sigma_{1}=-\frac{1}{2}\frac{ec\tau\nu(\varepsilon_{\rm
F})}{1-i\omega\tau}\langle{\bm{v}}_{\mathbf{p}}\cdot\bm{\mu}_{\mathbf{p}}\rangle.$
(25)

Here $\nu(\varepsilon_{\rm F})$ is the density of state at the Fermi level and
$\langle\cdot\rangle$ means averaging over the Fermi surface: $\langle
A\rangle\equiv\nu^{-1}(\varepsilon_{\rm
F})\int_{\mathbf{p}}\\!A\delta(\varepsilon_{\mathbf{p}}-\varepsilon_{\rm F})$.

We will call the scalar product of the velocity ${\bm{v}}_{\mathbf{p}}$ and
the in-plane orbital magnetic moment $\bm{\mu}_{\mathbf{p}}$ the “magnetic
helicity” of the quasiparticle carrying momentum ${\mathbf{p}}$. The average
of the magnetic helicity around the Fermi surface
$\langle{\bm{v}}_{\mathbf{p}}\cdot\bm{\mu}_{\mathbf{p}}\rangle$ encodes the
correlation between between the direction of motion of the quasiparticles and
their magnetic moment. Such correlation leads to the appearance of a magnetic
moment density when an electric current flows through the system. Equation
(25) implies that the twisted conductivity is related to this locking between
velocity and magnetic moment. This locking is also the origin of the
gyrotropic magnetic effect in 3D Ma and Pesin (2015); Zhong  _et al._ (2016)
and 2D Wang  _et al._ (2020).

Note that the ratio of $\sigma_{1}$ and $\sigma$ does not depend on the mean
free time $\tau$,

$\frac{\sigma_{1}}{\sigma}=-\frac{c}{e}\frac{\langle{\bm{v}}_{p}\cdot\bm{\mu}_{\mathbf{p}}\rangle}{\langle{\bm{v}}_{\mathbf{p}}^{2}\rangle}\,,$
(26)

and can be read out just from the wave functions of modes at the Fermi
surface. This has the consequence that if $\sigma\ll\frac{c}{2\pi}\approx
137\frac{e^{2}}{h}$, circular dichroism, defined as
$\text{CD}=(\mathcal{A}_{+}-\mathcal{A}_{-})/[2(\mathcal{A}_{+}+\mathcal{A}_{-})]$
wherer $\mathcal{A}_{+}$ and $\mathcal{A}_{-}$ are the absorption coefficient
of the corresponding helicities, is linear in frequency

$\text{CD}=\frac{\omega}{e}\frac{\langle{\bm{v}}_{p}\cdot\bm{\mu}_{\mathbf{p}}\rangle}{\langle{\bm{v}}_{\mathbf{p}}^{2}\rangle}\,.$
(27)

In the regime $\omega\tau\ll 1$, $\mathop{\mathrm{Im}}\sigma_{1}\sim\omega$,
so the Faraday rotation angle is of order $\omega^{2}$. In the opposite regime
$\omega\tau\gg 1$, the conductivities become purely imaginary and inversely
proportional to the frequency $\omega$. In particular, the twisted
conductivity $\mathop{\mathrm{Im}}\sigma_{1}\sim\omega^{-1}$. This means that
the Faraday rotation angle is constant in this regime

$\theta_{\rm F}=-\frac{2\pi e}{c}\nu(\varepsilon_{\rm
F})\langle{\bm{v}}_{\mathbf{p}}\cdot\bm{\mu}_{\mathbf{p}}\rangle.$ (28)

To illustrate the result, consider the case of twisted bilayer graphene. Here
$\nu(\varepsilon_{\rm F})$ contains a factor of 8 from the valley, spin, and
layer degeneracies. For large twisting angle $\theta$ (larger than the angle
1.1∘ where flat bands appear) and small Fermi momentum $p_{\rm F}\ll
k_{\theta}=2k_{\rm D}\sin\frac{\theta}{2}$ where $k_{\rm D}$ is the Dirac
momentum, the quasiparticles are localized on one layer, with small admixture
from the other layer. The magnetic moment of the quasiparticle is mostly
perpendicular to the momentum, but there is a small component along the
quasiparticle’s momentum, leading to a nonzero average magnetic helicity. The
angle of Faraday rotation in this case is (see the Supplementary Materials)

$\theta_{\rm
F}=-24\alpha\frac{v_{0}}{c}\frac{k_{\theta}d}{\hbar}\left(\frac{p_{\rm
F}}{k_{\theta}}\right)^{2}\frac{w_{\rm AA}w_{\rm
AB}}{(v_{0}k_{\theta})^{2}}\,,$ (29)

where $\alpha\approx\frac{1}{137}$ is the fine structure constant, $v_{0}$ is
the velocity of the Dirac fermion in single-layer graphene, $w_{\rm AA}$ and
$w_{\rm AB}$ are the interlayer coupling parameters, and $d$ is the distance
between the layers. For example, at twist angle $\theta=2^{\circ}$, using the
standard numerical values for the parameters (see, e.g., Ref. Tarnopolsky  _et
al._ (2019)) one gets $\theta_{\rm F}\sim 10^{-5}$ when $p_{F}$ is of the same
order of magnitude as $k_{\theta}$. The sign of $\theta_{\rm F}$ is such that
the plane of polarization is rotated in the same direction as the direction,
with respect to which the graphene layer farther from the source is rotated
with respect to the one closer to the source. Analogously, we obtain
$\text{CD}\sim 10^{-4}$ for $\omega\sim v_{0}k_{\theta}$.

_Superconducting case._ —We now consider the case of a superconducting thin
layer. Even when we do not know the microscopic mechanism for
superconductivity in twisted bilayer graphene and other twisted materials, we
can treat the problem of scattering of long-wavelength electromagnetic waves
(with frequency less than the superconducting gap, thus typically in the
microwave range) phenomenologically using the Ginzburg-Landau theory. Keeping
only the phase $\varphi$ of the order parameter, the Ginzburg-Landau energy
functional has the form Mineev and Samokhin (1994); Edelstein (1996)

$H=\\!\int\\!d^{3}x\,\delta(z)\biggl{[}\frac{\hbar^{2}n_{s}}{2m}(D_{a}\varphi)^{2}+\kappa
n_{s}B_{a}D_{a}\varphi\biggr{]},$ (30)

where the covariant derivatives are
$D_{t}\varphi=\partial_{t}\varphi+\frac{2e}{\hbar}A_{0}$,
$D_{a}\varphi=\partial_{a}\varphi+\frac{2e}{\hbar c}A_{a}$, and the index $a$
runs $x,y$. The term proportional to $\kappa$ in Eq. (30) is the so-called
“Lifshitz invariant” term. One can imagine that this term arises from a
nonminimal coupling of the order parameter $\psi$ with the electromagnetic
field: $\frac{i}{2}\kappa B_{a}(\psi^{\dagger}D_{a}\psi-
D_{a}\psi^{\dagger}\psi)$. This type of coupling has been considered
previously in the treatment of non-centrosymmetric superconductors Mineev and
Samokhin (1994); Edelstein (1996). In that context, the term is associated
with spin-orbit coupling.

For twisted bilayer graphene and other twisted materials, from symmetry
arguments one expects the Lifshitz invariant to appear independent of the
mechanism of superconductivity. Physically, the Lifshitz invariant term
implies that a moving condensate possesses a nonzero density of magnetic
moment. This is expected if the superconductivity is formed by Cooper pairing
of quasiparticles with nonzero average magnetic helicity, which, as we recall,
parametrizes the locking between the velocity and the magnetic moment.

One can couple this action with the electromagnetic field in the bulk and
solve the combined system of equations for the gauge field and the phase
$\varphi$. Equivalently, one can simply compute the twist conductivity from
Eq. (30). Differentiating the action with respect to $A_{a}$, and setting all
fields to be spatially homogeneous along the directions of the sheet, in
particular $\partial_{a}\varphi=0$ (which can be done when the electromagnetic
wave falls perpendicularly onto the sheet), one obtains a generalized London
equation

$J^{a}=-c\frac{\delta H}{\delta
A_{a}}=\delta(z)\left[-\frac{(2e)^{2}n_{s}}{mc}A_{a}-2e\frac{\kappa
n_{s}}{\hbar}B_{a}\right]\\\ +2e\frac{\kappa
n_{s}}{\hbar}\epsilon^{ab}\partial_{z}[\delta(z)A_{b}],$ (31)

from which one finds

$\displaystyle\sigma=\frac{i(2e)^{2}n_{s}}{m\omega}\,,\quad\sigma_{1}=\frac{ic}{\omega}2e\frac{\kappa
n_{s}}{\hbar}\,.$ (32)

Since $\sigma_{1}\sim\omega^{-1}$, the Faraday rotation angle is frequency-
independent. This can be seen by applying simple field-theoretical power
counting to the Lifshitz invariant term: as $B_{a}$ has dimension 2 and
$D_{a}\varphi$ has dimension one, $\kappa n_{s}$ must be dimensionless. Since
the Faraday rotation angle is dimensionless and proportional to $\kappa
n_{s}$, there should be no frequency dependence. This is exactly the same
behavior as for metals in the regime $\omega\tau\gg 1$.

The uniqueness of the Lifshitz invariant allows one to compute also the change
of the plane of polarization when the incoming light falls onto the sheet at
any angle. This calculation is done in the Supplementary Materials.

The Lifshitz invariant also leads to another effect—the chiral magnetic
Josephson effect in which the Josephson junction is built up from two chiral
superconductors linked by a uniaxial ferromagnet. Josephson current appears
even with zero phase difference dues to a phase offset Chernodub  _et al._ ;
Buzdin (2008). This phase offset originates from the parity breaking term and
is proportional to $\kappa$. This effect is similar to the chiral magnetic
effect in hot QCD Kharzeev (2006).

_Insulating states_.—In the insulating state, there is no low-energy degree of
freedom living on the thin layer, so the effective action is just a local
function of the electromagnetic field and its derivatives. Instead of the
Lifshitz invariant, the leading term which breaks the 2D reflection symmetry,
but preserve the combination of 2D reflection and $z\to-z$ symmetry, is

$\delta
S=\int\\!dt\,d^{3}x\,\delta(z)\beta\epsilon^{ab}E_{a}\partial_{z}E_{b},$ (33)

where $\beta$ is some constant. The coefficient $\beta$ has dimension $-2$,
thus one concludes that the Faraday rotation angle behaves like
$\beta\omega^{2}$ for waves with frequency less than the gap. This prediction
is independent of the nature of the insulating state, whether one is dealing
with a band insulator or a strongly correlated insulator. Calculations similar
to the one we have done for the superconducting case give $\theta_{\rm
F}=4\pi\beta(\omega/c)^{2}$.

To give a physical interpretation to the term (33) in the effective action,
let us imagine immersing a finite piece of thin twisted insulator into an
uniform static parallel electric field $E_{a}$. For simplicity let us also
assume that the boundary is also gapped. The effective action of this finite
piece of material is described by the same Eq. (33) with $\beta$ replaced by
the function $\beta(x)$, where $\beta(x)=\beta$ inside the piece and 0
outside. The charge density induced by a constant electric field field is

$\rho=\frac{\delta S}{\delta
A_{0}}=\delta^{\prime}(z)\epsilon^{ab}E_{a}\partial_{b}\beta(x).$ (34)

This charge distribution corresponds to zero total charge and dipole moment,
but leads to nonzero value of the electric quadrupole moment

$Q^{az}=\\!\int\\!dz\,d^{2}x\,zx^{a}\rho(z,x)=-\beta S\epsilon^{ab}E_{b},$
(35)

where $S$ is the total area of the piece of material. (This quadrupole moment
is well defined as it is along directions perpendicular to the direction of
the dipole moment.) Thus $\beta$ is the coefficient determining the quadrupole
polarization induced by an electric field, which provides an independent way
of measuring $\beta$.

_Conclusion._ —We have constructed the minimal theory describing the chiral
response of a thin sheet of a twisted material, applicable to twisted bilayer
graphene and related materials. In addition to the 2D electrical conductivity
of the sheet, one needs to introduce another transport coefficient, which we
term the “twisted conductivity.” For a metal, the twisted conductivity is
proportional to the average magnetic helicity (i.e., scalar product of the
velocity and the in-plane orbital magnetic moment) of the quasiparticles
around the Fermi line. For a superconductor, the chiral response is shown to
be related to the Lifshitz invariant, which parametrizes the magnetic moment
density created by a supercurrent. For the insulating phase, the chiral
response is related to the quadrupole moment induced by an in-plane electric
field.

The discussion in this paper is general and insensitive to the details of
microscopic physics. It would be interested to see if the chiral response of
the correlated phases is sensitive to the mechanism underlying these phases.

_Acknowledgments._ —The authors thank Stephen Carr, Van-Nam Do, Joel Moore,
Ashvin Vishwanath, and Grigory Tarnopolsky for discussions and comments on an
earlier version of this manuscript. DTS is supported, in part, by the by the
U.S. DOE grant No. DE-FG02-13ER41958, a Simons Investigator grant and by the
Simons Collaboration on Ultra-Quantum Matter from the Simons Foundation. DXN
was supported by Brown Theoretical Physics Center.

## References

  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proc. Natl. Acad. Sci. 108, 12233 (2011).
  * Lopes dos Santos  _et al._ (2007) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Graphene Bilayer with a Twist: Electronic Structure, Phys. Rev. Lett. 99, 256802 (2007).
  * Suárez Morell  _et al._ (2010) E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Flat bands in slightly twisted bilayer graphene: Tight-binding calculations, Phys. Rev. B 82, 121407(R) (2010).
  * Lopes dos Santos  _et al._ (2012) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Continuum model of the twisted graphene bilayer, Phys. Rev. B 86, 155449 (2012).
  * Cao  _et al._ (2018a) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018a).
  * Cao  _et al._ (2018b) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018b).
  * (7) S. Carr, D. Massatt, M. Luskin, and E. Kaxiras, Duality between atomic configurations and Bloch states in twistronic materials, arXiv:1803.01242 .
  * Kariyado and Vishwanath (2019) T. Kariyado and A. Vishwanath, Flat band in twisted bilayer Bravais lattices, Phys. Rev. Res. 1, 033076 (2019).
  * Shen  _et al._ (2020) C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, K. Watanabe, T. Taniguchi, R. Yang, Z. Y. Meng, D. Shi, O. V. Yazyev, and G. Zhang, Correlated states in twisted double bilayer graphene, Nature Physics 16, 520 (2020).
  * Liu  _et al._ (2020) X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, Y. Ronen, H. Yoo, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Tunable spin-polarized correlated states in twisted double bilayer graphene, Nature 583, 221 (2020).
  * Cao  _et al._ (2020) Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable correlated states and spin-polarized phases in twisted bilayer–bilayer graphene, Nature 583, 215 (2020).
  * (12) K.-T. Tsai, X. Zhang, Z. Zhu, Y. Luo, S. Carr, M. Luskin, E. Kaxiras, and K. Wang, Correlated Superconducting and Insulating States in Twisted Trilayer Graphene Moiré of Moiré Superlattices, arXiv:1912.03375 .
  * Kim  _et al._ (2016) C.-J. Kim, A. Sánchez-Castillo, Z. Ziegler, Y. Ogawa, C. Noguez, and J. Park, Chiral atomically thin films, Nature Nanotechnol. 11, 520 (2016).
  * Addison  _et al._ (2019) Z. Addison, J. Park, and E. J. Mele, Twist, slip, and circular dichroism in bilayer graphene, Phys. Rev. B 100, 125418 (2019).
  * Suárez Morell  _et al._ (2017) E. Suárez Morell, L. Chico, and L. Brey, Twisting dirac fermions: circular dichroism in bilayer graphene, 2D Materials 4, 035015 (2017).
  * Stauber  _et al._ (2018) T. Stauber, T. Low, and G. Gómez-Santos, Chiral Response of Twisted Bilayer Graphene, Phys. Rev. Lett. 120, 046801 (2018).
  * Ochoa and Asenjo-Garcia (2020) H. Ochoa and A. Asenjo-Garcia, Flat Bands and Chiral Optical Response of Moiré Insulators, Phys. Rev. Lett. 125, 037402 (2020).
  * Wang  _et al._ (2020) Y.-Q. Wang, T. Morimoto, and J. E. Moore, Optical rotation in thin chiral/twisted materials and the gyrotropic magnetic effect, Phys. Rev. B 101, 174419 (2020).
  * Balents (2019) L. Balents, General continuum model for twisted bilayer graphene and arbitrary smooth deformations, SciPost Phys. 7, 48 (2019).
  * Nair  _et al._ (2008) R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim, Fine Structure Constant Defines Visual Transparency of Graphene, Science 320, 1308 (2008).
  * Mak  _et al._ (2008) K. F. Mak, M. Y. Sfeir, Y. Wu, C. H. Lui, J. A. Misewich, and T. F. Heinz, Measurement of the Optical Conductivity of Graphene, Phys. Rev. Lett. 101, 196405 (2008).
  * Li and Heinz (2018) Y. Li and T. F. Heinz, Two-dimensional models for the optical response of thin films, 2D Materials 5, 025021 (2018).
  * Griffiths (1993) D. J. Griffiths, Boundary conditions at the derivative of a delta function, J. Phys. A 26, 2265 (1993).
  * Ma and Pesin (2015) J. Ma and D. A. Pesin, Chiral magnetic effect and natural optical activity in metals with or without Weyl points, Phys. Rev. B 92, 235205 (2015).
  * Zhong  _et al._ (2016) S. Zhong, J. E. Moore, and I. Souza, Gyrotropic Magnetic Effect and the Magnetic Moment on the Fermi Surface, Phys. Rev. Lett. 116, 077201 (2016).
  * Tarnopolsky  _et al._ (2019) G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Origin of Magic Angles in Twisted Bilayer Graphene, Phys. Rev. Lett. 122, 106405 (2019).
  * Mineev and Samokhin (1994) V. P. Mineev and K. V. Samokhin, Helical phases in superconductors, Sov. Phys. JETP 78, 747 (1994).
  * Edelstein (1996) V. M. Edelstein, The Ginzburg–Landau equation for superconductors of polar symmetry, J. Phys.: Condens. Matter 8, 339 (1996).
  * (29) M. N. Chernodub, J. Garaud, and D. E. Kharzeev, Chiral Magnetic Josephson junction: a base for low-noise superconducting qubits?, arXiv:1908.00392 .
  * Buzdin (2008) A. Buzdin, Direct Coupling Between Magnetism and Superconducting Current in the Josephson ${\varphi}_{0}$ Junction, Phys. Rev. Lett. 101, 107005 (2008).
  * Kharzeev (2006) D. Kharzeev, Parity violation in hot QCD: Why it can happen, and how to look for it, Phys. Lett. B 633, 260 (2006).

— Supplementary Material —
Electrodynamics of Thin Sheets of Twisted Material
Dung Xuan Nguyen and Dam Thanh Son

##  I From 3D to 2D

In this Section, we consider the problem of the scattering of electromagnetic
waves falling perpendicularly onto a slab of a finite width. Then taking the
limit where the thickness of the slab is much smaller than the wavelength of
the incoming radiation, we reproduce the formulas for the Faraday rotation
angle and circular dichroism obtained in the main text.

Consider a layer of material, whose response to an external electric field is
given by a nonlocal conductivity $\sigma(\omega;z,z^{\prime})$

$j_{a}(\omega;z)=\int\\!dz^{\prime}\,\sigma_{ab}(\omega;z,z^{\prime})E_{b}(\omega,z^{\prime}).$
(S1)

We assume that the layer has finite thickness $d$, so $\sigma_{ab}=0$ when
$|z|>d/2$ or $|z^{\prime}|>d/2$, or both. We will omit the argument $\omega$
in further formulas. We want to derive the formula for the transmission and
reflection coefficient in the limit of small $k$, $kd\ll 1$.

First, we can decompose $\sigma_{ab}$ into symmetric and antisymmetric parts

$\sigma_{ab}=\sigma_{S}\delta_{ab}+\sigma_{A}\epsilon_{ab}.$ (S2)

In the helicity basis

$j_{\pm}=j_{x}\pm ij_{y},\qquad E_{\pm}=E_{x}\pm iE_{y},$ (S3)

Eq. (S1) can be written as

$j_{\pm}(z)=\int\\!dz^{\prime}\,\sigma_{\pm}(z,z^{\prime})E_{\pm}(z^{\prime}),$
(S4)

where we defined

$\sigma_{\pm}=\sigma_{S}\mp i\sigma_{A}.$ (S5)

The symmetry with respect to reflection around the $y$ axis combined with the
exchange of the upper and lower layers $z\to-z$ (which is the two-fold
rotation around the $x$ axis) implies

$\sigma_{+}(z,z^{\prime})=\sigma_{-}(-z,-z^{\prime}).$ (S6)

Time reversal invariance implies
$\sigma_{ab}(\omega;z,z^{\prime})=\sigma_{ba}(\omega;z^{\prime},z)$, which
means

$\sigma_{+}(z,z^{\prime})=\sigma_{-}(z^{\prime},z).$ (S7)

Combining two symmetries, we find

$\sigma_{\pm}(z,z^{\prime})=\sigma_{\pm}(-z^{\prime},-z).$ (S8)

From now on we focus on the positive helicity, and drop the helicity index $+$
in formulas. Let us also define

$\sigma(z)=\int\\!dz^{\prime}\,\sigma(z,z^{\prime}).$ (S9)

From Eq. (S8) we find

$\int\\!dz\,\sigma(z,z^{\prime})=\sigma(-z^{\prime}).$ (S10)

The Maxwell equation can be written as

$\partial_{z}^{2}A(z)+k^{2}A(z)=-ik\\!\int\\!dz^{\prime}\,\tilde{\sigma}(z,z^{\prime})A(z^{\prime}),\quad\tilde{\sigma}(z,z^{\prime})=\frac{4\pi}{c}\sigma(z,z^{\prime}).$
(S11)

We will be looking for solution to this equation with the asymptotics

$A(z)\to e^{ikz},\quad z\to+\infty,$ (S12)

which means that in the opposite limit the asymptotics of $A(z)$ is

$A(z)=\frac{1}{T}e^{ikz}+\frac{R}{T}e^{-ikz},$ (S13)

where $T$ and $R$ are the transmission and reflection amplitudes,
respectively. Using Green’s function, Eq. (S11) can be written as

$A(z)=e^{ikz}-i\\!\int\\!dz_{1}\,dz_{2}\,\theta(z_{1}-z)\sin
k(z_{1}-z)\tilde{\sigma}(z_{1},z_{2})A(z_{2}),$ (S14)

from which we find

$\displaystyle\frac{1}{T}$
$\displaystyle=1+\frac{1}{2}\\!\int\\!dz_{1}\,dz_{2}\,e^{-ikz_{1}}\tilde{\sigma}(z_{1},z_{2})A(z_{2}),$
(S15a) $\displaystyle\frac{R}{T}$
$\displaystyle=-\frac{1}{2}\\!\int\\!dz_{1}\,dz_{2}\,e^{ikz_{1}}\tilde{\sigma}(z_{1},z_{2})A(z_{2}).$
(S15b)

The Maxwell equation (S11) can be solved by iteration, where the solution is
presented as an infinite series

$A(z)=\sum_{n=0}^{\infty}A^{(n)}(z),$ (S16)

where we define

$\displaystyle A^{(0)}(z)$ $\displaystyle=e^{ikz},$ (S17) $\displaystyle
A^{(n)}(z)$ $\displaystyle=-i\\!\int\\!dz_{1}\,dz_{2}\,\theta(z_{1}-z)\sin
k(z_{1}-z)\tilde{\sigma}(z_{1},z_{2})A^{(n-1)}(z_{2}).$ (S18)

In the limit of small $k$ the expansion (S16) is an expansion over $k$. First
consider $z$ inside the slab, $|z|<d/2$, then $\sin k(z-z_{1})\sim kd$ since
$z_{1}$ has to be inside the slab, we have

$A^{(n)}(z)\sim kA^{(n-1)}(z)\Rightarrow
A^{(n)}(z)\sim(kd)^{n},\qquad|z|<\frac{d}{2}\,.$ (S19)

Let us first compute $T$ and $R$ to order $k^{0}$. We replace in Eq. (S15)
$A(z)\to A^{(0)}(z)=e^{ikz}$. To order $k^{0}$ we can replace $e^{\pm
ikz_{1}}$ and $e^{ikz_{2}}$ by $1$ since $z_{1},z_{2}\lesssim d$. We find

$\displaystyle\frac{1}{T}$
$\displaystyle=1+\frac{1}{2}\\!\int\\!dz_{1}\,dz_{2}\,\tilde{\sigma}(z_{1},z_{2}),$
(S20) $\displaystyle\frac{R}{T}$
$\displaystyle=-\frac{1}{2}\int\\!dz_{1}\,dz_{2}\,\tilde{\sigma}(z_{1},z_{2}).$
(S21)

Comparing with Eqs. (17) we can identify

$\sigma=\int\\!dz_{1}\,dz_{2}\,\sigma(z_{1},z_{2}).$ (S22)

To compute $T$ and $R$ to order $k^{1}$ we need to replace, in Eqs. (S15)
$A(z)$ by $A^{(0)}(z)+A^{(1)}(z)$, where $A^{(0)}$ is computed to next-to-
leading order in $k$, and $A^{(1)}$ to leading order,

$\displaystyle A^{(0)}(z)$ $\displaystyle=1+ikz+O(k^{2}),$ (S23)
$\displaystyle A^{(1)}(z)$
$\displaystyle=-ik\\!\int\\!dz_{1}\,\theta(z_{1}-z)(z_{1}-z)\tilde{\sigma}(z_{1})+O(k^{2}).$
(S24)

We find

$\displaystyle\frac{1}{T}$
$\displaystyle=1+\frac{1}{2}\int\\!dz\,\tilde{\sigma}(z)-ik\\!\int\\!dz\,z\tilde{\sigma}(z)-\frac{ik}{2}\\!\int\\!dz\,dz^{\prime}\,\theta(z+z^{\prime})(z+z^{\prime})\tilde{\sigma}(z)\tilde{\sigma}(z^{\prime}),$
(S25) $\displaystyle\frac{R}{T}$
$\displaystyle=-\frac{1}{2}\int\\!dz\,\tilde{\sigma}(z)+\frac{ik}{2}\\!\int\\!dz\,dz^{\prime}\,\theta(z+z^{\prime})(z+z^{\prime})\tilde{\sigma}(z)\tilde{\sigma}(z^{\prime}).$
(S26)

Let us now divide $\tilde{\sigma}(z)$ into symmetric and antisymmetric parts

$\sigma(z)=\sigma_{s}(z)+\sigma_{a}(z),\qquad\sigma_{s}(z)=\sigma_{s}(-z),\quad\sigma_{a}(z)=-\sigma_{a}(-z)$
(S27)

Then we find

$\displaystyle\frac{1}{T}$
$\displaystyle=\left[1+\frac{1}{2}\int\\!dz\,\tilde{\sigma}_{s}(z)\right]\left[1-ik\\!\int\\!dz\,z\tilde{\sigma}_{a}(z)\right]-\frac{ik}{2}\\!\int\\!dz\,dz^{\prime}\,\theta(z+z^{\prime})(z+z^{\prime})[\tilde{\sigma}_{s}(z)\tilde{\sigma}_{s}(z^{\prime})+\tilde{\sigma}_{a}(z)\tilde{\sigma}_{a}(z^{\prime})],$
(S28) $\displaystyle\frac{R}{T}$
$\displaystyle=-\frac{1}{2}\int\\!dz\,\tilde{\sigma}(z)\left[1-ik\\!\int\\!dz\,z\tilde{\sigma}_{a}(z)\right]+\frac{ik}{2}\\!\int\\!dz\,dz^{\prime}\,\theta(z+z^{\prime})(z+z^{\prime})[\tilde{\sigma}_{s}(z)\tilde{\sigma}_{s}(z^{\prime})+\tilde{\sigma}_{a}(z)\tilde{\sigma}_{a}(z^{\prime})].$
(S29)

If one now identifies

$\displaystyle\tilde{\sigma}$
$\displaystyle=\int\\!dz\,\tilde{\sigma}(z)-ik\\!\int\\!dz\,dz^{\prime}\,\theta(z+z^{\prime})(z+z^{\prime})[\tilde{\sigma}_{s}(z)\tilde{\sigma}_{s}(z^{\prime})+\tilde{\sigma}_{a}(z)\tilde{\sigma}_{a}(z^{\prime})],$
(S30) $\displaystyle\tilde{\sigma}_{1}$
$\displaystyle=i\\!\int\\!dz\,z\tilde{\sigma}_{a}(z),$ (S31)

then equations (S28) and (S29) can be written as

$\displaystyle\frac{1}{T}$
$\displaystyle=\left(1+\frac{\tilde{\sigma}}{2}\right)(1-k\tilde{\sigma}_{1}),$
(S32) $\displaystyle\frac{R}{T}$
$\displaystyle=-\frac{\tilde{\sigma}}{2}(1-k\tilde{\sigma}_{1}),$ (S33)

which coincide, to order $O(k)$, with Eq. (17). For negative helicity $-$, it
is easy to find from Eq. (S6) that $\tilde{\sigma}$ remains the same, but
$\tilde{\sigma}_{1}$ flips sign.

##  II Twisted bilayer graphene

We compute the orbital magnetic moment of the quasiparticles twisted bilayer
graphene. We use the continuum theory Bistritzer and MacDonald (2011);
Tarnopolsky  _et al._ (2019). We follow the notation of Ref. Tarnopolsky  _et
al._ (2019). The system consists of two layers; the upper lay is rotated
clockwise by an angle $\theta>0$ with respect to the lower layer. The single-
particle Hamiltonian is

$H=\begin{pmatrix}-iv_{0}\bm{\sigma}_{\theta/2}\cdot\bm{\nabla}&T({\mathbf{x}})\\\
T^{\dagger}({\mathbf{x}})&-iv_{0}\bm{\sigma}_{-\theta/2}\cdot\bm{\nabla}\end{pmatrix}.$
(S34)

Here we denoted

$\bm{\sigma}_{\theta/2}=e^{-i\theta\sigma_{z}/4}\begin{pmatrix}\sigma_{x}\\\
\sigma_{y}\end{pmatrix}e^{i\theta\sigma_{z}/4}.$ (S35)

We also used

$T({\mathbf{x}})=\sum_{a=1}^{3}T_{a}e^{-i{\mathbf{q}}_{a}\cdot{\mathbf{x}}},$
(S36)

where we defined the momentum $\mathbf{q}_{a}$ as

${\mathbf{q}}_{1}=\begin{pmatrix}0\\\
-1\end{pmatrix}k_{\theta},\qquad{\mathbf{q}}_{2}=\begin{pmatrix}{\sqrt{3}}/2\\\
1/2\end{pmatrix}k_{\theta},\qquad{\mathbf{q}}_{3}=\begin{pmatrix}-{\sqrt{3}}/2\\\
1/2\end{pmatrix}k_{\theta},\qquad$ (S37)

and

$T_{1}=\begin{pmatrix}w_{\rm AA}&w_{\rm AB}\\\ w_{\rm AB}&w_{\rm
AA}\end{pmatrix},\quad T_{2}=\begin{pmatrix}w_{\rm AA}&w_{\rm AB}e^{-2\pi
i/3}\\\ w_{\rm AB}e^{2\pi i/3}&w_{\rm AA}\end{pmatrix},\quad
T_{3}=T_{2}^{\dagger}.$ (S38)

The matrix $T_{a}$ can be written as

$T_{a}=w_{\rm AA}+w_{\rm
AB}(\mathbf{\hat{z}}\times\hat{\mathbf{q}}_{a})\cdot\bm{\sigma},$ (S39)

where $\hat{\mathbf{q}}_{a}={\mathbf{q}}_{a}/k_{\theta}$ is the unit vector
directed along the direction of ${\mathbf{q}}_{a}$. It is convenient to rotate
the basis to transform the Hamiltonian into

$H\to\begin{pmatrix}e^{i\frac{\theta}{4}\sigma_{z}}&0\\\
0&e^{-i\frac{\theta}{4}\sigma_{z}}\end{pmatrix}H\begin{pmatrix}e^{-i\frac{\theta}{4}\sigma_{z}}&0\\\
0&e^{i\frac{\theta}{4}\sigma_{z}}\end{pmatrix}=\begin{pmatrix}-iv_{0}\bm{\sigma}\cdot\bm{\nabla}&T_{\theta}({\mathbf{x}})\\\
T_{\theta}^{\dagger}({\mathbf{x}})&-iv_{0}\bm{\sigma}\cdot\bm{\nabla}\end{pmatrix},$
(S40)

where

$T_{\theta}=\sum_{a}T_{\theta}^{a}e^{-i{\mathbf{q}}_{a}\cdot{\mathbf{x}}},\quad
T_{\theta}^{a}=w_{\rm AA}e^{i\theta\sigma_{z}/2}+w_{\rm
AB}(\mathbf{\hat{z}}\times\hat{\mathbf{q}}_{a})\cdot\bm{\sigma}.$ (S41)

To find the orbital magnetic moment of a state $|\psi\rangle$, one can turn on
a small in-plane magnetic field ${\mathbf{B}}$ by turning on opposite gauge
potentials on the two layers,

${\mathbf{A}}_{u}=-\frac{d}{2}(\mathbf{\hat{z}}\times{\mathbf{B}}),\quad{\mathbf{A}}_{d}=\frac{d}{2}(\mathbf{\hat{z}}\times{\mathbf{B}}),$
(S42)

where $u$ and $d$ indicates upper and lower layers, and $d$ is the distance
between the two layers. The interaction of the system with the gauge field
${\mathbf{A}}_{u}$ and ${\mathbf{A}}_{d}$ is given by

$\delta
H=-\frac{1}{c}({\mathbf{J}}_{u}\cdot{\mathbf{A}}_{u}+{\mathbf{J}}_{d}\cdot{\mathbf{A}}_{d})=\frac{1}{2c}({\mathbf{J}}_{u}-{\mathbf{J}}_{d})\cdot(\mathbf{\hat{z}}\times{\mathbf{B}}).$
(S43)

By comparing this with $\delta H=-\bm{\mu}\cdot{\mathbf{B}}$, we find

$\bm{\mu}=\frac{d}{2c}\mathbf{\hat{z}}\times({\mathbf{J}}_{u}-{\mathbf{J}}_{d})=-\frac{ed}{2c}\mathbf{\hat{z}}\times(\bm{\sigma}_{u}-\bm{\sigma}_{d}).$
(S44)

In the regime $w_{\rm AA},w_{\rm AB}\ll v_{0}k_{\theta}$, we can use
perturbation theory over $T$. To first order in perturbation the eigenstate of
a Hamiltonian $H=H_{0}+H_{1}$ is

$|0\rangle^{\prime}=|0\rangle-\sum_{n\neq 0}\frac{|n\rangle\langle
n|H_{1}|0\rangle}{E_{n}-E_{0}}+\mathcal{O}(T^{2}).$ (S45)

If $|0\rangle$ is a state on the upper layer with $p\approx 0$, then
$|n\rangle$ are states on the lower layer with momentum
${\mathbf{p}}+{\mathbf{q}}_{a}\approx q_{a}$, $a=1,2,3$. Thus we need to
consider only 8 bands $|0,\alpha\rangle$ and $|a,\alpha\rangle$ with
$\alpha=1,2$ being the spinor index. The unperturbed Hamiltonian is

$\langle
0,\alpha|H_{0}|0,\beta\rangle=v_{0}{\mathbf{p}}\cdot\bm{\sigma}_{\alpha\beta},\quad\langle
a,\alpha|H_{0}|b,\beta\rangle=v_{0}({\mathbf{q}}_{a}+{\mathbf{p}})\cdot\bm{\sigma}_{\alpha_{\beta}}\delta_{ab},$
(S46)

while the perturbation part of the Hamiltonian is

$\langle 0,\alpha|H_{1}|a,\beta\rangle=(T^{\theta}_{a})_{\alpha\beta}.$ (S47)

The expectation value of the magnetic moment is then

$\bm{\mu}=-\frac{d}{2c}\mathbf{\hat{z}}\times\left(\langle
0|\bm{\sigma}_{u}|0\rangle-\sum_{mn}\frac{\langle 0|H_{1}|m\rangle\langle
m|\bm{\sigma}_{d}|n\rangle\langle
n|H_{1}|0\rangle}{(E_{m}-E_{0})(E_{n}-E_{0})}\right).$ (S48)

Note that $\langle 0|\bm{\sigma}|0\rangle=\hat{\mathbf{p}}$, and therefore
does not contribute to the magnetic helicity. Thus, we only need to compute
the second term inside the bracket In each sector with fix $a$, $H_{0}$ has
two energy eigenvectors

$H_{0}|a,\pm\rangle=\pm v_{0}p_{a}|a,\pm\rangle,\qquad
p_{a}\equiv|{\mathbf{q}}_{a}+{\mathbf{p}}|.$ (S49)

The sum over intermediate states in $a$-sector gives

$S_{a}\equiv\sum_{\alpha=\pm}\frac{|a,\alpha\rangle\langle
a,\alpha|}{E_{\alpha}-E_{0}}=\frac{1}{2v_{0}(p_{a}-p)}(1+\bm{\sigma}\cdot\hat{\mathbf{p}}_{a})-\frac{1}{2v_{0}(p_{a}+p)}(1-\bm{\sigma}\cdot\hat{\mathbf{p}}_{a}).$
(S50)

Putting everything into the expression for the orbital magnetic moment, we
then find

${\bm{v}}_{\mathbf{p}}\cdot\bm{\mu}_{\mathbf{p}}=\frac{d}{2c}{\bm{v}}_{\mathbf{p}}\cdot\sum_{a=1}^{3}\langle
0|T_{a}^{\theta}S_{a}(\mathbf{\hat{z}}\times\bm{\sigma})S_{a}(T_{a}^{\theta})^{\dagger}|0\rangle.$
(S51)

For $p\ll k_{\theta}$, this can be evaluated to give

$\langle{\bm{v}}_{\mathbf{p}}\cdot\bm{\mu}_{\mathbf{p}}\rangle=\frac{6w_{\rm
AA}w_{\rm
AB}}{(v_{0}k_{\theta})^{2}}\frac{p}{k_{\theta}}\frac{edv_{0}^{2}}{2c}+\mathcal{O}(p^{2}).$
(S52)

The presence of a slip (vector $\mathbf{d}$ in Ref. Bistritzer and MacDonald
(2011)) changes the overall phase of the matrix $T({\mathbf{x}})$ and hence
according to Eq. (S51) does not change the result for the magnetic helicity.

##  III Faraday rotation of the superconducting phase with a finite incident
angle $\alpha$

###  III.1 Setup and equation of motions

We consider a thin chiral superconductor with thickness $2\varepsilon$. The
scattering light lives on the $(y,z)$ plane and makes an angle $\alpha$ with
the $z$ axis as shown in Fig S1. We recall the action of a thin chiral
superconductor interacting with the electromagnetic field

$\displaystyle\mathcal{S}=\int\\!d^{4}x\left\\{-\frac{1}{16\pi}F_{\mu\nu}F^{\mu\nu}+\delta(z)\left[\frac{f^{2}}{2}(D_{t}\varphi)^{2}-\frac{n_{s}}{2m}(D_{a}\varphi)^{2}-n_{s}\kappa
B_{a}D_{a}\varphi\right]\right\\},$ (S53)

Figure S1: Light scattering diagram

where we used $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ and
$\partial_{0}=\frac{1}{c}\partial_{t}$. We chose subscripts $a,b$ to denote
inplane spacial directions. We chose the Coulomb gauge $A_{0}=0$, instead of
$\partial_{a}\varphi=0$ for $\alpha=0$ as in the main text, the equation of
motion of $\varphi$ is non-trivial. The components of incident, reflection and
transmission lights has the following ansatz

$\displaystyle A_{x}^{I}e^{-i(\omega t-\mathbf{k}\cdot\mathbf{r})},\quad
A_{y}^{I}e^{-i(\omega t-\mathbf{k}\cdot\mathbf{r})},\quad
A_{z}^{I}e^{-i(\omega t-\mathbf{k}\cdot\mathbf{r})},$ (S54) $\displaystyle
A_{x}^{R}e^{-i(\omega t-\mathbf{k}^{\prime}\cdot\mathbf{r})},\quad
A_{y}^{R}e^{-i(\omega t-\mathbf{k}^{\prime}\cdot\mathbf{r})},\quad
A_{z}^{R}e^{-i(\omega t-\mathbf{k}^{\prime}\cdot\mathbf{r})},$ (S55)
$\displaystyle A_{x}^{T}e^{-i(\omega t-\mathbf{k}\cdot\mathbf{r})},\quad
A_{y}^{R}e^{-i(\omega t-\mathbf{k}\cdot\mathbf{r})},\quad
A_{z}^{R}e^{-i(\omega t-\mathbf{k}\cdot\mathbf{r})}.$ (S56)

where we denote $\mathbf{k}=(0,k_{y},k_{z})$ and
$\mathbf{k}^{\prime}=(0,k_{y},-k_{z})$, with
$k_{y}=k\sin\alpha,k_{z}=k\cos\alpha$. Our task is finding $A^{R}_{i}$ and
$A^{T}_{i}$ in term of $A^{I}_{i}$. The equation of motion (e.o.m) of $A_{0}$
gives

$0=\frac{1}{c}\partial_{t}\partial_{i}A_{i}+2\delta(z)f^{2}\partial_{t}\varphi.$
(S57)

Equation (S57) is nothing but the Gauss law which implies $\nabla_{i}A_{i}=0$
outside the superconductor. The e.o.m of $\varphi$ is

$0=\delta(z)\left[-f^{2}\partial_{t}^{2}\varphi+\frac{n_{s}}{m}\partial_{a}D_{a}\varphi+n_{s}\kappa\partial_{a}B_{a}\right].$
(S58)

Finally the e.o.m of $A_{i}$ are

$0=-\frac{1}{4\pi
c^{2}}\partial_{t}^{2}A_{z}-\frac{\epsilon^{ab}}{4\pi}\partial_{a}B_{b}-n_{s}\kappa\delta(z)\epsilon^{ab}\partial_{a}D_{b}\varphi,$
(S59) $0=-\frac{1}{4\pi
c^{2}}\partial_{t}^{2}A_{x}-\frac{(\partial_{y}B_{z}-\partial_{z}B_{y})}{4\pi}-\delta(z)\left(\frac{2n_{s}e}{m\hbar
c}D_{x}\varphi+2\frac{n_{s}e}{\hbar c}\kappa
B_{x}\right)+n_{s}\kappa\partial_{z}\left(\delta(z)D_{y}\varphi\right),$ (S60)
$0=-\frac{1}{4\pi
c^{2}}\partial_{t}^{2}A_{y}-\frac{(\partial_{z}B_{x}-\partial_{x}B_{z})}{4\pi}-\delta(z)\left(\frac{2n_{s}e}{m\hbar
c}D_{y}\varphi+2\frac{n_{s}e}{\hbar c}\kappa
B_{y}\right)-n_{s}\kappa\partial_{z}\left(\delta(z)D_{x}\varphi\right).$ (S61)

The combination of Gauss law (S57) and the e.o.m for $A_{i}$ gives us

$\frac{1}{c^{2}}\partial_{t}^{2}A_{i}-\bigtriangleup A_{i}=0$ (S62)

outside the superconductor. It is satisfied automatically with photon’s
dispersion relation $\omega=c|\mathbf{k}|$. In current setup, one can replace
$\partial_{x}\rightarrow 0$. We also replace $\partial_{y}\rightarrow ik_{y}$
and $\partial_{t}\rightarrow-i\omega$ to convert above e.o.m to the momentum
space. We employ the long-wave length limit approximations

$A_{i}(0)=\frac{1}{2}(A_{i}(\varepsilon)+A_{i}(-\varepsilon)),\qquad
A^{\prime}_{i}(0)=\frac{1}{2}(A^{\prime}_{i}(\varepsilon)+A^{\prime}_{i}(-\varepsilon)),$
(S63)

with the notation $A^{\prime}_{i}=\partial_{z}A_{i}$. Integrating equations
(S57)–(S61) from $-\varepsilon$ to $\varepsilon$, we obtain

$-i\omega\left((A_{z}(\varepsilon)-A_{z}(-\varepsilon))+2c4\pi
f^{2}\varphi\right)=0,$ (S64)
$f^{2}\omega^{2}\varphi-\frac{n_{s}}{m}\left[k_{y}^{2}\varphi-\frac{i}{c}k_{y}\left(A_{y}(\varepsilon)+A_{y}(-\varepsilon)\right)\right]+ik_{y}n_{s}\kappa\frac{1}{2}\left(A^{\prime}_{x}(\varepsilon)+A^{\prime}_{x}(-\varepsilon)\right)=0,$
(S65) $-ik_{y}(A_{y}(\varepsilon)-A_{y}(-\varepsilon))+ik_{y}\frac{4\pi
n_{s}e}{\hbar c}\kappa(A_{x}(\varepsilon)+A_{x}(-\varepsilon))=0,$ (S66)
$A^{\prime}_{x}(\varepsilon)-A^{\prime}_{x}(-\varepsilon)-4\pi\left\\{\frac{2n_{s}e^{2}}{m\hbar^{2}c^{2}}(A_{x}(\varepsilon)+A_{x}(-\varepsilon))-\frac{n_{s}e}{\hbar
c}\kappa\left[-ik_{y}(A_{z}(\varepsilon)+A_{z}(-\varepsilon))+(A^{\prime}_{y}(\varepsilon)+A^{\prime}_{y}(-\varepsilon))\right]\right\\}=0,$
(S67)
$A^{\prime}_{y}(\varepsilon)-A^{\prime}_{y}(-\varepsilon)-ik_{y}\left(A_{z}(\varepsilon)-A_{z}(-\varepsilon)\right)-4\pi\left\\{\frac{2n_{s}e}{m\hbar
c}\left[ik_{y}\varphi+\frac{e}{\hbar
c}(A_{y}(\varepsilon)+A_{y}(-\varepsilon))\right]+\frac{n_{s}e}{\hbar
c}\kappa(A^{\prime}_{x}(\varepsilon)+A^{\prime}_{x}(-\varepsilon))\right\\}=0.$
(S68)

Multiplying equations (S60)-(S61) by $z$ then integrating from $-\varepsilon$
to $\varepsilon$, we obtain

$-(A_{x}(\varepsilon)-A_{x}(-\varepsilon))+4\pi
n_{s}\kappa\left[-ik_{y}\varphi-\frac{e}{\hbar
c}(A_{y}(\varepsilon)+A_{y}(-\varepsilon))\right]=0,$ (S69)

$-(A_{y}(\varepsilon)-A_{y}(-\varepsilon))+4\pi\frac{n_{s}e}{\hbar
c}\kappa(A_{x}(\varepsilon)+A_{x}(-\varepsilon))=0.$ (S70)

We now have 9 equations including (S64)-(S70),and the Gauss law for
$z<-\varepsilon$ and $z>\varepsilon$

$A_{y}^{R}k_{y}-A_{z}^{R}k_{z}=0,\quad A_{y}^{T}k_{y}+A_{z}^{T}k_{z}=0,$ (S71)

to find 7 unknowns $\varphi,A^{R}_{i},A^{T}_{i}$ from $A^{I}_{i}$ and
$k_{y},k_{z}$ with the condition $A_{y}^{I}k_{y}+A_{z}^{I}k_{z}=0$ implied by
the Gauss law. One can check that 2 of above 9 equations are redundant. We can
solve uniquely $A_{i}^{R},A_{i}^{T}$ from $A_{i}^{I}$. We will summarize the
results in the next subsection.

###  III.2 Transmission of linear polarized light

Due to the Gauss law, we define in-plane vector potential as

$A_{in}^{I}=\frac{A_{y}^{I}}{\cos\alpha}=-\frac{A_{z}^{I}}{\sin\alpha},\quad
A_{in}^{R}=\frac{A_{y}^{R}}{\cos\alpha}=\frac{A_{z}^{R}}{\sin\alpha},\quad
A_{in}^{T}=\frac{A_{y}^{T}}{\cos\alpha}=-\frac{A_{z}^{T}}{\sin\alpha}.$ (S72)

We also define the out of plane vector potential as

$A_{out}^{I}=-A_{x}^{I},\quad A_{out}^{R}=A_{x}^{R},\quad
A_{out}^{T}=-A_{x}^{T}.$ (S73)

The polarization function of the transmission light is

$\mathcal{P}^{T}=\frac{A_{out}^{T}}{A_{in}^{T}}.$ (S74)

From the definition (S74), we see that real (pure imaginary) $\mathcal{P}$
correspond to linear (elliptical) polarized light. Especially, $\mathcal{P}=i$
($\mathcal{P}=-i$) corresponds to right-hand circle (left-hand circle)
polarized light. We consider the linear polarized incident light with either
in-plane polarization or out-of-plane. We calculate the polarization function
using equations in the previous section and quote the result at the leading
orders in the momentum $k$:

In-plane incident light ($A_{in}^{I}\neq 0,A_{out}^{I}=0$)

We obtain the result of polarization functions as

$\mathcal{P}^{T}=\frac{\kappa
mc\left(-64\pi^{2}\frac{f^{2}}{c^{2}}\frac{n_{s}e^{2}}{m\hbar^{2}c^{2}}\cos(\alpha)+ik\sin^{2}(\alpha)\left(\frac{4\pi
f^{2}}{c^{2}}\left(16\pi^{2}\kappa^{2}\frac{n_{s}^{2}e^{2}}{\hbar^{2}c^{2}}+1\right)-\frac{4\pi
n_{s}e^{2}}{m\hbar^{2}c^{2}}\right)\right)}{2\frac{f^{2}}{c^{2}}\left(16\pi^{2}\kappa^{2}\frac{n_{s}^{2}e^{2}}{\hbar^{2}c^{2}}-1\right)+2\frac{n_{s}e^{2}}{m\hbar^{2}c^{2}}\sin^{2}(\alpha)}+\mathcal{O}(k^{2}).$
(S75)

Out-of-plane incident light ($A_{in}^{I}=0,A_{out}^{I}\neq 0$)

We obtain the result of polarization functions as

$1/\mathcal{P}^{T}=\frac{\frac{\kappa m\hbar
c}{e\cos^{2}(\alpha)}\left(64\pi^{2}\frac{f^{2}}{c^{2}}\frac{n_{s}e^{2}}{m\hbar^{2}c^{2}}\cos(\alpha)+ik\sin^{2}(\alpha)\left(4\pi\frac{f^{2}}{c^{2}}\left(16\pi^{2}\kappa^{2}\frac{n_{s}^{2}e^{2}}{\hbar^{2}c^{2}}+1\right)-\frac{4\pi
n_{s}e^{2}}{m\hbar^{2}c^{2}}\right)\right)}{2\frac{f^{2}}{c^{2}}\left(16\pi^{2}\kappa^{2}\frac{n_{s}^{2}e^{2}}{\hbar^{2}c^{2}}-1\right)}\,.$
(S76)

From equation (S75) and (S76), at the incident angle $\alpha=0$, we obtain the
Faraday rotation angle of the linear polarized light (both in-plane and out-
of-plane)

$\theta_{F}=\arctan\left(\frac{8\pi\kappa\frac{n_{s}e}{\hbar
c}}{1-16\pi^{2}\kappa^{2}\frac{n_{s}^{2}e^{2}}{\hbar^{2}c^{2}}}\right)\approx
8\pi\frac{\kappa n_{s}e}{\hbar c},$ (S77)

which is the same as the Faraday rotation angle in the main text for the
superconducting phase by combining Eqs. (32) and (18) . With a general value
of $\alpha$, the transmission light is neither linearly nor elliptically
polarized.