# Neutral magic-angle bilayer graphene: Condon instability and chiral
resonances

T. Stauber1,2, M. Wackerl2, P. Wenk2, D. Margetis3, J. González4, G. Gómez-
Santos5, and J. Schliemann2 1Instituto de Ciencia de Materiales de Madrid,
CSIC, E-28049 Madrid, Spain
2Institut für Theoretische Physik, Universität Regensburg, Germany
3Institute for Physical Science and Technology, and Department of Mathematics,
and Center for Scientific Computation and Mathematical Modeling, University of
Maryland, College Park, Maryland 20742, USA
4Instituto de Estructura de la Materia, CSIC, E-28006 Madrid, Spain
5Departamento de Física de la Materia Condensada, Instituto Nicolás Cabrera
and Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid,
E-28049 Madrid, Spain

###### Abstract

We discuss the full optical response of twisted bilayer graphene at the
neutrality point close to the magic angle within the continuum model (CM).
Firstly, we identify three different channels consistent with the underlying
$D_{3}$ symmetry, yielding the total, magnetic, and chiral response. Secondly,
we numerically calculate the full optical response in the immediate vicinity
of the magic angle $\theta_{m}$ which provides a direct mapping of the CM onto
an effective two-band model. We, further, show that the ground-state of the CM
in the immediate vicinity of $\theta_{m}$ is unstable towards transverse
current fluctuations, a so-called Condon instability. Thirdly, due to the
large counterflow, the acoustic plasmonic excitations with typical wave
numbers have larger energies than the optical ones and their energy density
may be largely enhanced at certain frequencies which we denominate as chiral
resonances. Finally, we discuss symmetry relations for the optical response
and their consequences for the chiral response.

##  I Introduction

Twisted bilayer grapheneLopes dos Santos  _et al._ (2007); Shallcross  _et
al._ (2008); Suárez Morell  _et al._ (2010); Schmidt  _et al._ (2010); Li  _et
al._ (2010); de Laissardière  _et al._ (2010); Bistritzer and MacDonald
(2011); Lopes dos Santos  _et al._ (2012); Dean  _et al._ (2013) has attracted
much attention due to a plethora of new correlated phases such as correlated
insulators,Cao  _et al._ (2018a) unconventional superconductivityCao  _et al._
(2018b); Yankowitz  _et al._ (2019) or anomalous quantum Hall
ferromagnetism.Sharpe  _et al._ (2019); Polshyn  _et al._ (2020) Most features
are related to the emergence of a flat band which is related to the vanishing
of the Fermi velocity at the so-called magic angle $\theta_{m}$. In addition,
a prominent counterflow can be found where the current has opposite direction
with respect to the two layers and which becomes balanced at
$\theta_{m}$.Bistritzer and MacDonald (2011) A third feature is the pronounced
circular dichroismKim  _et al._ (2016); Morell  _et al._ (2017) at frequencies
close to the van Hove singularities.

Fermi velocity, counterflow and circular dichroism are related to the total
(electric), magnetic and chiral response. These were first introduced in Refs.
Stauber  _et al._ , 2018a, b and discussed in detail for large twist angles.
Here, we shall calculate these quantities for twist angles around the magic
angle.

We will also address flat band plasmonics in twisted bilayer grapheneHu  _et
al._ (2017) that is related to localized collective modes.Stauber and Kohler
(2016) This topic is an area of active interest.Lewandowski and Levitov
(2019); Novelli  _et al._ (2020); Hesp  _et al._ (2021); Huang  _et al._
(2022); Kuang  _et al._ (2021) Finally, we give a general discussion and new
insights on the chiral optical response and how it is related to symmetries.

The flat bandsBistritzer and MacDonald (2011); Tarnopolsky  _et al._ (2019);
Song  _et al._ (2019); Hejazi  _et al._ (2019); Koshino  _et al._ (2018); Park
_et al._ (2020) and the optical responseStauber  _et al._ (2013); Vela  _et
al._ (2018); Novelli  _et al._ (2020); Dai  _et al._ (2021); Han and Lai
(2022) of twisted bilayer graphene have been investigated in numerous articles
so far. However, a detailed discussion on the scaling behavior of the response
function for small frequencies, as $\omega\to 0$, around the immediate
vicinity of the magic angle, is missing until now. Another topic is related to
the Condon instabilityCondon and Walstedt (1968) that has recently been
discussed in related systems. Andolina  _et al._ (2020); Nataf  _et al._
(2019); Guerci  _et al._ (2020); Sánchez  _et al._ (2021) We will argue that
the Condon instability arises in the continuum modelLopes dos Santos  _et al._
(2007); Bistritzer and MacDonald (2011) of twisted bilayer graphene at the
neutrality point in the immediate vicinity of the magic angle.

Apart from the above, the large counterflow or magnetic response has also been
discussed in several papers.Bistritzer and MacDonald (2011); Stauber  _et al._
(2018a); Guerci  _et al._ (2021) Nevertheless, the importance of acoustic
plasmonic modes has not received sufficient attention, so far. We believe that
our results will be relevant to flat-band plasmonics in twisted bilayer
graphene, especially at the chiral resonance where the energy density can be
largely enhanced.

Finally, the chirality in graphene might be used to design novel cavities that
lead to strong chiral light-matter interaction.Stauber  _et al._ (2020a) In
order to understand the underlying physics we point out some new aspects
related to particle-hole symmetry. This leads us to distinguish between
electron and hole transitions where the initial states are energetically
closer and further away from the neutrality point than the final states,
respectively.

The paper is organized as follows. In Sec. II, we define the continuum model
for twisted bilayer graphene. In Sec. III, we introduce the minimal model for
the full linear response that defines the total, magnetic, and chiral
responses. In Sec. IV, we present our numerical results for the dissipative
and reactive response. In Sec. V, we carry out this task, albeit in the
immediate vicinity of the magic angle. By this procedure, we obtain a scaling
relation that eventually leads to the prediction of a Condon instability. In
Sec. VI, we discuss flat-band plasmonics without excess charges for twist
angles near the magic angle, and highlight the fact that genuinely acoustic
plasmons might be dominating the plasmonic properties around the magic angle.
In Sec. VII, we outline the symmetry conditions for chiral response. Finally,
we conclude our paper with Sec. VIII, a summary and outlook. Supplementary
Information on the numerical recipe of how to calculate the dissipative
response in the clean limit as well as on an analytical calculation of the
optical conductivity in the immediate vicinity of the magic angle is also
provided.

##  II Hamiltonian

The local Hamiltonian of a twisted bilayer graphene can be approximated byMoon
and Koshino (2013)

$\displaystyle
H=\left(\begin{matrix}H_{0}^{-\theta/2}&V^{\dagger}({\bm{r}})\\\
V({\bm{r}})&H_{0}^{\theta/2}\end{matrix}\right)\;,$ (1)

where $H_{0}^{\gamma}=-i\hbar v_{F}{\bm{\tau}}^{\gamma}\cdot\partial_{\bm{r}}$
denotes the Hamiltonian of the separate layers with
$(\bm{\tau}^{\gamma}_{x},\bm{\tau}^{\gamma}_{y})=e^{\rm
i\gamma\bm{\tau}_{z}/2}(\bm{\tau}_{x},\bm{\tau}_{y})e^{-\rm
i\gamma\bm{\tau}_{z}/2}$, $\bm{\tau}_{x,y,z}$ being the Pauli matrices. The
interlayer coupling $V({\bm{r}})$ also denotes a $2\times 2$-matrix and
defines the coupling between the layers. It explicitly depends on the stacking
order, but a common approximation is that all components are defined by only
one common function $u({\bm{r}})$.Lopes dos Santos  _et al._ (2007); Moon and
Koshino (2013) Expanding $u({\bm{r}})$ into the first three Fourier-components
of the moiré lattice and representing the Hamilton operator by plane waves,
one arrives at the non-interacting Hamiltonian used for calculating the
response to total fieldsLopes dos Santos  _et al._ (2007); Bistritzer and
MacDonald (2011)

$\displaystyle\mathcal{H}=$ $\displaystyle\hbar
v_{F}\sum_{\bm{k};\alpha,\beta}c_{\bm{k},\alpha,1}^{\dagger}\;\bm{\tau}_{\alpha\beta}^{-\theta/2}\cdot\bm{k}\;c_{\bm{k},\beta,1}$
$\displaystyle+$ $\displaystyle\hbar
v_{F}\sum_{\bm{k};\alpha,\beta}c_{\bm{k},\alpha,2}^{\dagger}\;\bm{\tau}_{\alpha\beta}^{+\theta/2}\cdot\bm{k}\;c_{\bm{k},\beta,2}$
$\displaystyle+$
$\displaystyle\frac{t_{\perp}}{3}\sum_{\bm{k};\alpha,\beta;\bm{G}}(c_{\bm{k}+\bm{G},\alpha,1}^{\dagger}\;T_{\alpha\beta}(\bm{G})\;c_{\bm{k},\beta,2}+H.c.)\;,$
(2)

where the separation between twisted cones is
$\Delta\bm{K}=2|\bm{K}|\sin(\theta/2)\left(0,1\right)$ with
$\bm{K}=\tfrac{4\pi}{3a_{g}}\left(1,0\right)$. Interlayer hopping is
restricted to wavevectors $\bm{G}=\bm{0},\bm{G}_{1},\bm{G}_{2}$ with
$\bm{G}_{1}=|\Delta\bm{K}|\left(\tfrac{-\sqrt{3}}{2},-\tfrac{3}{2}\right)$,
$\bm{G}_{2}=|\Delta\bm{K}|\left(\tfrac{\sqrt{3}}{2},-\tfrac{3}{2}\right)$, and

$T(\bm{0})=\left(\begin{matrix}\kappa&1\\\
1&\kappa\end{matrix}\right);\;\;T(\bm{G}_{1})=T^{*}(\bm{G}_{2})=\left(\begin{matrix}\kappa
e^{\rm i2\pi/3}&1\\\ e^{-\rm i2\pi/3}&\kappa e^{\rm
i2\pi/3}\end{matrix}\right).$ (3)

Calculations are performed with $t=2.78\,\text{eV}$ and
$t_{\perp}=0.33\,\text{eV}$, being $\hbar v_{F}=\tfrac{\sqrt{3}}{2}ta_{g}$ the
Fermi velocity with graphene lattice constant
$a_{g}=2.46\,\mathring{\text{A}}$; the interlayer distance has been taken as
$a=3.5\,\mathring{\text{A}}$. In the first part of the work, we discuss the
symmetric model with $\kappa=1$ and in the second part of the work, the
asymmetric model introduced in Ref. Koshino  _et al._ , 2018 with $\kappa=0.8$
is used that accounts for out-of-plane relaxation, see also Ref. Guinea and
Walet, 2019.

Let us finally note that besides the parameter $\kappa$, the above model is
only characterized by one dimensionless parameter
$\alpha_{\theta_{i}}=\frac{\sqrt{A_{i}/3}}{2\pi}\frac{t_{\perp}}{t}$ which
combines $t_{\perp}$ and the twist angle $\theta_{i}$ parametrized by $i$ via
$A_{i}=3i^{2}+3i+1$ with $\cos(\theta_{i})=1-\tfrac{1}{2A_{i}}$.Bistritzer and
MacDonald (2011) This can readily be seen from the Hamiltonian of Eq. (1) by
introducing the dimensionless coordinates
$\tilde{\bm{r}}=|\Delta\bm{K}|{\bm{r}}$ such that the new interlayer coupling
between the layers is independent of the twist angle.Lopes dos Santos  _et
al._ (2007) In principle, $i\in\mathbb{N}$ denotes a commensurate twist angle,
but the expressions can be generalised to arbitrary real numbers
$i\in\mathbb{R}$.Bistritzer and MacDonald (2011)

##  III Linear response

To describe chiral effects without breaking time-reversal or rotational
($C_{3}$) symmetry, we have to treat an effectively three-dimensional system.
The minimal model thus consists of treating each layer of the twisted bilayer
separately. With the Kubo formula
$j_{\alpha}^{\ell}=-\chi_{j^{\ell}_{\alpha}\,j^{\ell^{\prime}}_{\beta}}A_{\beta}^{\ell^{\prime}}$,
where $A_{\alpha}^{\ell}$ denotes the gauge field and summation over repeated
indices is implied, the $4\times 4$ local (${\bm{q}}=0$) conductivity tensor
then is

$\sigma^{\ell,\ell^{\prime}}_{\alpha\beta}(\omega)=i\frac{\chi_{j^{\ell}_{\alpha}\,j^{\ell^{\prime}}_{\beta}}(\omega)}{\omega+i0^{+}},$
(4)

with axis indices $\alpha,\beta=x,y$ and plane indices
$\ell,\ell^{\prime}=1,2$.

The retarded current-current response is given by

$\chi_{j^{\ell}_{\alpha}\,j^{\ell^{\prime}}_{\beta}}(\omega)=g_{s}g_{v}\int_{1.BZ}\frac{d^{2}\bm{k}}{(2\pi)^{2}}\sum_{n,m}\frac{n_{F}(\epsilon_{m,\bm{k}})-n_{F}(\epsilon_{n,\bm{k}})}{\hbar\omega+i0^{+}-\epsilon_{n,\bm{k}}+\epsilon_{m,\bm{k}}}\langle
m,\bm{k}|j^{\ell}_{\alpha}|n,\bm{k}\rangle\langle
n,\bm{k}|j^{\ell^{\prime}}_{\beta}|m,\bm{k}\rangle\;.$ (5)

Here, $g_{s}=g_{v}=2$ are the spin and valley degeneracies. The states
$|m,\bm{k}\rangle$ are eigenstates of $\mathcal{H}$ in subband $m$ and of
momentum $\bm{k}$ in the first Brillouin zone of the superstructure. Their
eigenenergies are $\epsilon_{n,\bm{k}}$ and $n_{F}$ is the Fermi function. For
single layer graphene the current operator is $\bm{j}=-ev_{F}\bm{\tau}$ and
also for twisted bilayer graphene with the Hamiltonian of Eq. (II), the
general current operator is independent of ${\bm{k}}$.

The full current response due to an applied in-plane electric field that
satisfies rotational (or $C_{3}$) and time-reversal symmetry reads

$\left(\begin{matrix}j_{x}^{1}\\\ j_{y}^{1}\\\ j_{x}^{2}\\\
j_{y}^{2}\end{matrix}\right)=\left(\begin{matrix}\sigma_{0}&0&\sigma_{1}&\sigma_{xy}\\\
0&\sigma_{0}&-\sigma_{xy}&\sigma_{1}\\\
\sigma_{1}&-\sigma_{xy}&\sigma_{0}&0\\\
\sigma_{xy}&\sigma_{1}&0&\sigma_{0}\end{matrix}\right)\left(\begin{matrix}E_{x}^{1}\\\
E_{y}^{1}\\\ E_{x}^{2}\\\ E_{y}^{2}\end{matrix}\right)\;.$ (6)

The conductivities $\sigma_{\mu}=i\frac{\chi_{\mu}(\omega)}{\omega+i0^{+}}$
with $\mu=0,1,xy$ are defined via the following current-current response
functions:

$\displaystyle\chi_{0}$
$\displaystyle=\chi_{j^{1}_{x},j^{1}_{x}}=\chi_{j^{2}_{y},j^{2}_{y}}$ (7)
$\displaystyle\chi_{1}$
$\displaystyle=\chi_{j^{1}_{x},j^{2}_{x}}=\chi_{j^{1}_{y},j^{2}_{y}}=\chi_{j^{2}_{x},j^{1}_{x}}=\chi_{j^{2}_{y},j^{1}_{y}}$
(8) $\displaystyle\chi_{xy}$
$\displaystyle=\chi_{j^{1}_{x},j^{2}_{y}}=-\chi_{j^{1}_{y},j^{2}_{x}}=-\chi_{j^{2}_{x},j^{1}_{y}}=\chi_{j^{2}_{y},j^{1}_{x}}$
(9)

Note that the electric field may be different at the two layers and that the
above symmetries allow for different in-plane conductivities in layer 1 and 2
which may arise due to a perpendicular gate voltage.Stauber  _et al._ (2018a)
Also the influence of a perpendicular magnetic field can be included.Margetis
and Stauber (2021)

In the following we will discuss the response functions that transform with
respect to the irreducible representations of the underlying lattice symmetry
group $D_{3}$ which consists of two one-dimensional and one two-dimensional
representation. The response functions of the total (electronic) current
density ${\bm{j}}_{\rm{tot}}={\bm{j}}^{1}+{\bm{j}}^{2}$ and of the magnetic
current density ${\bm{j}}_{\rm{mag}}={\bm{j}}^{1}-{\bm{j}}^{2}$ transform as
the one-dimensional representations $A_{1}$ and $A_{2}$, respectively. These
current densities are induced by an in-plane electric and magnetic field,
respectively.Stauber  _et al._ (2018a); Sánchez  _et al._ (2021) The chiral
response involves the two current densities $j_{x}^{1}$ and $j_{y}^{2}$ which
transform as the two-dimensional representation $E$. This defines the total,
magnetic and chiral response, respectively, as

$\displaystyle\sigma_{\rm{tot}}$ $\displaystyle=2(\sigma_{0}+\sigma_{1})\;,$
(10) $\displaystyle\sigma_{\rm{mag}}$
$\displaystyle=2(\sigma_{0}-\sigma_{1})\;,$ (11)
$\displaystyle\sigma_{\rm{\rm{chi}}}$ $\displaystyle=\sigma_{xy}\;.$ (12)

This also defines the current-current response
$\chi_{\nu}=-i(\omega+i0^{+})\sigma_{\nu}$ and the Drude weight
$D_{\nu}=\lim_{\omega\to 0}\omega{\rm{Im}}\sigma_{\nu}(\omega)$ with
$\nu=\rm{tot},\rm{mag},\rm{\rm{chi}}$. The subindices $xy$ and $\rm{\rm{chi}}$
can be used interchangeably.

The above definitions also allow to deduce the exact symmetry relations when
the twist angle is reversed, i.e., when the opposite enantiomer is
considered:Stauber  _et al._ (2018a)

$\displaystyle\sigma_{0}(\theta)$ $\displaystyle=\sigma_{0}(-\theta)$ (13)
$\displaystyle\sigma_{1}(\theta)$ $\displaystyle=\sigma_{1}(-\theta)$ (14)
$\displaystyle\sigma_{xy}(\theta)$ $\displaystyle=-\sigma_{xy}(-\theta)$ (15)

It thus suffices to only consider the response for one twist-direction.

##  IV Optical response around the magic angle

The current response consists of a dissipative (imaginary) and reactive (real)
response. Numerically, the dissipative part is in principle equivalent to the
evaluation of a generalized density of states. Thus, we have

$\displaystyle{\rm{Im}}\chi^{\ell,\ell^{\prime}}_{\alpha\beta}(\omega)=\pi\frac{g_{s}g_{v}}{A}\sum_{\bm{k}}\sum_{n,m}\mathcal{O}_{m,n,{\bm{k}};\alpha\beta}^{\ell,\ell^{\prime}}\left[n_{F}(\epsilon_{m,\bm{k}})-n_{F}(\epsilon_{n,\bm{k}})\right]\delta(\hbar\omega-\epsilon_{n,\bm{k}}+\epsilon_{m,\bm{k}})\;,$
(16)

where we introduced the transition matrix element

$\mathcal{O}_{m,n,{\bm{k}};\alpha\beta}^{\ell,\ell^{\prime}}=\langle
m,\bm{k}|j^{\ell}_{\alpha}|n,\bm{k}\rangle\langle
n,\bm{k}|j^{\ell^{\prime}}_{\beta}|m,\bm{k}\rangle\;.$ (17)

Note that the matrix elements can always be considered as real since the final
spectral density must be real when time-reversal symmetry is not broken (which
is the case here).

Since the current response function is an analytic function in the upper
$\omega$-complex plane, see Eq. (5), the real part is obtained from the Cauchy
or Kramers-Kronig relation. By calculation of the principal value integral,
this procedure gives

$\displaystyle{\rm{Re}}\chi^{\ell,\ell^{\prime}}_{\alpha\beta}(\omega)=\frac{1}{\pi}\mathcal{P}\int_{-\infty}^{\infty}d\omega^{\prime}\frac{{\rm{Im}}\chi^{\ell,\ell^{\prime}}_{\alpha\beta}(\omega^{\prime})}{\omega^{\prime}-\omega}\;,$
(18)

which can be written as an integral over only positive frequencies using
${\rm{Im}}\chi^{\ell,\ell^{\prime}}_{\alpha\beta}(\omega)=-{\rm{Im}}\chi^{\ell,\ell^{\prime}}_{\alpha\beta}(-\omega)$.

Due to its linear dispersion, the continuum model does not possess a nominal
diamagnetic term. However, since the integral of Eq. (18) extends over all
frequencies, one needs to invoke a high-frequency cut-off and the resulting
term can be viewed as an effective diamagnetic term. For frequencies beyond
the cut-off, we assume that we have the current response of decoupled graphene
layers. Details of the regularization procedure can be found in Ref. Stauber
_et al._ , 2013.

Figure 1:  Dissipative response of the total (left), magnetic (center), and
chiral (right) current for the symmetric continuum model $\kappa=1$ at various
twist angles around the magic angle $\theta_{m}\simeq 1.03^{\circ}$. The real
conductivities Re$\sigma_{\rm{tot}}$, Re$\sigma_{\rm{mag}}$, and
Re$\sigma_{\rm{chi}}$ are given in units of the universal absorption of
single-layer graphene, $\sigma_{G}=\frac{e^{2}}{4\hbar}$. The inset of the
left panel shows the universal absorption of two graphene layers,
$2\sigma_{G}$, independent of the twist angle.

Furthermore, the complex conductivity can be obtained by first considering the
dissipative (real) part and then the reactive (imaginary) part. The only
difference is the Drude term that needs to be added to the dissipative part as
follows:

$\displaystyle{\rm{Re}}\sigma_{\alpha\beta}^{\ell,\ell^{\prime}}(\omega)=\pi
D_{\alpha\beta}^{\ell,\ell^{\prime}}\delta(\omega)+\frac{\pi}{\omega}\frac{g_{s}g_{v}}{A}\sum_{\bm{k}}\sum_{n,m}\mathcal{O}_{m,n,{\bm{k}};\alpha,\beta}^{\ell,\ell^{\prime}}\left[n_{F}(\epsilon_{m,\bm{k}})-n_{F}(\epsilon_{n,\bm{k}})\right]\delta(\omega-\epsilon_{n,\bm{k}}+\epsilon_{m,\bm{k}})\;,$
(19)

where

$\displaystyle D_{\alpha\beta}^{\ell,\ell^{\prime}}=\lim_{\omega\to
0}\chi_{j^{\ell}_{\alpha}\,j^{\ell^{\prime}}_{\beta}}(\omega)$ (20)

denotes the Drude weight matrix. According to the above equation, also the
total, magnetic, and chiral Drude weight can be defined according to Eqs.
(10-12).

For the dissipative response, we will only discuss the regular term of the
real conductivities which is characterized by plateaus in the low-frequency
limit $\omega\to 0^{+}$, i.e., we neglect the Drude weight which only
contributes for $\omega=0$ in the absence of intrinsic damping. These plateaus
are denoted as

$\displaystyle\sigma_{\rm{tot}}^{0}$ $\displaystyle=\lim_{\omega\to
0^{+}}{\rm{Re}}\sigma_{\rm{tot}}(\omega)\;,$ (21)
$\displaystyle\sigma_{\rm{mag}}^{0}$ $\displaystyle=\lim_{\omega\to
0^{+}}{\rm{Re}}\sigma_{\rm{mag}}(\omega)\;,$ (22)
$\displaystyle\sigma_{\rm{\rm{chi}}}^{0}$ $\displaystyle=\lim_{\omega\to
0^{+}}{\rm{Re}}\sigma_{xy}(\omega)\;.$ (23)

Let us finally note that due to the local response with ${\bm{q}}=0$, all
intraband contributions are contained in the Drude term, and the interband
term is often referred to as the regular contribution. In the Kramers-Kronig
relation, though, only the interband (regular) term enters due to
${\rm{Im}}\chi_{\alpha\beta}(\omega=0)=0$.

Figure 2:  Reactive response of the total (left), magnetic (center), and
chiral (right) current for the symmetric continuum model with $\kappa=1$ in
Eq. (3) at various twist angles around the magic angle $\theta_{m}\simeq
1.03^{\circ}$. The real current susceptibilities Re$\chi_{\rm{tot}}$,
Re$\chi_{\rm{mag}}$, and Re$\chi_{\rm{chi}}$ are given in units of
$t\frac{e^{2}}{\hbar^{2}}$ with $t=2.78$eV.

###  IV.1 Dissipative response

In Fig. 1, the dissipative response of the electric, magnetic and chiral
currents at the neutrality point is shown in terms of the real part of the
conductivity for angles around the magic angle $\theta\simeq 1.03^{\circ}$. It
is numerically obtained from Eq. (19) following the recipe outlined in the
Supplementary Information, see also Ref. Wiesenekker  _et al._ , 1988.

The total optical response for $\omega\to 0$ is characterized for all twist
angles by the universal conductivity of two uncoupled graphene layers,
$2\sigma_{G}$, with $\sigma_{G}=\frac{g_{s}g_{v}}{16}\frac{e^{2}}{\hbar}$; see
inset of the left panel of Fig. 1. This is an important consistency check as
in the low-frequency regime, the conductivity is entirely determined by the
Dirac cone.Nair  _et al._ (2008); Kuzmenko  _et al._ (2008); Falkovsky and
Pershoguba (2007); Stauber  _et al._ (2008) The $1/\omega$-prefactor of the
sums on the right-hand side of Eq. (19) is thus compensated by the weight of
the corresponding Fermi line whose circumference is also proportional to
$\omega$. In addition, two other universal plateaus emerge at larger
frequencies; see the left panel of Fig. 1. These features will be discussed in
Sec. V.

In the center panel of Fig. 1, the dissipative magnetic conductivity is shown.
As in the case of the total response, there are plateau for $\omega\to 0$ that
strongly increase around the magic angle, reaching values larger than
$10^{6}\sigma_{G}$. This might also be the origin of the large orbital
$g$-factor seen experimentally in twisted bilayer graphene.Li  _et al._
(2020); Sharpe  _et al._ (2021); Tschirhart  _et al._ (2021)

In the right panel of Fig. 1, the dissipative chiral conductivity is shown.
Again, there are plateaus marked by the Dirac regime whose values change sign
at $\theta\simeq 1.08^{\circ}$. Interestingly, this is the angle where the
spectrum displays an approximate $C_{6}$-symmetry at each valley which renders
the chiral Drude weight zero even for relatively large finite chemical
potential $|\mu|\lesssim 75$meV as discussed in Ref. Stauber  _et al._ ,
2020b.

###  IV.2 Reactive response

In Fig. 2, the reactive response of the total, magnetic and chiral current at
the neutrality point is shown for angles around the magic angle. It is
obtained from Eq. (16) via the Kramers-Kronig relation of Eq. (18). For this,
the dissipative part needs to be determined up to a frequency
$\omega_{\Lambda}$ for which
${\rm{Re}}\sigma_{\rm{tot}}(\omega\gtrsim\omega_{\Lambda})\approx
2\sigma_{G}$,
${\rm{Re}}\sigma_{\rm{mag}}(\omega\gtrsim\omega_{\Lambda})\approx
2\sigma_{G}$,
${\rm{Re}}\sigma_{\rm{chi}}(\omega\gtrsim\omega_{\Lambda})\approx 0$.Stauber
_et al._ (2013) These high-frequency limits represent the response of two
uncoupled layers and are also a consequence of the optical sum-rule.Sabio  _et
al._ (2008)

In the left panel of Fig. 2, the real part of the total current response is
shown. It must be zero for $\omega\to 0$ as there is no excess charge in the
system,Stauber  _et al._ (2020b) and we can adjust small numerical
errors.111Due to our numerical procedure, there is some uncertainty in
defining the cutoff-frequency and values between
$D_{\rm{tot}}=-0.001t\frac{e^{2}}{\hbar^{2}}$ ($\theta=1.3^{\circ}$) and
$D_{\rm{tot}}=0.004t\frac{e^{2}}{\hbar^{2}}$ ($\theta=1.0^{\circ}$) are
obtained. These shifts are also introduced to $\chi_{\rm{mag}}$ even though
this does hardly have an effect as the absolute values are much higher.

In the center panel of Fig. 2, the real part of the magnetic current response
is shown. We note that there is a non-monotonic behavior with respect to the
twist angle, i.e., even though the dissipative magnetic response is peaked
around the magic angle $\theta_{m}\simeq 1.03^{\circ}$, see center panel of
Fig. 1, the reactive response is not peaked at magic angle, but reaches a
maximum around $\theta\simeq 1.1^{\circ}$. This is due to the fact that for
these angles, the magnetic response reaches very high values at finite
frequencies with $0.1\leq\omega\leq 1$meV that yield the large response due to
the integration of Eq. (18). In Sec. V, though, we will argue that there is a
finite domain of twist angles in the immediate vicinity of the magic angle for
which the magnetic current response becomes maximal and even diverges.

In the right panel of Fig. 2, the real part of the chiral current response is
shown. It must be zero for $\omega\to 0$ as there is no excess charge in the
system,Stauber  _et al._ (2020b) and we can adjust small numerical
errors.222Due to our numerical procedure, there is some uncertainty in
defining the cutoff-frequency and values between
$D_{\rm{chi}}=-0.0025t\frac{e^{2}}{\hbar^{2}}$ ($\theta=1.1^{\circ}$) and
$D_{\rm{chi}}=0.0025t\frac{e^{2}}{\hbar^{2}}$ ($\theta=1.0^{\circ}$) are
obtained. For $\theta=1.0^{\circ}$, the maximal values can be as large as
$\chi_{\rm{chi}}=0.15t\frac{e^{2}}{\hbar^{2}}$ at $\hbar\omega\cong 1.9\mu$eV.

###  IV.3 Discussion on the Condon instability

There has been considerable interest in finding systems with a symmetry-broken
ground-state due to photon-condensation, the so-called Condon
instability.Andolina  _et al._ (2020); Nataf  _et al._ (2019); Guerci  _et
al._ (2020) In bilayer systems, this instability can also be discussed by
calculating the magnetic response $D_{\rm{mag}}$. Within the random-phase
approximation, the response must reach a critical value $D_{\rm{mag}}^{C}$
with

$\displaystyle\frac{\mu_{0}a}{4}D_{\rm{mag}}^{C}=-1\;,$ (24)

where $\mu_{0}$ denotes the magnetic permeability.Sánchez  _et al._ (2021)

For AA-stacked graphene, this limit is reached due to the logarithmic
divergence of the magnetic susceptibility.Sánchez  _et al._ (2021) However,
the response of twisted bilayer graphene is generally too weak to reach the
instability, i.e., including damping, one obtains
$D_{\rm{mag}}=-6.6t\frac{e^{2}}{\hbar^{2}}$.Stauber  _et al._ (2018a, b) Our
refined calculations without damping now yield a significantly lower bound for
$\theta=1.1^{\circ}$ with $D_{\rm{mag}}=-36t\frac{e^{2}}{\hbar^{2}}$. In the
above units, this translates to $D_{\rm{mag}}\approx 0.008(\mu_{0}a)^{-1}$ and
we have $D_{\rm{mag}}/D_{\rm{mag}}^{C}\approx 0.002$. This is still far away
from a possible Condon transition. However, in Sec. V, we will find a Condon
instability in the immediate vicinity of the magic angle by employing a
scaling approach.

Figure 3:  Band-structure of the two flat bands around charge neutrality of
the continuum model with $\kappa=1$ in Eq. (3) for various twist angles around
the magic angle $\theta_{m}\sim 1.032^{\circ}$. In the left panel, the band-
structure with smallest band-width is shown. In the center panel, one can
observe the avoided crossings along the $KM$-direction, whereas in the right
panel the avoided crossings are along the $\Gamma K$-direction; see insets.
The arrows indicate the avoided crossings for $\theta=1.05^{\circ}$ and
$\theta=1.02^{\circ}$.

We can compare our results also with previously reported values for the
magnetic susceptibility.Guerci  _et al._ (2020) The static magnetic
susceptibility $\chi_{\rm{mag}}^{0}$ is directly related to the magnetic Drude
weight at the neutrality point and given by
$\chi_{\rm{mag}}^{0}=\frac{a^{2}}{4}D_{\rm{mag}}$.Stauber  _et al._ (2018b)
With $D_{\rm{mag}}=-6.6t\frac{e^{2}}{\hbar^{2}}$, this yields
$\chi_{\rm{mag}}^{0}=0.02\frac{\mu_{B}}{\rm{nm}^{2}T}$ with $\mu_{B}$ the Bohr
magneton. This value, obtained for $\kappa=1$, is slightly larger than the one
reported in Ref. Guerci  _et al._ , 2021 for the continuum model with
$\kappa=0.2$.

With $D_{\rm{mag}}\approx-36t\frac{e^{2}}{\hbar^{2}}$, see central panel of
Fig. 2, we obtain for the static magnetic susceptibility an even larger value
of $\chi_{\rm{mag}}^{0}=0.12\frac{\mu_{B}}{\rm{nm}^{2}T}$. This amounts to
$18\mu_{B}$ per moiré cell only due to the orbital motion of counter-
propagating electrons. This purely quantum mechanical effect is remarkable as
no charge excitations are involved.

##  V Optical response at the magic angle

In this section, we discuss the optical response in the immediate vicinity of
the magic angle, i.e., we will highly zoom into this region of possible twist
angles. As we will see, for any angle one can find an energy regime which is
still characterized by the Dirac cone, i.e., one will never be exactly at the
magic angle just as one can never approach an irrational number. Furthermore,
other plateaus develop which can be anticipated from the band structure which
shall be discussed before we describe the scaling relations.

###  V.1 Bands around the magic angle

In Fig. 3, the band structure for the symmetric model is shown for different
twist angles around the magic angle $\theta_{m}\simeq 1.032^{\circ}$ where the
linear dispersion (Fermi velocity) at the $K$-point vanishes.Song  _et al._
(2019); Hejazi  _et al._ (2019); Koshino (2019) Notice that this does not
coincide with the smallest band-width condition which would yield a magic
twist-angle of $\theta_{m}^{*}\simeq 1.11^{\circ}$. This can be appreciated on
the left panel of Fig. 3 where a new regime starts with an accidental crossing
on the $\Gamma M$-direction.

In the center and right panel of Fig. 3, we see the evolution towards the
magic angle from above and below, respectively. Most notably, there is an
avoided crossing that is moving closer to the $K$-point when approaching the
magic angle which is highlighted in the insets. For even smaller angles, a
stable band-inversion emerges with the avoided crossing moving outward and
eventually inward again to form the second magic angle. The evolution in
$\theta$ around the second magic angle at $\theta\sim 0.49^{\circ}$, however,
is qualitatively different.

###  V.2 Scaling in the immediate vicinity of the magic angle

The universal conductivity of graphene for small frequencies,
$\sigma_{G}=\frac{e^{2}}{4\hbar}$, is due to the perfect cancellation between
the transition-matrix element and the Fermi velocity.Nair  _et al._ (2008);
Kuzmenko  _et al._ (2008); Peres and Stauber (2008) This is also the case for
the total conductivity of twisted bilayer graphene for transitions around the
Dirac cones. Considering different quantities such as the magnetic absorption
related to ${\rm{Re}}\sigma_{\rm{mag}}$ or the chiral absorption related to
${\rm{Re}}\sigma_{\rm{chi}}$ will not show this cancellation and we expect the
following relations for $\omega\to 0$:Bistritzer and MacDonald (2011)

$\displaystyle\sigma_{\rm{mag}}^{0}=\sigma_{G}\left(\frac{v_{\rm{mag}}}{v_{F}}\right)^{2}\;,\;\sigma_{\rm{chi}}^{0}=\sigma_{G}\left(\frac{v_{\rm{chi}}}{v_{F}}\right)^{2}$
(25)

Above, we defined suitable velocities that characterize the magnetic and
chiral excitations.

As the Fermi velocity vanishes at the magic angle, Eq. (25) suggest that the
magnetic and chiral absorption diverge. Our numerical calculations confirm
precisely this, as can be seen in Fig. 4, where we show the three response
functions for twist angles below the magic angle $\theta_{m}\sim
1.032^{\circ}$. Whereas the absorption shows universal behavior, the magnetic
as well as the chiral absorption diverge.

In order to discuss the scaling behavior of these quantities, we introduce the
effective parameter

$\displaystyle\alpha=\frac{\alpha_{\theta_{i}}-\alpha_{m}}{\alpha_{m}}\approx\frac{\theta_{m}-\theta_{i}}{\theta_{i}}\;,$
(26)

with $\alpha_{\theta_{i}}=\frac{\sqrt{A_{i}/3}}{2\pi}\frac{t_{\perp}}{t}$ and
$\alpha_{m}\sim 0.605$ (for $i=31.54$).

First, we investigate the scaling of the Dirac regime $\epsilon_{D}$ that is
defined by the abrupt increase of the absorption from $2\sigma_{G}$ to
$12\sigma_{G}$. As is shown in the inset of the left panel of Fig. 4, there is
a linear behavior of $\ln\epsilon_{D}$ as function of $\ln\alpha$, leading to
$\epsilon_{D}\simeq 0.88\alpha^{\gamma_{\epsilon}}\mu$eV with
$\gamma_{\epsilon}\simeq 1.35\pm 0.04$.

Figure 4:  The optical response for twist angles below the magic angle
$\theta_{m}=1.032^{\circ}$. The insets shows the logarithm of the Dirac regime
(left panel) and the optical response functions $\sigma_{\rm{mag}}$ (center
panel) and $\sigma_{\rm{chi}}$ (right panel) at $\omega=0$ as function of the
effective parameter ln$\alpha$ defined in Eq. (26).

Along the same lines, we obtain the scaling behavior of the magnetic and
chiral absorption plateau for $\hbar\omega<\epsilon_{D}$ as

$\displaystyle\sigma_{\rm{mag}}^{0}=570\sigma_{G}\alpha^{-\gamma_{\rm{mag}}}\;,\;\sigma_{\rm{chi}}^{0}=0.48\sigma_{G}\alpha^{-\gamma_{\rm{chi}}}\;.$
(27)

with $\gamma_{\rm{mag}}\simeq 1.62\pm 0.03$ and $\gamma_{\rm{chi}}\simeq
1.41\pm 0.02$.

It is generally argued that the Fermi velocity scales linearly in
$\alpha$.Lopes dos Santos  _et al._ (2007); Bistritzer and MacDonald (2011);
Lopes dos Santos  _et al._ (2012); G. Catalina, B. Amorim, E. V. Castro, J. M.
V. P. Lopes and Peres (2019); Watson and Luskin (2021); Becker  _et al._
(2021) This implies that $v_{\rm{mag}}=24v_{F}\alpha^{1-\gamma_{\rm{mag}}/2}$.
In addition, the chiral velocity must also tend to zero at the magic angle as
$v_{\rm{chi}}=0.7v_{F}\alpha^{1-\gamma_{\rm{chi}}/2}$.

###  V.3 Condon instability at the magic angle

As mentioned in Sec. IV.3, in AA-stacked bilayer graphene there is a Condon
instability at $T\sim 0$. Since in twisted bilayer graphene, the electronic
wave functions at the magic angle are highly localized around the AA-stacked
islands,de Laissardière  _et al._ (2010) there might be the possibility of a
Condon instability in twisted bilayer graphene around $\theta_{m}$.Guerci  _et
al._ (2021)

The imaginary part of the conductivity or magnetic Drude weight is obtained
from the Kramers-Kronig relation which can be split into the following two
contributions:

$\displaystyle D_{\rm{mag}}$
$\displaystyle=\frac{2}{\pi}\left[2\omega_{\Lambda}\sigma_{G}-\int_{0}^{\omega_{\Lambda}}d\omega{\rm{Re}}\sigma_{\rm{mag}}(\omega)\right]$
(28)
$\displaystyle=\frac{2}{\pi}\left[2\omega_{\Lambda}\sigma_{G}-\left(\int_{0}^{\omega_{D}^{\rm{mag}}}+\int_{\omega_{D}^{\rm{mag}}}^{\omega_{\Lambda}}\right)d\omega{\rm{Re}}\sigma_{\rm{mag}}(\omega)\right]$
$\displaystyle=D_{\rm{mag}}^{*}+D_{\rm{mag}}^{reg}\;,$

where $\omega_{\Lambda}$ denotes the high-frequency cutoff and
$\omega_{D}^{\rm{mag}}>0$ is the smallest frequency after the van Hove
singularity for which
$\sigma_{\rm{mag}}^{0}=\sigma_{\rm{mag}}(\omega_{D}^{\rm{mag}})$, i.e., for
$\theta=1.03^{\circ}$, this gives $\hbar\omega_{D}^{\rm{mag}}\approx
0.01\mu$eV.

The second term $D_{\rm{mag}}^{reg}$ is assumed to be regular. The possible
divergent contribution at the magic angle, $D_{\rm{mag}}^{*}$, can be
estimated as follows:

$\displaystyle
D_{\rm{mag}}^{*}=-\frac{2}{\pi}\int_{0}^{\omega_{D}^{\rm{mag}}}d\omega\sigma_{\rm{mag}}(\omega)\sim-\alpha^{-\gamma_{\rm{mag}}+\gamma_{\epsilon}^{\rm{mag}}}$
(29)

The exponent $\gamma_{\epsilon}^{\rm{mag}}$ is again obtained from a linear
fit of a log-log plot and we obtain $\gamma_{\epsilon}^{\rm{mag}}\simeq
1.41\pm 0.04$. We thus find a divergence at the magic angle that scales like
$D_{\rm{mag}}^{*}\sim-\alpha^{-\delta_{\rm{mag}}}$ with
$\delta_{\rm{mag}}=0.21\pm 0.05$. Since the Condon instability is marked by
$D_{\rm{mag}}\sim D_{\rm{mag}}^{C}=\frac{4}{\mu_{0}a}$, see Eq. (24), there
will be a symmetry-broken ground-state with orbital magnetic domains at the
magic angle.

The presence of an instability due to transverse current fluctuations in a
non-interacting model of Eq. (II) is a remarkable result and we are not aware
of any other non-interacting model that exhibits a symmetry-broken ground-
state other than $AA$-stacked bilayer graphene.Sánchez  _et al._ (2021) Let us
finally note that the total chiral Drude weight $D_{\rm{chi}}$ has to vanish
at the neutrality point due to gauge symmetry.Stauber  _et al._ (2018b, 2020a)

###  V.4 Mapping to effective model

The absorption spectrum in the immediate vicinity of the magic angle can
approximately be understood from the universal conductivity formulaStauber
_et al._ (2015) of a general dispersion
$\epsilon_{\bm{k}}\propto|{\bm{k}}|^{\nu}$

$\displaystyle\sigma(\omega)=\frac{g_{s}g_{v}g_{\ell}g_{C_{3}}}{16}\nu\frac{e^{2}}{\hbar}=g_{\ell}g_{C_{3}}\sigma_{G}\;.$
(30)

Above, we introduced the usual spin, valley, and layer-degree of freedom, but
also a possible $g_{C_{3}}$ degeneracy which takes the value 3 in case of an
explicit 3-fold degeneracy (otherwise it is 1). In the following, we will
discuss the results in units of the universal conductivity of graphene
$\sigma_{G}=\frac{g_{s}g_{v}}{16}\frac{e^{2}}{\hbar}$. Notice that we
introduce here explicitly the degeneracy factors which are usually set to
$g_{s}=g_{v}=2$.

At low frequencies, there will in principle always be a regime where the
absorption is governed by the universal absorption of Dirac Fermions with
$\nu=1$ and we have $\sigma(\omega)=2\sigma_{G}$. For twist angles in the
immediate vicinity of the magic angle, the plateau of a single quadratic
dispersion relation with $\nu=2$ is obtained with
$\sigma(\omega)=4\sigma_{G}$, seen in the left panel of Fig. 4 for
$\theta=1.03^{\circ}$ for $0.01\mu$eV$\lesssim\epsilon\lesssim 1\mu$eV.

Between these plateaus, a new plateau emerges with $\sigma_{0}=12\sigma_{G}$,
because a new absorption channel opens at the frequency of the avoided
crossing as seen in the inset of the center and right panels of Fig. 3 and
marked by arrows for $\theta=1.05^{\circ}$ and $\theta=1.02^{\circ}$. Even
though the band minima are elongated, as a first approximation they can be
assumed to be a quadratic dispersion and due to the $C_{3}$-symmetry, there
are three of them for each Dirac point. We thus numerically obtain
$\sigma(\omega)=12\sigma_{G}$.333This plateau $12\sigma_{G}$ is only obtained
for twist angles which are already sufficiently close to the magic angle,
i.e., the band structure for $\theta=1.05^{\circ}$ shows an avoided crossing,
but does not reach this plateau, yet.

However, we have been neglecting the contribution of the central Dirac cone
and the above qualitative discussion can be made quantitative by considering
the following two-band model which was first introduced in Refs. Bena and
Simon, 2011; Montambaux, 2012; Hejazi  _et al._ , 2019 for $M=0$:

$\displaystyle
H_{\bm{k}}=\frac{1}{2m}\begin{bmatrix}M&\varpi^{2}+\eta\varpi^{*}\\\
{\varpi^{*}}^{2}+\eta\varpi&-M\end{bmatrix}$ (31)

where $\varpi=\hbar(k_{x}-ik_{y})$. The model has eigenenergies
$2m\epsilon_{\bm{k}}=\pm\sqrt{M^{2}+k^{4}+2k^{3}\cos(3\theta)\eta+k^{2}\eta^{2}}$
displaying trigonal warping and zeros at $|{\bm{k}}|=\eta$. For $M=0$, there
are three nodal points which lie in the directions $\theta=\frac{2\pi n}{3}$
($\eta<0$) and $\theta=\pi-\frac{2\pi n}{3}$ ($\eta>0$) with $n\in\mathbb{N}$.
This transition can be also seen in the center and right panel of Fig. 3,
where the avoided crossing changes from the $KM$-direction (right from the
$K$-point) to the $\Gamma M$-direction (left from the $K$-point), related by a
$60^{\circ}$-rotation.

As shown in the Supplementary Information, the above model with $M=0$ yields
$\sigma=12\sigma_{G}$ for small frequencies and $\sigma=4\sigma_{G}$ for large
frequencies. The reason for not obtaining the Dirac regime
$\sigma=2\sigma_{G}$ is because the model of Eq. (31) with $M=0$ does not
exhibit a gap at the three nodal points with $|{\bm{k}}|=\eta$.

This can partially be remedied by introducing a ${\bm{k}}$-dependent mass term
with $M=|\hbar{\bm{k}}|^{2}$ such that the gap or Dirac-regime energy is given
by $\epsilon_{D}=\frac{\hbar^{2}}{m}\eta^{2}$. From the numerical approach we
obtain $\epsilon_{D}=0.88\alpha^{\gamma_{\epsilon}}\mu$eV. At the magic angle
$\theta_{m}\sim 1.03^{\circ}$, we can further extract the mass-term since
$\eta=0$. Remarkably, we get $m\approx m_{0}$ where $m_{0}$ is the mass of
free electrons. This allows us to connect $\eta$ to $\alpha$ of Eq. (26):

$\displaystyle\hbar\eta=\sqrt{0.88m\mu\rm{eV}}\alpha^{\gamma_{\epsilon}/2}$
(32)

Notice that with the discussion of the dimensionless energy scale
$\tilde{\omega}$ defined in the Supplementary Information, we would obtain the
same scaling relation. Eq. (32) together with $m\approx m_{0}$ provides a
direct mapping between the continuum model of twisted bilayer graphene and the
model of Eq. (31) in the immediate vicinity of the magic angle in the flat-
band regime.

Figure 5:  The real part of the current susceptibility Re$\chi_{\nu}(\omega)$
with $\nu=\rm{tot},\rm{mag},\rm{chi}$ of the asymmetric continuum model with
$\kappa=0.8$ in Eq. (II) at the neutrality point in units of
$t\frac{e^{2}}{\hbar^{2}}$ for temperatures $T=0,10$K. The optical gap is
indicated by the white area. Left panel: twist angle $\theta=1.3^{\circ}$.
Center panel: twist angle $\theta=1.2^{\circ}$. Right panel: twist angle
$\theta=1.1^{\circ}$.

##  VI Flat-band plasmonics

Since twisted bilayer graphene consists of two layers, there will be two
plasmonic modes. For layers far away, these modes hardly hybridize, but for an
interlayer distance $a=3.5\,\mathring{\text{A}}$, anti-bonding and bonding
modes emerge. Due to the long-ranged Coulomb interaction, the dispersions show
square-root and linear behavior in the momentum $q$ and define the so-called
optical (charge even) and acoustic (charge odd) branches, respectively. In the
local approximation, they are generally given by

$\displaystyle\omega_{+}^{2}$
$\displaystyle=\frac{\chi_{\rm{tot}}(\omega_{+})q}{2\epsilon_{0}\epsilon}\;,$
(33) $\displaystyle\omega_{-}^{2}$
$\displaystyle=\frac{\chi_{\rm{mag}}(\omega_{-})aq^{2}}{2\epsilon_{0}}\;,$
(34)

which define self-consistent equations for the plasmonic frequencies
$\omega_{+}$ and $\omega_{-}$ with momentum $q$, respectively. Note that the
optical mode depends on the dielectric environment through
$\epsilon=(\epsilon_{up}+\epsilon_{down})/2$, but the acoustic mode does
not.Stauber (2014)

The plasmon dispersion does not depend on the chiral Drude weight since the
non-retarded approximation does not allow for a coupling of longitudinal and
transverse modes.Stauber  _et al._ (2018a); Lin  _et al._ (2020) Nevertheless,
the optical (acoustic) mode, usually defined by electric (magnetic) dipole
oscillations, is now accompanied by parallel magnetic (electric) dipole
oscillations. With the magnetic dipole related to the magnetic current as
$2{\bm{m}}=a{\bm{j}}_{\rm{mag}}\times{\bm{e}}_{z}$, this is expressed by the
following relations:

$\displaystyle{\bm{e}}_{\bm{q}}\cdot{\bm{m}}$
$\displaystyle=-aX_{\rm{tot}}{\bm{e}}_{\bm{q}}\cdot{\bm{j}}_{\rm{tot}}\;,$
(35) $\displaystyle 4a{\bm{e}}_{{\bm{q}}_{\perp}}\cdot{\bm{j}}_{\rm{tot}}$
$\displaystyle=-X_{\rm{mag}}{\bm{e}}_{{\bm{q}}_{\perp}}\cdot{\bm{m}}\;,$ (36)

with

$\displaystyle
X_{\rm{tot}}=\frac{\chi_{\rm{chi}}}{\chi_{\rm{tot}}}\;,\;X_{\rm{mag}}=\frac{\chi_{\rm{chi}}}{\chi_{\rm{mag}}}\;.$
(37)

The total current is related to the electric dipole,
${\bm{j}}_{\rm{tot}}=-i\omega{\bm{p}}$, and we have
${\bm{p}}\parallel{\bm{m}}\parallel{\bm{q}}$ for the optical mode and
${\bm{p}}\parallel{\bm{m}}\perp{\bm{q}}$ for the acoustic mode.

The above relations are obtained from the transport equations of Eq. (6) and
hold also in the static limit, i.e., the total Drude weight $D_{\rm{tot}}$ and
chiral Drude weight $D_{\rm{chi}}$ are Fermi-line properties as discussed in
Ref. Stauber  _et al._ , 2020b. Similar conclusions have been drawn in Refs.
He  _et al._ , 2020; Antebi  _et al._ , 2022.

Let us finally note, that Eqs. (33) and (34) can be generalized to a non-local
approximation by the replacements
$\chi_{\rm{tot}}(\omega)\to\chi_{\rm{tot}}(\omega,{\bm{q}})$ and
$\chi_{\rm{mag}}(\omega)\to\chi_{\rm{mag}}(\omega,{\bm{q}})$ that leads to
flat plasmonic bands.Stauber and Kohler (2016)

###  VI.1 Poynting vector

Even though the optical and acoustic plasmon dispersions only depend on
$\chi_{\rm{tot}}$ and $\chi_{\rm{mag}}$, respectively, the Poynting vector
depends also on the chiral response $\chi_{\rm{chi}}$. To show this, let both
modes be induced by the sheet current $j_{\parallel}$ parallel to the plasmon
momentum ${\bm{q}}$, i.e., decomposing the Fourier components of the current
into longitudinal and transverse parts, we have
${\bm{j}}=j_{\parallel}{\bm{e}}_{\bm{q}}+j_{\perp}{\bm{e}}_{{\bm{q}}_{\perp}}$
for layer $\ell=1$.

For the optical mode, the sheet currents of the two layers are parallel and
for the acoustic mode, the sheet currents of the two layers are anti-parallel.
In the instantaneous approximation, the self-fields are purely longitudinal,
and we have ${\bm{q}}_{\perp}\cdot{\bm{E}}^{\ell}=0$ as well as
${\bm{q}}\cdot{\bm{E}}^{1}={\bm{q}}\cdot{\bm{E}}^{2}$ for the optical mode and
${\bm{q}}\cdot{\bm{E}}^{1}=-{\bm{q}}\cdot{\bm{E}}^{2}$ for the acoustic mode.
This yields the relation between the longitudinal and transverse current as
$j_{\perp}/j_{\parallel}=-2X_{\rm{tot}}$ and
$j_{\perp}/j_{\parallel}=2X_{\rm{mag}}$ for the two modes, respectively.

We then get in the non-retarded limit, close to the sheet and up to second
order in $qa$, the following expressions for the Poynting vectors of the
optical (tot) and acoustic (mag) mode (see also Ref. Stauber  _et al._ ,
2020a):

$\displaystyle\mathcal{P}_{\rm{tot}}$
$\displaystyle=\mathcal{P}_{0}\left(\begin{matrix}1+\gamma_{\rm{tot}}(qa)^{2}\\\
2sgn(z)X_{\rm{tot}}qa\\\ 0\end{matrix}\right)\;,$ (38)
$\displaystyle\mathcal{P}_{\rm{mag}}$
$\displaystyle=\mathcal{P}_{0}\left(\begin{matrix}4X_{\rm{mag}}^{2}\frac{k_{0}^{2}}{q^{2}}+\gamma_{\rm{mag}}(qa)^{2}\\\
-2sgn(z)X_{\rm{mag}}qa\\\ 0\end{matrix}\right)\;,$ (39)

with $\mathcal{P}_{0}=\frac{qj_{\parallel}^{2}}{2\epsilon_{0}\omega}$ and
$\gamma_{\nu}=[1+(4X_{\nu}^{2}-1)\frac{k_{0}^{2}}{q^{2}}]/4$, where
$k_{0}=\omega/c$ is the wavelength of light in free space and
$\nu=\rm{tot},\rm{mag}$. This shows that the chirality modifies the plasmonic
energy flux. Also note that the Poynting vector of the acoustic mode to lowest
order in $aq$ and in the non-chiral limit $X_{\rm{mag}}=0$ becomes zero since
this mode consists of perfectly cancelling counterpropagating current
densities.

Let us now discuss the limiting case $X_{\nu}\ll 1$ for $k_{0}/q\to 0$ and
$aq\to 0$. We then have for the Poynting vectors of the optical (tot) and
acoustic (mag) mode the following expressions:

$\displaystyle\mathcal{P}_{\rm{tot}}$
$\displaystyle=\mathcal{P}_{0}\left(\begin{matrix}1\\\
2sgn(z)X_{\rm{tot}}qa\\\ 0\end{matrix}\right)\;,$ (40)
$\displaystyle\mathcal{P}_{\rm{mag}}$
$\displaystyle=\mathcal{P}_{0}\left(\begin{matrix}0\\\
-2sgn(z)X_{\rm{mag}}qa\\\ 0\end{matrix}\right)\;.$ (41)

From the different longitudinal component of $\mathcal{P}_{\rm{tot}}$ and
$\mathcal{P}_{\rm{mag}}$, we infer that the reflection properties of the
optical and acoustic mode must be fundamentally different. In the case of the
acoustic mode, the chiral nature of the plasmon should be enhanced and show
unique (quite likely circular) features in typical SNOM-experiments such as
the ones of Ref. Hesp  _et al._ , 2021.

###  VI.2 Chiral resonance

From the definition of $X_{\rm{tot}}$, we infer that there is a diverging
regime for $\chi_{\rm{tot}}=0$. This regime seems to be necessarily realized
at the neutrality point for $\omega\to 0$, since the total Drude weight has to
vanish, $D_{\rm{tot}}=0$. However, in the d.c. limit also the chiral Drude
weight needs to vanish, $D_{\rm{chi}}=0$, again due to gauge
invariance.Stauber  _et al._ (2020a) At the neutrality point, no deflection is
thus expected even for the acoustic mode. At finite chemical potential,
though, Bloch electrons are deviated without a magnetic field as has recently
been discussed by several authors.Stauber  _et al._ (2018a, b); Bahamon  _et
al._ (2020); He  _et al._ (2020); Polshyn  _et al._ (2020); Sharpe  _et al._
(2021); Antebi  _et al._ (2022)

At finite frequencies, we expect sweet spots whenever $\chi_{\nu}(\omega)\to
0$ with $\nu=\rm{tot},\rm{mag}$. These frequencies lead to $X_{\nu}\to\infty$
which we will denominate as chiral resonances. At these frequencies, also the
plasmonic modes seem to eventually disappear, see Eqs. (33) and (34). However,
a coupling between the optical and acoustic mode will emerge and the plasmon
dispersion will then also depend on $\chi_{\rm{chi}}$.Lin  _et al._ (2020);
Margetis and Stauber (2021)

Chiral resonances also occur if $\chi_{\rm{chi}}\gg 1$. This is, e.g., the
case for a twist angle $\theta=1.3^{\circ}$ at $\hbar\omega\approx 50$meV
where $X_{\rm{mag}}\approx 10$, see left panel of Fig. 5. At this frequency,
the Poynting vector is largely enhanced at small wave numbers. Other sweet
spots may be limited to low temperatures, e.g., for twist angle
$\theta=1.2^{\circ}$ at $\hbar\omega\approx 22$meV, see central panel of Fig.
5. At these chiral resonance, we also assume a coupling between the optical
and acoustic mode.

###  VI.3 Chiral plasmons at the neutrality point

A Dirac system does not host plasmons at the neutrality point. Even though
electron-hole transitions may lead to positive and negative charge densities,
the charge response is always negative such that the RPA-condition for
plasmonic excitations can never be fulfilled.

This changes in moiré systems, where flat bands emerge. The moiré potential
that confines the electrons in the AA-stacked region then acts as restoring
force such that the electronic and hole charge density can oscillate in-phase.
From a technical point of view, this can be deduced from the highly peaked
absorption due to the flat bands as this may lead to a positive charge
response due to the Kramer-Kronig relation. As the collective motion is
composed of localized electrons, also the plasmonic bands are usually
flat.Stauber and Kohler (2016); Lewandowski and Levitov (2019); Novelli  _et
al._ (2020); Hesp  _et al._ (2021); Kuang  _et al._ (2021); Huang  _et al._
(2022)

One crucial condition for long-lived plasmons is the presence of an optical
gap which emerges in the continuum model by considering relaxation
effects.Koshino  _et al._ (2018) Now, if the absorption is sufficiently
peaked, a positive reactive part of the charge excitations can leak inside the
optical gap even though there are no nominal charges in the system. This
implies the possibility of a mode (”plasmon”) as a pole in RPA response. The
resulting response functions are shown in Fig. 5 for different twist angles
and temperatures with $\kappa=0.8$.

The features of the plasmonic excitations can be summarized as follows: (i)
Optical plasmons can exist right above the optical gap and persist for
temperatures up to $T\sim 50$K for $\theta<1.3^{\circ}$. This is similar to
the optical plasmons in flat bands with excess charge.Lewandowski and Levitov
(2019); Khaliji  _et al._ (2020) (ii) Acoustic plasmons can exist almost in
the whole optical window. Most notably, the magnetic Drude weight carries by
far the largest optical weight and we expect excitations with frequencies
larger than that of the corresponding optical plasmon for $qa\lesssim 100$. At
the chiral resonance for which $X_{\rm{mag}}(\omega)$ reaches a maximum, these
modes are characterized by a largely enhanced energy density $w\sim
X_{\rm{mag}}^{2}$ as can be deduced from the continuity equation and Eq. (39).

Let us finally highlight that both plasmon modes are intrinsically chiral
since $\chi_{\rm{chi}}$ is finite throughout the protected window. This is due
to the broken particle-hole symmetry as will be discussed in the Sec. VII.

##  VII Chiral response at the neutrality point

Chiral response in twisted bilayer graphene has been observed experimentally
in Ref. Kim  _et al._ , 2016 and is thus manifested in misaligned van-der-
Waals heterostructures. In Ref. Morell  _et al._ , 2017, it was shown that
neglecting the relative rotation of the pseudospin-orientation between the two
layers renders the chiral response. The difference in pseudospin orientation,
which is a consequence of the real space chiral symmetry, is thus responsible
for the chiral response in the non-interacting continuum model.

In this section, we will directly link the chiral response to particle-hole
symmetry and argue how a slight particle-hole asymmetry will lead to a finite
chiral response characterized by van Hove singularities. Our results should
also be interesting in view of other mechanisms causing particle-hole
breaking, such as non-local tunnelingXie and MacDonald (2021) or
Hartree(-Fock) renormalizationBultinck  _et al._ (2020); Zhang  _et al._
(2020); Lian  _et al._ (2021); Bernevig  _et al._ (2021); Xie  _et al._
(2021); Rademaker  _et al._ (2019); Seo  _et al._ (2019); González and Stauber
(2021); Xie and MacDonald (2021) of the bands.

###  VII.1 Symmetries of response functions

The continuum model displays particle-hole symmetry if the pseudo-spin
rotation is neglected
$\bm{\tau}_{\alpha}^{\gamma}\rightarrow\bm{\tau}_{\alpha}$.Moon and Koshino
(2013) This can be seen by the following anti-unitary transformation
$\mathcal{U}=\mathcal{S}\mathcal{P}\mathcal{K}$. The unitary operator
$\mathcal{S}$ reverts the sign of $k_{x}$,
$\mathcal{S}|k_{x},k_{y},\alpha,\ell\rangle=|-k_{x},k_{y},\alpha,\ell\rangle$,
the unitary operator $\mathcal{P}$ adds a $\pi$-phase to states in layer
$\ell=2$,
$\mathcal{P}|k_{x},k_{y},\alpha,2\rangle=-|k_{x},k_{y},\alpha,2\rangle$, and
the complex-conjugate $\mathcal{K}$ effectively changes the sign of $k_{y}$.
We thus have $\mathcal{U}\mathcal{H}\mathcal{U}^{-1}=-\mathcal{H}$.

We can now discuss the effect of $\mathcal{U}$ on the general response
function. For this, we suppress the index ${\bm{k}}$ and write

$\chi_{\mathcal{A}\mathcal{B}}=\sum_{n,m}\frac{n_{F}(\epsilon_{m})-n_{F}(\epsilon_{n})}{\omega+i0^{+}-\epsilon_{n}+\epsilon_{m}}\langle
m|\mathcal{A}|n\rangle\langle n|\mathcal{B}|m\rangle\;.$ (42)

Using the eigenbasis $\\{\tilde{n}\\}$ of $\mathcal{H}$, with
$|\tilde{n}\rangle=\mathcal{U}|n\rangle$ and
$\mathcal{H}|\tilde{n}\rangle=\epsilon_{{\tilde{n}}}|\tilde{n}\rangle$ where
$\epsilon_{{\tilde{n}}}=-\epsilon_{n}$, one can calculate any response as

$\chi_{\mathcal{A}\mathcal{B}}=\sum_{{\tilde{n}},{\tilde{m}}}\frac{n_{F}(\epsilon_{{\tilde{m}}})-n_{F}(\epsilon_{{\tilde{n}}})}{\omega+i0^{+}-\epsilon_{{\tilde{n}}}+\epsilon_{{\tilde{m}}}}\langle{\tilde{m}}|\mathcal{A}|{\tilde{n}}\rangle\langle{\tilde{n}}|\mathcal{B}|{\tilde{m}}\rangle\;.$
(43)

We can then write

$\displaystyle\frac{n_{F}(\epsilon_{{\tilde{m}}}-\mu)-n_{F}(\epsilon_{{\tilde{n}}}-\mu)}{\omega+i0^{+}-\epsilon_{{\tilde{n}}}+\epsilon_{{\tilde{m}}}}=\frac{n_{F}(\epsilon_{n}+\mu)-n_{F}(\epsilon_{m}+\mu)}{\omega+i0^{+}-\epsilon_{m}+\epsilon_{n}}\;,$
(44)

where we have explicitly included the chemical potential $\mu$ in the argument
of the Fermi function. We now have for the antiunitary transformation
$\langle{\tilde{n}}|\phi\rangle=[\langle
n|(\mathcal{U}^{\dagger}|\phi\rangle]^{*}$. Therefore, we have
$\langle{\tilde{m}}|\mathcal{A}|{\tilde{n}}\rangle=\langle
m|\tilde{\mathcal{A}}|n\rangle^{*}=\langle n|\tilde{\mathcal{A}}|m\rangle$
with $\tilde{\mathcal{A}}$ defined below. The particle-hole symmetry
$\mathcal{U}$ thus leads to the following relation:

$\displaystyle\chi_{\mathcal{A}\mathcal{B}}(\mu)=\chi_{\tilde{\mathcal{A}}\tilde{\mathcal{B}}}(-\mu)\;,$
(45)

with $\tilde{\mathcal{A}}=\mathcal{U}\mathcal{A}\mathcal{U}^{-1}$ and
$\tilde{\mathcal{B}}=\mathcal{U}\mathcal{B}\mathcal{U}^{-1}$. We now see,
because of ${\bm{\tau}}_{x}=\tilde{\bm{\tau}}_{x}$ and
${\bm{\tau}}_{y}=-\tilde{\bm{\tau}}_{y}$, that the response obeys the
following relations:

$\displaystyle\sigma_{0}(\mu)$ $\displaystyle=\sigma_{0}(-\mu)$ (46)
$\displaystyle\;\sigma_{1}(\mu)$ $\displaystyle=\sigma_{1}(-\mu)$ (47)
$\displaystyle\;\sigma_{xy}(\mu)$ $\displaystyle=-\sigma_{xy}(-\mu)\;.$ (48)

For $\mu=0$, we thus have $\sigma_{xy}=0$ for all temperatures and frequencies
as claimed.

###  VII.2 Electron and hole transitions

To make the discussion more illustrative, we switch to the particle-hole
picture by defining $\epsilon_{n}^{e}=\epsilon_{n}$ if $\epsilon_{n}>0$ and
$\epsilon_{n}^{h}=-\epsilon_{n}$ if $\epsilon_{n,\bm{k}}<0$. We only consider
vertical transitions and a general transition $n\to m$ at half-filling with
$\mu=0$ is now characterized by the initial and final energies,
$\epsilon_{n}^{h}\to\epsilon_{m}^{e}$.

For the electron-hole symmetric model, there are transitions with
$\epsilon_{n}^{h}=\epsilon_{n}^{e}$. However, this symmetry is usually
slightly broken and generally one finds
$\epsilon_{n}^{h}\neq\epsilon_{m}^{e}$. We can thus classify all (relevant)
transitions by either electron transitions if
$\epsilon_{n}^{e}>\epsilon_{m}^{h}$ or by hole transitions if
$\epsilon_{n}^{e}<\epsilon_{m}^{h}$.

Let us now denote response functions consisting of electronic (hole)
transitions as $\chi^{e(h)}$. The particle-hole transformation $\mathcal{U}$
further relates
$\epsilon_{{\tilde{n}}}^{e}=\epsilon_{{\tilde{n}}}=-\epsilon_{n}=\epsilon_{n}^{h}$
and
$\epsilon_{{\tilde{m}}}^{h}=-\epsilon_{{\tilde{m}}}=\epsilon_{m}=\epsilon_{m}^{e}$.
We now see, because of ${\bm{\tau}}_{x}=\tilde{\bm{\tau}}_{x}$ and
${\bm{\tau}}_{y}=-\tilde{\bm{\tau}}_{y}$, that the response of electron
transitions and hole transitions obeys the following relations:

$\displaystyle\chi_{xx}^{e}=\chi_{xx}^{h}\;,\;\chi_{xy}^{e}=-\chi_{xy}^{h}$
(49)

Figure 6:  Left panel: Illustration of the detailed balance relation of a
particle-hole symmtric mode. Via the anti-unitary transformation
$\mathcal{U}$, the transitions from $\epsilon_{n}\to\epsilon_{m}$ at momentum
$(k_{x},k_{y})$ are directly related to the transitions from
$\epsilon_{\tilde{m}}\to\epsilon_{\tilde{n}}$ at momentum $(-k_{x},k_{y})$.
Any hole transition ($\epsilon_{n}^{h}>\epsilon_{m}^{e}$) is automatically
related to an electron transition
($\epsilon_{{\tilde{n}}}^{e}>\epsilon_{{\tilde{m}}}^{h}$) since
$\epsilon_{{\tilde{n}}}^{e}=\epsilon_{n}^{h}$ and
$\epsilon_{{\tilde{m}}}^{h}=\epsilon_{m}^{e}$. Center panel: Chiral response
Re$\sigma_{\rm{chi}}(\omega)$ of the asymmetric continuum model with
$\kappa=0.8$ in Eq. (II) at the neutrality point with twist angle
$\theta=1.1^{\circ}$ for temperatures $T=0,10,300$K. The inset highlights the
chiral response around $\hbar\omega=95$meV. Right panel: Corresponding band
structure and density of states (DOS) on logarithmic scale. The transitions
related to the van Hove singularities around $\hbar\omega=25$meV and
$\hbar\omega=95$meV are indicated by red (electronic transition) and blue
(hole-like transition) arrows.

Numerically, we find that the dominant chiral electron (hole) transitions
between different bands and with small energy denominator are negative
(positive). However, for larger energy denominators, we also find chiral
electronic (hole) transitions which have the opposite sign. Furthermore, the
sign of the chiral response due to electron (hole) transitions between the
same bands can change. The momenta of electron and hole transitions then
normally form a well-defined boundary in the Brillouin-zone. For transitions
within the flat bands, however, we also found fractal boundaries.

###  VII.3 Detailed balance

The transformation $\mathcal{U}$ links the momentum $(k_{x},k_{y})$ to
momentum $(-k_{x},k_{y})$. Eq. (44) guarantees that the transition $n\to m$ at
momentum $(k_{x},k_{y})$ from $\epsilon_{n}$ to $\epsilon_{m}$ and at chemical
potential $-\mu$ carries the same weight as the transition
$\tilde{m}\to\tilde{n}$ at momentum $-(k_{x},k_{y})$ from
$\epsilon_{\tilde{m}}$ to $\epsilon_{\tilde{n}}$ and at chemical potential
$\mu$. Since also the matrix elements have the same (absolute) value, we thus
obtain a detailed balance relation for the above transitions at the neutrality
point $\mu=0$. This is illustrated in the left panel of Fig. 6.

With
$\chi_{\alpha\beta}=\sum_{m,n;k_{x},k_{y}}\chi_{\alpha\beta}(m,n;k_{x},k_{y})$,
we can link a single electron transition to a single hole transition as
follows:

$\displaystyle\chi_{xx}(m,n;k_{x},k_{y})$
$\displaystyle=\chi_{xx}(n,m;-k_{x},k_{y})\;,$ (50)
$\displaystyle\chi_{xy}(m,n;k_{x},k_{y})$
$\displaystyle=-\chi_{xy}(n,m;-k_{x},k_{y})\;.$ (51)

This detailed balance between the electron transition at $(k_{x},k_{y})$ and
the corresponding hole transition at $(-k_{x},k_{y})$ eventually leads to a
vanishing chiral response at half-filling.

We can also define a different particle-hole transition as was proposed by
Moon and Koshino.Moon and Koshino (2013) Together with time-reversal and
rotational symmetry, this leads to

$\displaystyle\chi_{xx}(m,n;k_{x},k_{y})$
$\displaystyle=\chi_{xx}(n,m;-k_{x},-k_{y})\;,$ (52)
$\displaystyle\chi_{xy}(m,n;k_{x},k_{y})$
$\displaystyle=-\chi_{xy}(n,m;-k_{x},-k_{y})\;.$ (53)

###  VII.4 Dissipative chiral response close to the magic angle

We will now discuss the chiral response of the full model of Eq. (II) at the
neutrality point. Crucially, the rotation in pseudospin-space needs to be
included to break particle-hole symmetry as discussed before. However, the
approximate electron-hole symmetry suggested by $\mathcal{U}$ will still
relate sublattice and layer, leading to a coherence of the wave function
between these two degrees of freedom which must not be related to the
underlying lattice (spatial) symmetry.Stauber  _et al._ (2020b); Ochoa and
Asenjo-Garcia (2020)

Since electron-hole symmetry is slightly broken, we can label all transitions
as either electron or hole transitions. The electronic wave function is not
strongly affected by this small perturbation and due to continuity arguments,
around certain regions in ${\bm{k}}$-space, electron and hole transitions must
still have well-defined, but opposite signs.

Apart from the transition-matrix element, the response is also determined by
the electronic dispersion. In any Bloch-band, there is at least one van Hove
singularity and in principle, we expect an enhanced optical response if either
the initial or final state is located at one singular ${\bm{k}}$-point.
However, the transition-matrix element might be suppressed due to symmetries
and precisely the approximate particle-hole symmetry suppresses the optical
transitions of the total current at the $M$-point.Moon and Koshino (2013) This
is not the case, though, for the magnetic and chiral transitions and we thus
expect a large response due to the large van Hove singularity which can also
be located around the $K$ or $\Gamma$-point.

In the electron-hole symmetric model, van Hove singularities necessarily
appear in the occupied and unoccupied bands at
$\epsilon_{vH}^{h}=\epsilon_{vH}^{e}$. Slightly breaking this symmetry will
lead to a splitting with $\epsilon_{vH}^{h}\neq\epsilon_{vH}^{e}$. Possible
transitions are now of electron and hole nature that have opposite chiral
response, but do not cancel each other anymore. Also the band-edges of the
electronic and hole bands will slightly shift due to the broken symmetry,
given rise to either pure electron or hole transitions. To conclude, we expect
prominent features coming from singularities of the band structure, either
discontinuities or logarithmic divergencies, where the electronic and hole
transitions are not compensated by each other.

This can be seen in the center panel of Fig. 6 where the dissipative response
of twisted bilayer with twist angle $\theta=1.1^{\circ}$ and $\kappa=0.8$ is
shown. There are always two peaks that come in pairs, a negative peak and a
positive peak associated with either electron or hole transitions.

The first pair originates from transitions within the flat bands and is
strongly temperature dependent, i.e., practicable absent at room temperature.
The second and third pair are related to transitions from the flat to the
first remote band and associated to van Hove singularities located at the
$\Gamma$ and $K$-point, respectively. They thus do not as strongly depend on
temperature and in both cases, the negative (positive) response is related to
electron (hole) transitions. The response of the third pair is highlighted in
the inset of the center panel of Fig. 6 for the sake of clarity.

In the right panel of Fig. 6, the band structure is shown and the electron
(red arrow) and hole (blue arrow) transitions are shown for the second and
third pairs. Generally, we expect strong chiral response at energies involving
a large density of states. These energies can be identified from the density-
of-states (DOS), shown next to the band structure. However, the larger the
transition energy becomes, the weaker the response is.

##  VIII Summary and Outlook

Technically speaking, we have investigated the full optical response of magic
angle graphene at the neutrality point consisting of the total, magnetic and
chiral response. The dissipative response is obtained without the usual
damping term by analytically integrating the delta-function on a linearized
grid. The reactive response is then obtained via the Kramers-Kronig relation
applying a suitable cutoff for large frequencies. By this, we obtain accurate
results close to the magic angle even for low energies.

Generally speaking, we have investigated the continuum model introduced in
Refs. Lopes dos Santos  _et al._ , 2007; Bistritzer and MacDonald, 2011 which
resembles the standard model to address general topics related to van-der-
Waals heterostructures. This model is believed to be well-understood, but here
we showed that the ground-state of the non-interacting continuum model at the
neutrality point is unstable in the immediate vicinity of the magic angle with
respect to transverse current fluctuations. We thus predict a so-called Condon
instabilityGuerci  _et al._ (2021) using a novel scaling approach.

The Condon instability at the magic angle is supposedly interesting only from
a theoretical point of view. However, we also presented new results with high
potential for technological impact. We pointed out that the plasmonic bonding
mode (acoustic or magnetic plasmon) should be larger in energy than the
ordinary plasmonic anti-bonding mode (optical or electric plasmon).
Furthermore, the energy density of this acoustic mode can be largely enhanced
at a certain frequency which we label as chiral resonance. This novel
resonance has not been discussed in the literature so far and should be
present for a wide range of twist angles and temperatures.

Another interesting aspect concerns an effective model to describe the physics
around the magic angle,Hejazi  _et al._ (2019) initially proposed in Refs.
Bena and Simon, 2011; Montambaux, 2012 in a different context. This model
makes use of an effective parameter that stands for the twist angle, and we
now provided a direct mapping to the standard continuum model of twisted
bilayer graphene, i.e., to the real twist angle. We also included a momentum-
dependent mass-term that makes sure that the universal conductivity of
$\frac{e^{2}}{2\hbar}$ is reached for $\omega\to 0$.

Lastly, we discussed the chiral aspects of the continuum model and outlined in
detail the implications of an approximate particle-hole symmetry. We
distinguished between electron and hole transitions that give equal
contributions to the chiral response, but which cancel exactly. Since
particle-hole symmetry is generally broken, we show that the finite chiral
response usually comes in pairs consisting of a positive and negative signal
since electron and hole transitions have opposite chirality, respectively.

To conclude, we hope that our results on the Condon instability will stimulate
new analytical studies of the continuum model at the magic angle regime. We
further hope that our results on the acoustic plasmonic excitations with its
chiral features will stimulate experiments which pave the way towards
technological use of this phenomenon.

##  IX Acknowledgments

This work was supported by the mobility program Salvador Madariaga under
PRX19/00024 and by the projects No. PGC2018-096955-B-C42, No.
PID2020-113164GB-i00, and No. CEX2018-000805-M financed by MCIN/
AEI/10.13039/501100011033. The access to computational resources of CESGA
(Centro de Supercomputación de Galicia) is also gratefully acknowledged. The
work of T.S. and of J.S. was further supported by Deutsche
Forschungsgemeinschaft via SFB 1277. D.M. wishes to thank Dr. A. B. Watson for
an inspiring discussion on the twisted bilayer graphene near the magic angle.

##  Appendix A Numerical integration of a generalized density of states

In this appendix, we describe the numerical recipe how to obtain the optical
response functions without introducing the usual damping term. The main
numerical task in our approach is the numerical evaluation of two-dimensional
integrals that involve a delta-function. If we determine this integral up to
large frequencies $\omega$, we can take advantage of the Kramers-Kronig
relation in order to obtain the reactive part of the response function. The
remaining one-dimensional integral over frequencies does usually not pose any
difficulties and the recipe concerning the cut-off procedure has been outlined
in Ref. Stauber  _et al._ , 2013.

We will calculate the response without disorder, i.e., we will take the delta-
function literally and perform the integration analytically after having
discretized the Brillouin zone’s . This can be done by introducing a
triangular grid on the Brillouin zone and assuming a linear interpolation. In
Fig. 7, we show the discretization used for the calculations. We have checked
that the final result does not crucially depend on the discretization. Another
optimization is obtained by assuming a quadratic interpolation between the
three base-points.Wiesenekker  _et al._ (1988) We have checked that for large
grids used here, this does not lead to significant improvements, also nicely
explained in Ref. Pedersen  _et al._ , 2008.

We discretize the Brillouin zone by $N$ with $n,m=0,...,N$ in the following
way:

$\displaystyle{\bm{k}}=\frac{n}{N}{\bm{G}}_{1}+\frac{m}{N}{\bm{G}}_{2}\;,$
(54)

with the lattice vectors
$\bm{G}_{1}=|\Delta\bm{K}|\left(\tfrac{-\sqrt{3}}{2},-\tfrac{3}{2}\right)$,
$\bm{G}_{2}=|\Delta\bm{K}|\left(\tfrac{\sqrt{3}}{2},-\tfrac{3}{2}\right)$, see
Fig. 7 A). In our calculations, we chose discretizations up to $N\cong 10000$;
for twist angles in the immediate vicinity of the magic angle even as large as
$N\cong 20000$.

We shall calculate the following generalized density of states with $g$
denoting a degeneracy factor:

$\displaystyle\rho(\epsilon)=\frac{g}{\mathcal{A}}\sum_{\bm{k}}f_{\bm{k}}\delta(\epsilon-\epsilon_{\bm{k}})$
(55)

As we assume periodic boundary conditions, the sample area is given by
$\mathcal{A}=N^{2}A_{c}$ where $A_{c}$ is the area of the unit cell. In the
case of twisted bilayer graphene, we have
$A_{c}=\frac{\sqrt{3}}{2}a_{g}^{2}A_{i}$ as the area of the moiré supercell
with $a_{g}=2.46\,\mathring{\text{A}}$, $A_{i}=3i^{2}+3i+1$ and
$\cos(\theta_{i})=1-\tfrac{1}{2A_{i}}$.

We will now consider each of the $N^{2}$ mini-rhombi individually which are
characterized by the vertices ${\bm{k}}_{i}$, $\epsilon_{i}$, and optionally
$f_{i}$ with $i=A,B,C,D$. First, we will divide the mini-rombus in two and
consider first the triangle defined by $i=A,B,C$ and afterwards the triangle
defined by $i=B,C,D$.

To outline the algorithm, we will only consider the first triangle and further
assume that $\epsilon_{A}\leq\epsilon_{B}\leq\epsilon_{C}$ which can always be
achieved by relabelling the vertices. We now interpolate linearly between the
three vertices such that any momentum ${\bm{k}}$ and energy $\epsilon$ inside
the triangle can be parameterized by two parameters $t,s\in[0,1]$ (due to the
prior ordering):

$\displaystyle\left(\begin{matrix}k_{x}\\\ k_{y}\\\
\epsilon_{\bm{k}}\end{matrix}\right)=\left(\begin{matrix}k_{B,x}-k_{A,x}\\\
k_{B,y}-k_{A,y}\\\
\epsilon_{{\bm{k}}_{B}}-\epsilon_{{\bm{k}}_{A}}\end{matrix}\right)t+\left(\begin{matrix}k_{C,x}-k_{A,x}\\\
k_{C,y}-k_{A,y}\\\
\epsilon_{{\bm{k}}_{C}}-\epsilon_{{\bm{k}}_{A}}\end{matrix}\right)s+\left(\begin{matrix}k_{A,x}\\\
k_{A,y}\\\ \epsilon_{{\bm{k}}_{A}}\end{matrix}\right)\;,$ (56)

We can now write the integral that contains the contribution $\rho_{\Delta}$
to $\rho$ over the triangle with respect to the two variables $t$ and $s$. The
integration limits corresponding to the vertices $[A,B,C]$ are now given with
respect to the axis defined by $t,s$, i.e., $[(0,0),(1,0),(0,1)]$. Neglecting
for the moment the weight function $f_{\bm{k}}$ and setting $g=1$, we arrive
at the following expression:

$\displaystyle\rho_{\triangle}(\epsilon)$
$\displaystyle=\frac{1}{(2\pi)^{2}}\int_{\triangle}d^{2}k\;\delta(\epsilon-\epsilon_{\bm{k}})\;,$
(57)
$\displaystyle=\frac{\tilde{J}}{(2\pi)^{2}}\int_{0}^{1}dt\int_{0}^{1-t}ds\;\delta(s-s(\epsilon,t))\;,$
(58)

where we introduced the Jacobian
$J=|(k_{B,x}-k_{A,x})(k_{C,y}-k_{A,y})-(k_{C,x}-k_{A,x})(k_{B,y}-k_{A,y})|$
with $\tilde{J}=J/(\epsilon_{C}-\epsilon_{A})$ and
$s(\epsilon,t)=\frac{\epsilon-\epsilon_{A}}{\epsilon_{C}-\epsilon_{A}}-\frac{\epsilon_{B}-\epsilon_{A}}{\epsilon_{C}-\epsilon_{A}}t$.

The integral depends on the value of $\epsilon$ relative to the energies
$\epsilon_{i}$ and we obtain

$\displaystyle\rho_{\Delta}(\epsilon)=\frac{\tilde{J}}{(2\pi)^{2}}\Big{[}$
$\displaystyle\frac{\epsilon-\epsilon_{A}}{\epsilon_{B}-\epsilon_{A}}\theta(\epsilon-\epsilon_{A})\theta(\epsilon_{B}-\epsilon)$
$\displaystyle+$
$\displaystyle\frac{\epsilon_{C}-\epsilon}{\epsilon_{C}-\epsilon_{B}}\theta(\epsilon-\epsilon_{B})\theta(\epsilon_{C}-\epsilon)\Big{]}\;.$
(59)

The total density of states is then obtained by the sum
$\rho(\epsilon)=\sum_{\Delta}\rho_{\Delta}(\epsilon)$.

The weight function $f_{\bm{k}}$ can now be included by linear interpolation.
With

$\displaystyle f(\epsilon)=$ $\displaystyle
f_{B}\frac{\epsilon-\epsilon_{A}}{\epsilon_{B}-\epsilon_{A}}\theta(\epsilon-\epsilon_{A})\theta(\epsilon_{B}-\epsilon)$
$\displaystyle+$ $\displaystyle
f_{A}\frac{\epsilon_{B}-\epsilon}{\epsilon_{B}-\epsilon_{A}}\theta(\epsilon-\epsilon_{A})\theta(\epsilon_{B}-\epsilon)$
$\displaystyle+$ $\displaystyle
f_{C}\frac{\epsilon-\epsilon_{B}}{\epsilon_{C}-\epsilon_{B}}\theta(\epsilon-\epsilon_{B})\theta(\epsilon_{C}-\epsilon)$
$\displaystyle+$ $\displaystyle
f_{B}\frac{\epsilon_{C}-\epsilon}{\epsilon_{C}-\epsilon_{B}}\theta(\epsilon-\epsilon_{B})\theta(\epsilon_{C}-\epsilon)\;,$
(60)

and reincorporation of the degeneracy factor $g$, the generalized density of
states is thus approximated by

$\displaystyle\rho_{\Delta}(\epsilon)=\frac{g\tilde{J}}{(2\pi)^{2}}f(\epsilon)\Big{[}$
$\displaystyle\frac{\epsilon-\epsilon_{A}}{\epsilon_{B}-\epsilon_{A}}\theta(\epsilon-\epsilon_{A})\theta(\epsilon_{B}-\epsilon)$
$\displaystyle+$
$\displaystyle\frac{\epsilon_{C}-\epsilon}{\epsilon_{C}-\epsilon_{B}}\theta(\epsilon-\epsilon_{B})\theta(\epsilon_{C}-\epsilon)\Big{]}\;.$
(61)

Apart from increasing the discretization, the numerical results can be further
smoothened by explicitly taking advantage of the rotational symmetry, i.e.,
$2{\bm{j}}_{\rm{tot}}\cdot{\bm{j}}_{\rm{tot}}=(j_{x}^{1}+j_{x}^{2})^{2}+(j_{y}^{1}+j_{y}^{2})^{2}$,
$2{\bm{j}}_{\rm{mag}}\cdot{\bm{j}}_{\rm{mag}}=(j_{x}^{1}-j_{x}^{2})^{2}+(j_{y}^{1}-j_{y}^{2})^{2}$,
and $2{\bm{j}}_{xy}\cdot{\bm{j}}_{xy}=j_{x}^{1}j_{y}^{2}-j_{x}^{2}j_{y}^{1}$.

Figure 7:  (A) The rhombic Brillouin zone defined by the reciprocal lattice
vectors ${\bm{G}}_{1}$ and ${\bm{G}}_{2}$. (B) Zoom-In of a small rhombus with
side length $|{\bm{G}}_{1}|/N$. The vertices are labeled by $i=A,B,C,D$ and
characterized by ${\bm{k}}_{i}$, $\epsilon_{i}$, and optionally also by
$f_{i}$.

##  Appendix B Real part of interband conductivity: Analytical derivations

In this appendix, we describe analytically the real part of the interband
conductivity for $\mu=0$ and $T=0$, by use of the two-band model introduced in
Refs. Bena and Simon, 2011; Montambaux, 2012. In this model, the Dirac cone
coexists with a parabolic profile in the Hamiltonian. We focus on the limits
of the interband conductivity as $\omega\to 0$ and $\omega\to\infty$.

###  B.1 Model Hamiltonian

The reduced, two-band Hamiltonian without a gap reads

$H_{red}=-\begin{pmatrix}0&\varpi^{2}+\eta\varpi^{*}\\\
{\varpi^{*}}^{2}+\eta\varpi&0\end{pmatrix}~{},$ (62)

where $\varpi=k_{x}\mp ik_{y}$ in the vicinity of $\bm{K}$
($\bm{K}^{\prime}$). Here, we have set $\hbar=1=2m$ for later algebraic
convenience. The parameter $\eta$ is assumed positive and small ($0<\eta\ll
1$). It expresses the relative strength of the Dirac cone. From now on, we
focus on the point $\bm{K}$. We will comment on the case with $\eta<0$ below.

This Hamiltonian yields the eigenenergies

$\epsilon_{\bm{k},\pm}=\pm|F_{\bm{k}}|~{},\quad
F_{\bm{k}}={\varpi^{*}}^{2}+\eta\varpi~{},$ (63)

and the normalized eigenvectors

$|\pm\rangle_{\bm{k}}=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\ \mp
e^{i\vartheta_{\bm{k}}}\end{pmatrix}~{},\quad\vartheta_{\bm{k}}={\rm
Arg}F_{\bm{k}}~{}.$ (64)

The eigenenergies are expressed explicitly by

$\epsilon_{\bm{k},\pm}=\pm\sqrt{k^{4}+2\eta
k^{3}\cos(3\theta)+\eta^{2}k^{2}}=\pm\epsilon(\bm{k};\eta)~{},$ (65)

in the polar coordinates $(k,\theta)$ with center at $\bm{K}$.

Evidently, the scaling of the momentum with $\eta$ according to
$\bm{k}=\eta\tilde{\bm{k}}$ results in
$\epsilon(\bm{k};\eta)=\eta^{2}\tilde{\epsilon}(\tilde{\bm{k}})$ where

$\tilde{\epsilon}(\bm{k})=\sqrt{k^{4}+2k^{3}\cos(3\theta)+k^{2}}=\epsilon(\bm{k};1)~{}.$
(66)

It is algebraically convenient to use the scaled momentum and eigenenergy
(see, however, Eq. (68)). For ease of notation, we henceforth drop the tildes
from $\tilde{\bm{k}}$ and $\tilde{\epsilon}$.

Next, we describe the local minima of $\epsilon_{\bm{k},+}=\epsilon(\bm{k})$.
By $\nabla_{\bm{k}}(\epsilon^{2})=0$ we obtain
$\bm{k}=k_{c}(\cos\theta_{c},\sin\theta_{c})$ where $k_{c}=0$, or $k_{c}=1$
with $\theta_{c}=\pi-2\pi n/3$, $n\in\mathbb{N}$. These points yield zero
bandgap. The other critical points of $\epsilon(\bm{k})$ correspond to saddle
points, with nonzero bandgap, and are disregarded. If $\eta<0$, the local
minima correspond to $k_{c}=0$, or $k_{c}=1$ with $\theta_{c}=2\pi n/3$ (by
$k_{x}\rightarrow-k_{x}$).

We turn our attention to the velocity matrix element needed for the interband
conductivity. By setting $F=F_{R}+iF_{I}$ ($F_{R}={\rm{Re}}F$ and
$F_{I}={\rm{Im}}F$), we have

$\displaystyle\langle-|\nabla_{\bm{k}}H_{red}|+\rangle$
$\displaystyle=-i\epsilon(\bm{k})\,\nabla_{\bm{k}}\vartheta_{\bm{k}}$
$\displaystyle=\frac{i}{\epsilon(\bm{k})}\\{(\nabla_{\bm{k}}F_{R})F_{I}-F_{R}(\nabla_{\bm{k}}F_{I})\\}~{}.$

Let’s compute the $x$-component, for example. We find

$\langle-|\partial_{k_{x}}H_{red}|+\rangle=~{}-2i\,\frac{\Lambda_{x}(\bm{k})}{\epsilon(\bm{k})}$

where

$\Lambda_{x}(\bm{k})=k_{y}\left\\{\left(k_{x}-\frac{1}{2}\right)^{2}+k_{y}^{2}-\frac{3}{4}\right\\}~{}.$
(67)

For $\bm{K}^{\prime}$, one simply has to replace $k_{x}$ by $-k_{x}$. Note
that $\Lambda_{x}(\bm{k})=0$ at the local minima of $\epsilon(\bm{k})$
determined above.

###  B.2 Integral of interband conductivity

The diagonal elements of the interband (regular) conductivity are computed
from the formula ($\alpha=x,\,y$)

$\sigma_{\alpha\alpha}^{R}(\omega)=4ig_{s}g_{v}g_{\ell}\sigma_{G}\
(\omega+i0^{+})\int\limits\frac{d^{2}{\bm{k}}}{(2\pi)^{2}}\
\frac{n_{F}(\epsilon(\bm{k};\eta))-n_{F}(-\epsilon(\bm{k};\eta))}{\epsilon(\bm{k};\eta)}\frac{|\langle-|\partial_{k_{\alpha}}H_{red}|+\rangle|^{2}}{4\epsilon(\bm{k};\eta)^{2}-(\omega+i0^{+})^{2}}~{};\quad\sigma_{G}=\frac{e^{2}}{4\hbar}.$
(68)

Here, the factor $g_{s}g_{v}g_{\ell}$ accounts for the layer-degree of
freedom, and the usual spin and valley degeneracies. In Eq. (68) we use the
unscaled momentum $\bm{k}$ and the eigenenergy $\epsilon(\bm{k};\eta)$ from
Eq. (65). We set $\mu=0$ and $T=0$, take $\alpha=x$, and change the
integration variable from $\bm{k}$ to $\eta\bm{k}$. Thus, we arrive at the
simplified integral

$\sigma_{xx}^{R}(\omega)=-8ig_{s}g_{v}g_{\ell}\sigma_{G}\tilde{\omega}\int\limits\frac{d^{2}{\bm{k}}}{(2\pi)^{2}}\
\frac{1}{\epsilon(\bm{k})^{3}}\frac{\Lambda_{x}(\bm{k})^{2}}{\epsilon(\bm{k})^{2}-(\tilde{\omega}+i0^{+})^{2}}$
(69)

where $\epsilon(\bm{k})$ is given by Eq. (66) and

$\tilde{\omega}=\frac{\omega}{2\eta^{2}}~{}.$ (70)

We will keep the symbol $\tilde{\omega}$ (with tilde) throughout.

Our task is to compute ${\rm{Re}}\sigma_{xx}^{R}$ by carrying out the
integration in _local_ polar coordinates by consideration of points
$\bm{k}=\bm{k}_{*}$ such that $\epsilon(\bm{k})=\tilde{\omega}$ (if
$\tilde{\omega}>0$). A difficulty is that these points may locally form non-
circular curves. The integral for ${\rm{Re}}\sigma_{xx}^{R}$ has significant
contributions from the vicinity of each curve. We study the following limits:
(i) $\tilde{\omega}\to 0$, when the curves of interest are formed near local
minima of $\epsilon(\bm{k})$; and (ii) $\tilde{\omega}\to+\infty$, when
$k_{*}$ is large.

####  B.2.1 Limit $\tilde{\omega}\to 0$

For each critical point of interest we set
$\bm{q}=\bm{k}-k_{c}(\cos\theta_{c},\sin\theta_{c})$, and find a suitable
expansion for the solutions $\bm{q}$ of $\epsilon(\bm{k})=\tilde{\omega}$ by
perturbations if $\tilde{\omega}\ll 1$. For this purpose, we invoke the local
polar coordinates $(q,\phi)$, where $\bm{q}=q(\cos\phi,\sin\phi)$
($0\leq\phi<2\pi$); and determine $q=|\bm{q}|$ as a function of $\phi$ and
$\tilde{\omega}$. Let $q_{*}(\phi)$ be such a solution. Subsequently, we
expand $\epsilon(\bm{k})^{2}$ near $q=q_{*}$.

First, consider $k_{c}=0$, which amounts to the center point ($\bm{K}$). After
some algebra, we obtain

$q_{*}=\tilde{\omega}\left\\{1-\tilde{\omega}\cos(3\phi)+O(\tilde{\omega}^{2})\right\\}~{},$
(71)

where $O(\tilde{\omega}^{2})$ denotes a correction of the order of
$\tilde{\omega}^{2}$. This formula entails an approximation of the form

$\displaystyle\epsilon(\bm{k})^{2}-\tilde{\omega}^{2}$ $\displaystyle\simeq
Q_{1}(\phi)\,(q-q_{*})+Q_{2}(\phi)\,(q-q_{*})^{2}$ (72)

where $Q_{1}(\phi)=2q_{*}[1+3q_{*}\cos(3\phi)]$ and $Q_{2}(\phi)=1$.

Second, consider $k_{c}=1$ with $\theta_{c}=\pi$, which amounts to the
critical point at $\bm{k}=(-1,0)$, for $n=0$. We find

$q_{*}\simeq\frac{\tilde{\omega}}{\sqrt{1+8\sin^{2}\phi}}\left\\{1+\tilde{\omega}\frac{\cos\phi(1+4\sin^{2}\phi)}{(1+8\sin^{2}\phi)^{\frac{3}{2}}}\right\\}~{}.$
(73)

This formula implies expansion (72) with

$\displaystyle
Q_{1}(\phi)=2q_{*}\left[(1+8\sin^{2}\phi)+3q_{*}(1+4\sin^{2}\phi)\right]$ (74)

and $Q_{2}(\phi)=1+8\sin^{2}\phi$.

Third, we consider the critical points with $k_{c}=1$ and $\theta_{c}=\pi-2\pi
n/3$ for $n=1,\,2$, i.e., at $\bm{k}=(1/2,\pm\sqrt{3}/2)$. We thus obtain the
following expansions for $q_{*}=q_{*}(\phi)$:

$\displaystyle q_{*}$
$\displaystyle\simeq\frac{\tilde{\omega}}{\sqrt{4\cos^{2}\phi+3\mp
2\sqrt{3}\sin(2\phi)}}\left\\{1-\tilde{\omega}\right.$
$\displaystyle\qquad\left.\times\frac{2\cos\phi\cos(2\phi)\pm\sqrt{3}\sin\phi}{[4\cos^{2}\phi+3\mp
2\sqrt{3}\sin(2\phi)]^{\frac{3}{2}}}\right\\}~{}.$ (75)

Each of these formulas implies expansion (72) with

$Q_{1}(\phi)=2q_{*}\left\\{4\cos^{2}\phi+3\mp
2\sqrt{3}\sin(2\phi)+3q_{*}[2\cos\phi\cos(2\phi)\pm\sqrt{3}\sin\phi]+2q_{*}^{2}\right\\}$
(76)

and $Q_{2}(\phi)=4\cos^{2}\phi+3\mp 2\sqrt{3}\sin(2\phi)$.

In all of the above cases, we have $Q_{1}(\phi)\neq 0$ for every $\phi$. The
expansions for $\epsilon(\bm{k})$ near local minima are uniform in $\phi$; and
capture the zero bandgap with a negligible correction of the order of
$\tilde{\omega}^{3/2}$ or smaller. This property can be used to show (as a
self-consistency check) that our leading-order result for
${\rm{Re}}\sigma_{xx}^{R}$, given below, has a negligible correction if
$\tilde{\omega}\ll 1$. We omit details on this here.

Next, by Eq. (69), we split the integral for $\sigma_{xx}^{R}(\omega)$ into
four contributions, one for each local minimum of $\epsilon(\bm{k})$. Using
the local polar coordinates $(q,\phi)$, we first carry out the integration in
$q$ by employing the formula

$\frac{1}{\epsilon(\bm{k})^{2}-(\tilde{\omega}+i0^{+})^{2}}=i\pi
Q_{1}(\phi)^{-1}\delta(q-q_{*}(\phi))+\mathcal{P}\left(\frac{1}{q-q_{*}}\right)$

for each contribution. In the above, $\mathcal{P}(\cdot)$ indicates the
principal-value integral. The two-dimensional integral for
${\rm{Re}}\sigma_{xx}^{R}$ immediately reduces to an integral with respect to
the polar angle $\phi$, from the delta function term.

Accordingly, we perform the remaining integration, with respect to $\phi$. For
$\tilde{\omega}\ll 1$, we write

${\rm{Re}}\sigma_{xx}^{R}(\omega)\simeq\frac{1}{2}g_{\ell}g_{s}g_{v}\sigma_{G}\left\\{I^{(c)}+9\sum_{n=0}^{2}I^{(n)}\right\\}~{},$
(77)

where $I^{(c)}$ and $I^{(n)}$ correspond to the center point ($k_{c}=0$) and
the points $k_{c}=1$ and $\theta_{c}=\pi-2\pi n/3$, respectively. We define
and compute the following requisite integrals:

$I^{(c)}=\frac{1}{2\pi}\int_{0}^{2\pi}d\phi\,\sin^{2}\phi=\frac{1}{2}~{},$
$I^{(0)}=\frac{1}{2\pi}\int_{0}^{2\pi}d\phi\,\frac{\sin^{2}\phi}{(1+8\sin^{2}\phi)^{2}}=\frac{1}{54}~{},$
$I^{(1)}+I^{(2)}=\frac{1}{2\pi}\int_{0}^{2\pi}d\phi\,\sin^{2}\phi\sum_{s=\pm}\frac{1}{[4\cos^{2}\phi+3+s2\sqrt{3}\sin(2\phi)]^{2}}=\frac{7}{27}~{}.$

Hence, we finally obtain

${\rm{Re}}\sigma_{xx}^{R}(\omega)=\frac{3}{2}g_{\ell}g_{s}g_{v}\sigma_{G}=12\sigma_{G}\qquad\mbox{as}\
\tilde{\omega}\to 0~{},$ (78)

for the TBG system. The anticipated correction to this result for small
nonzero $\tilde{\omega}$ is of the order of $\tilde{\omega}^{2}$.

####  B.2.2 Limit $\tilde{\omega}\to\infty$

In this case, we apply a procedure similar to the above. In particular, we
solve the equation $\epsilon(\bm{k})=\tilde{\omega}$ for large
$\tilde{\omega}$, to find $k=k_{*}(\phi)\gg 1$. Then we expand the difference
$\epsilon(\bm{k})^{2}-\tilde{\omega}^{2}$ in powers of $k-k_{*}$. The quantity
${\rm{Re}}\sigma_{xx}^{R}$ is determined by integration near the curve
$k=k_{*}(\phi)$.

In detail, by perturbations for $\tilde{\omega}\gg 1$, we obtain

$k_{*}=\sqrt{\tilde{\omega}}\left\\{1+\frac{1}{2}\tilde{\omega}^{-1/2}\cos(3\theta)+O(\tilde{\omega}^{-1})\right\\}~{}.$
(79)

Here, we use the polar angle $\theta$, where
$\bm{k}=k(\cos\theta,\sin\theta)$. The above formula implies expansion (72)
with $q=k$ and $\phi=\theta$ (since $\bm{k}=\bm{q}$ here), while

$Q_{1}(\theta)=2k_{*}\left[2k_{*}^{2}-3k_{*}\cos(3\theta)+1\right]$ (80)

and $Q_{2}(\theta)=6k_{*}^{2}-6k_{*}\cos(3\theta)+1$.

These considerations lead to the simplified integral

$\displaystyle{\rm{Re}}\sigma_{xx}^{R}(\omega)$ $\displaystyle\simeq
g_{\ell}g_{s}g_{v}\sigma_{G}\tilde{\omega}^{-2}\frac{1}{2\pi}\int_{0}^{2\pi}d\theta\,k_{*}(\theta)^{4}\sin^{2}\theta$
$\displaystyle\simeq\frac{1}{2}g_{\ell}g_{s}g_{v}\sigma_{G}=4\sigma_{G}\quad\mbox{as}\
\tilde{\omega}\to+\infty~{},$ (81)

for the TBG system.

## References

  * 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).
  * Shallcross  _et al._ (2008) S. Shallcross, S. Sharma, and O. A. Pankratov, “Quantum interference at the twist boundary in graphene,” Phys. Rev. Lett. 101, 056803 (2008).
  * 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 (2010).
  * Schmidt  _et al._ (2010) H. Schmidt, T. Lüdtke, P. Barthold, and R. J. Haug, “Mobilities and scattering times in decoupled graphene monolayers,” Phys. Rev. B 81, 121403 (2010).
  * Li  _et al._ (2010) Guohong Li, A. Luican, J. M. B. Lopes dos Santos, A. H. Castro Neto, A. Reina, J. Kong, and E. Y. Andrei, “Observation of van hove singularities in twisted graphene layers,” Nat. Phys. 6, 109–113 (2010).
  * de Laissardière  _et al._ (2010) G. Trambly de Laissardière, D. Mayou, and L. Magaud, “Localization of dirac electrons in rotated graphene bilayers,” Nano Letters 10, 804–808 (2010).
  * Bistritzer and MacDonald (2011) Rafi Bistritzer and Allan H. MacDonald, “Moiré bands in twisted double-layer graphene,” P. Natl. Acad. Sci. Usa. 108, 12233–12237 (2011).
  * 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).
  * Dean  _et al._ (2013) C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, K. L. Shepard, J. Hone, and P. Kim, “Hofstadter’s butterfly and the fractal quantum hall effect in moirésuperlattices,” Nature 497, 598 – 602 (2013).
  * Cao  _et al._ (2018a) Yuan Cao, Valla Fatemi, Ahmet Demir, Shiang Fang, Spencer L. Tomarken, Jason Y. Luo, Javier D. Sanchez-Yamagishi, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, Ray C. Ashoori, and Pablo Jarillo-Herrero, “Correlated insulator behaviour at half-filling in magic-angle graphene superlattices,” Nature 556, 80 – 84 (2018a).
  * Cao  _et al._ (2018b) Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature 556, 43 – 50 (2018b).
  * Yankowitz  _et al._ (2019) Matthew Yankowitz, Shaowen Chen, Hryhoriy Polshyn, Yuxuan Zhang, K. Watanabe, T. Taniguchi, David Graf, Andrea F. Young, and Cory R. Dean, “Tuning superconductivity in twisted bilayer graphene,” Science 363, 1059–1064 (2019).
  * Sharpe  _et al._ (2019) Aaron L. Sharpe, Eli J. Fox, Arthur W. Barnard, Joe Finney, Kenji Watanabe, Takashi Taniguchi, M. A. Kastner, and David Goldhaber-Gordon, “Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene,” Science 365, 605–608 (2019).
  * Polshyn  _et al._ (2020) H. Polshyn, J. Zhu, M. A. Kumar, Y. Zhang, F. Yang, C. L. Tschirhart, M. Serlin, K. Watanabe, T. Taniguchi, A. H. MacDonald, and A. F. Young, “Electrical switching of magnetic order in an orbital chern insulator,” Nature 588, 66–70 (2020).
  * Kim  _et al._ (2016) Cheol-Joo Kim, Sánchez-Castillo A., Zack Ziegler, Yui Ogawa, Cecilia Noguez, and Jiwoong Park, “Chiral atomically thin films,” Nat. Nanotechnol. 11, 520–524 (2016).
  * Morell  _et al._ (2017) E Suárez Morell, Leonor Chico, and Luis Brey, “Twisting dirac fermions: circular dichroism in bilayer graphene,” 2D Materials 4, 035015 (2017).
  * Stauber  _et al._ (2018a) T. Stauber, T. Low, and G. Gómez-Santos, “Chiral response of twisted bilayer graphene,” Phys. Rev. Lett. 120, 046801 (2018a).
  * Stauber  _et al._ (2018b) T. Stauber, T. Low, and G. Gómez-Santos, “Linear response of twisted bilayer graphene: Continuum versus tight-binding models,” Phys. Rev. B 98, 195414 (2018b).
  * Hu  _et al._ (2017) F. Hu, Suprem R. Das, Y. Luan, T.-F. Chung, Y. P. Chen, and Z. Fei, “Real-space imaging of the tailored plasmons in twisted bilayer graphene,” Phys. Rev. Lett. 119, 247402 (2017).
  * Stauber and Kohler (2016) Tobias Stauber and Heinerich Kohler, “Quasi-flat plasmonic bands in twisted bilayer graphene,” Nano Lett. 16, 6844–6849 (2016).
  * Lewandowski and Levitov (2019) Cyprian Lewandowski and Leonid Levitov, “Intrinsically undamped plasmon modes in narrow electron bands,” Proceedings of the National Academy of Sciences 116, 20869–20874 (2019).
  * Novelli  _et al._ (2020) Pietro Novelli, Iacopo Torre, Frank H. L. Koppens, Fabio Taddei, and Marco Polini, “Optical and plasmonic properties of twisted bilayer graphene: Impact of interlayer tunneling asymmetry and ground-state charge inhomogeneity,” Phys. Rev. B 102, 125403 (2020).
  * Hesp  _et al._ (2021) Niels C. H. Hesp, Iacopo Torre, Daniel Rodan-Legrain, Pietro Novelli, Yuan Cao, Stephen Carr, Shiang Fang, Petr Stepanov, David Barcons-Ruiz, Hanan Herzig Sheinfux, Kenji Watanabe, Takashi Taniguchi, Dmitri K. Efetov, Efthimios Kaxiras, Pablo Jarillo-Herrero, Marco Polini, and Frank H. L. Koppens, “Observation of interband collective excitations in twisted bilayer graphene,” Nature Physics 17, 1162–1168 (2021).
  * Huang  _et al._ (2022) Tianye Huang, Xuecou Tu, Changqing Shen, Binjie Zheng, Junzhuan Wang, Hao Wang, Kaveh Khaliji, Sang Hyun Park, Zhiyong Liu, Teng Yang, Zhidong Zhang, Lei Shao, Xuesong Li, Tony Low, Yi Shi, and Xiaomu Wang, “Observation of chiral and slow plasmons in twisted bilayer graphene,” Nature 605, 63–68 (2022).
  * Kuang  _et al._ (2021) Xueheng Kuang, Zhen Zhan, and Shengjun Yuan, “Collective excitations and flat-band plasmon in twisted bilayer graphene near the magic angle,” Phys. Rev. B 103, 115431 (2021).
  * Tarnopolsky  _et al._ (2019) Grigory Tarnopolsky, Alex Jura Kruchkov, and Ashvin Vishwanath, “Origin of magic angles in twisted bilayer graphene,” Phys. Rev. Lett. 122, 106405 (2019).
  * Song  _et al._ (2019) Zhida Song, Zhijun Wang, Wujun Shi, Gang Li, Chen Fang, and B. Andrei Bernevig, “All magic angles in twisted bilayer graphene are topological,” Phys. Rev. Lett. 123, 036401 (2019).
  * Hejazi  _et al._ (2019) Kasra Hejazi, Chunxiao Liu, Hassan Shapourian, Xiao Chen, and Leon Balents, “Multiple topological transitions in twisted bilayer graphene near the first magic angle,” Phys. Rev. B 99, 035111 (2019).
  * Koshino  _et al._ (2018) Mikito Koshino, Noah F. Q. Yuan, Takashi Koretsune, Masayuki Ochi, Kazuhiko Kuroki, and Liang Fu, “Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene,” Phys. Rev. X 8, 031087 (2018).
  * Park  _et al._ (2020) Youngju Park, Bheema Lingam Chittari, and Jeil Jung, “Gate-tunable topological flat bands in twisted monolayer-bilayer graphene,” Phys. Rev. B 102, 035411 (2020).
  * Stauber  _et al._ (2013) Tobias Stauber, Pablo San-Jose, and Luis Brey, “Optical conductivity, drude weight and plasmons in twisted graphene bilayers,” New J. Phys. 15, 113050 (2013).
  * Vela  _et al._ (2018) Adriana Vela, M. V. O. Moutinho, F. J. Culchac, P. Venezuela, and Rodrigo B. Capaz, “Electronic structure and optical properties of twisted multilayer graphene,” Phys. Rev. B 98, 155135 (2018).
  * Dai  _et al._ (2021) Zhen-Bing Dai, Yan He, and Zhiqiang Li, “Effects of heterostrain and lattice relaxation on the optical conductivity of twisted bilayer graphene,” Phys. Rev. B 104, 045403 (2021).
  * Han and Lai (2022) Chen-Di Han and Ying-Cheng Lai, “Optical response of two-dimensional dirac materials with a flat band,” Phys. Rev. B 105, 155405 (2022).
  * Condon and Walstedt (1968) J. H. Condon and R. E. Walstedt, “Direct evidence for magnetic domains in silver,” Phys. Rev. Lett. 21, 612–614 (1968).
  * Andolina  _et al._ (2020) G. M. Andolina, F. M. D. Pellegrino, V. Giovannetti, A. H. MacDonald, and M. Polini, “Theory of photon condensation in a spatially varying electromagnetic field,” Phys. Rev. B 102, 125137 (2020).
  * Nataf  _et al._ (2019) Pierre Nataf, Thierry Champel, Gianni Blatter, and Denis M. Basko, “Rashba cavity $\rm{QED}$: A route towards the superradiant quantum phase transition,” Phys. Rev. Lett. 123, 207402 (2019).
  * Guerci  _et al._ (2020) Daniele Guerci, Pascal Simon, and Christophe Mora, “Superradiant phase transition in electronic systems and emergent topological phases,” Phys. Rev. Lett. 125, 257604 (2020).
  * Sánchez  _et al._ (2021) M. Sánchez Sánchez, G. Gómez-Santos, and T. Stauber, “Collective magnetic excitations in $\rm{AA}$\- and $\rm{AB}$-stacked graphene bilayers,” Phys. Rev. B 104, 245412 (2021).
  * Guerci  _et al._ (2021) Daniele Guerci, Pascal Simon, and Christophe Mora, “Moiré lattice effects on the orbital magnetic response of twisted bilayer graphene and condon instability,” Phys. Rev. B 103, 224436 (2021).
  * Stauber  _et al._ (2020a) Tobias Stauber, Tony Low, and Guillermo Gómez-Santos, “Plasmon-enhanced near-field chirality in twisted van der waals heterostructures,” Nano Letters 20, 8711–8718 (2020a).
  * Moon and Koshino (2013) Pilkyung Moon and Mikito Koshino, “Optical absorption in twisted bilayer graphene,” Phys. Rev. B 87, 205404 (2013).
  * Guinea and Walet (2019) Francisco Guinea and Niels R. Walet, “Continuum models for twisted bilayer graphene: Effect of lattice deformation and hopping parameters,” Phys. Rev. B 99, 205134 (2019).
  * Margetis and Stauber (2021) Dionisios Margetis and Tobias Stauber, “Theory of plasmonic edge states in chiral bilayer systems,” Phys. Rev. B 104, 115422 (2021).
  * Wiesenekker  _et al._ (1988) G. Wiesenekker, G. te Velde, and E. J. Baerends, “Analytic quadratic integration over the two-dimensional brillouin zone,” J. Phys. C: Solid State Phys. 21, 4263–4283 (1988).
  * 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–1308 (2008).
  * Kuzmenko  _et al._ (2008) A. B. Kuzmenko, E. van Heumen, F. Carbone, and D. van der Marel, “Universal optical conductance of graphite,” Phys. Rev. Lett. 100, 117401 (2008).
  * Falkovsky and Pershoguba (2007) L. A. Falkovsky and S. S. Pershoguba, “Optical far-infrared properties of a graphene monolayer and multilayer,” Phys. Rev. B 76, 153410 (2007).
  * Stauber  _et al._ (2008) T. Stauber, N. M. R. Peres, and A. K. Geim, “Optical conductivity of graphene in the visible region of the spectrum,” Phys. Rev. B 78, 085432 (2008).
  * Li  _et al._ (2020) Si-Yu Li, Yu Zhang, Ya-Ning Ren, Jianpeng Liu, Xi Dai, and Lin He, “Experimental evidence for orbital magnetic moments generated by moiré-scale current loops in twisted bilayer graphene,” Phys. Rev. B 102, 121406 (2020).
  * Sharpe  _et al._ (2021) Aaron L. Sharpe, Eli J. Fox, Arthur W. Barnard, Joe Finney, Kenji Watanabe, Takashi Taniguchi, Marc A. Kastner, and David Goldhaber-Gordon, “Evidence of orbital ferromagnetism in twisted bilayer graphene aligned to hexagonal boron nitride,” Nano Letters 21, 4299–4304 (2021).
  * Tschirhart  _et al._ (2021) C. L. Tschirhart, M. Serlin, H. Polshyn, A. Shragai, Z. Xia, J. Zhu, Y. Zhang, K. Watanabe, T. Taniguchi, M. E. Huber, and A. F. Young, “Imaging orbital ferromagnetism in a moiré chern insulator,” Science 372, 1323–1327 (2021).
  * Stauber  _et al._ (2020b) T. Stauber, J. González, and G. Gómez-Santos, “Change of chirality at magic angles of twisted bilayer graphene,” Phys. Rev. B 102, 081404 (2020b).
  * Sabio  _et al._ (2008) J. Sabio, J. Nilsson, and A. H. Castro Neto, “$f$-sum rule and unconventional spectral weight transfer in graphene,” Phys. Rev. B 78, 075410 (2008).
  * Note (1) Due to our numerical procedure, there is some uncertainty in defining the cutoff-frequency and values between $D_{\rm{tot}}=-0.001t\frac{e^{2}}{\hbar^{2}}$ ($\theta=1.3^{\circ}$) and $D_{\rm{tot}}=0.004t\frac{e^{2}}{\hbar^{2}}$ ($\theta=1.0^{\circ}$) are obtained.
  * Note (2) Due to our numerical procedure, there is some uncertainty in defining the cutoff-frequency and values between $D_{\rm{chi}}=-0.0025t\frac{e^{2}}{\hbar^{2}}$ ($\theta=1.1^{\circ}$) and $D_{\rm{chi}}=0.0025t\frac{e^{2}}{\hbar^{2}}$ ($\theta=1.0^{\circ}$) are obtained.
  * Koshino (2019) Mikito Koshino, “Band structure and topological properties of twisted double bilayer graphene,” Phys. Rev. B 99, 235406 (2019).
  * Peres and Stauber (2008) N. M. R. Peres and T. Stauber, “Transport in a clean graphene sheet at finite temperature and frequency,” International Journal of Modern Physics B 22, 2529–2536 (2008).
  * G. Catalina, B. Amorim, E. V. Castro, J. M. V. P. Lopes and Peres (2019) G. Catalina, B. Amorim, E. V. Castro, J. M. V. P. Lopes and N. Peres, “Twisted bilayer graphene: Low-energy physics, electronic and optical properties,” in  _Handbook of Graphene: Volume 3_ (Wiley, Hoboken , NJ, 2019) pp. 177–231.
  * Watson and Luskin (2021) Alexander B. Watson and Mitchell Luskin, “Existence of the first magic angle for the chiral model of bilayer graphene,” J. Math. Phys. 62, 091502 (2021).
  * Becker  _et al._ (2021) Simon Becker, Mark Embree, Jens Wittsten, and Maciej Zworski, “Spectral characterization of magic angles in twisted bilayer graphene,” Phys. Rev. B 103, 165113 (2021).
  * Stauber  _et al._ (2015) T. Stauber, D. Noriega-Pérez, and J. Schliemann, “Universal absorption of two-dimensional systems,” Phys. Rev. B 91, 115407 (2015).
  * Note (3) This plateau $12\sigma_{G}$ is only obtained for twist angles which are already sufficiently close to the magic angle, i.e., the band structure for $\theta=1.05^{\circ}$ shows an avoided crossing, but does not reach this plateau, yet.
  * Bena and Simon (2011) Cristina Bena and Laurent Simon, “Dirac point metamorphosis from third-neighbor couplings in graphene and related materials,” Phys. Rev. B 83, 115404 (2011).
  * Montambaux (2012) Gilles Montambaux, “An equivalence between monolayer and bilayer honeycomb lattices,” The European Physical Journal B 85, 375 (2012).
  * Stauber (2014) Tobias Stauber, “Plasmonics in dirac systems: from graphene to topological insulators,” Journal of Physics: Condensed Matter 26, 123201 (2014).
  * Lin  _et al._ (2020) Xiao Lin, Zifei Liu, Tobias Stauber, Guillermo Gómez-Santos, Fei Gao, Hongsheng Chen, Baile Zhang, and Tony Low, “Chiral plasmons with twisted atomic bilayers,” Phys. Rev. Lett. 125, 077401 (2020).
  * He  _et al._ (2020) Wen-Yu He, David Goldhaber-Gordon, and K. T. Law, “Giant orbital magnetoelectric effect and current-induced magnetization switching in twisted bilayer graphene,” Nature Communications 11, 1650 (2020).
  * Antebi  _et al._ (2022) Ohad Antebi, Ady Stern, and Erez Berg, “In-plane orbital magnetization as a probe for symmetry breaking in strained twisted bilayer graphene,” Phys. Rev. B 105, 104423 (2022).
  * Bahamon  _et al._ (2020) Dario A. Bahamon, G. Gómez-Santos, and T. Stauber, “Emergent magnetic texture in driven twisted bilayer graphene,” Nanoscale 12, 15383–15392 (2020).
  * Khaliji  _et al._ (2020) Kaveh Khaliji, Tobias Stauber, and Tony Low, “Plasmons and screening in finite-bandwidth two-dimensional electron gas,” Phys. Rev. B 102, 125408 (2020).
  * Xie and MacDonald (2021) Ming Xie and A. H. MacDonald, “Weak-field hall resistivity and spin-valley flavor symmetry breaking in magic-angle twisted bilayer graphene,” Phys. Rev. Lett. 127, 196401 (2021).
  * Bultinck  _et al._ (2020) Nick Bultinck, Eslam Khalaf, Shang Liu, Shubhayu Chatterjee, Ashvin Vishwanath, and Michael P. Zaletel, “Ground state and hidden symmetry of magic-angle graphene at even integer filling,” Phys. Rev. X 10, 031034 (2020).
  * Zhang  _et al._ (2020) Yi Zhang, Kun Jiang, Ziqiang Wang, and Fuchun Zhang, “Correlated insulating phases of twisted bilayer graphene at commensurate filling fractions: A hartree-fock study,” Phys. Rev. B 102, 035136 (2020).
  * Lian  _et al._ (2021) Biao Lian, Zhi-Da Song, Nicolas Regnault, Dmitri K. Efetov, Ali Yazdani, and B. Andrei Bernevig, “Twisted bilayer graphene. iv. exact insulator ground states and phase diagram,” Phys. Rev. B 103, 205414 (2021).
  * Bernevig  _et al._ (2021) B. Andrei Bernevig, Biao Lian, Aditya Cowsik, Fang Xie, Nicolas Regnault, and Zhi-Da Song, “Twisted bilayer graphene. v. exact analytic many-body excitations in coulomb hamiltonians: Charge gap, goldstone modes, and absence of cooper pairing,” Phys. Rev. B 103, 205415 (2021).
  * Xie  _et al._ (2021) Fang Xie, Aditya Cowsik, Zhi-Da Song, Biao Lian, B. Andrei Bernevig, and Nicolas Regnault, “Twisted bilayer graphene. vi. an exact diagonalization study at nonzero integer filling,” Phys. Rev. B 103, 205416 (2021).
  * Rademaker  _et al._ (2019) Louk Rademaker, Dmitry A. Abanin, and Paula Mellado, “Charge smoothening and band flattening due to hartree corrections in twisted bilayer graphene,” Phys. Rev. B 100, 205114 (2019).
  * Seo  _et al._ (2019) Kangjun Seo, Valeri N. Kotov, and Bruno Uchoa, “Ferromagnetic mott state in twisted graphene bilayers at the magic angle,” Phys. Rev. Lett. 122, 246402 (2019).
  * González and Stauber (2021) J. González and T. Stauber, “Magnetic phases from competing hubbard and extended coulomb interactions in twisted bilayer graphene,” Phys. Rev. B 104, 115110 (2021).
  * 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).
  * Pedersen  _et al._ (2008) Thomas G. Pedersen, Christian Flindt, Jesper Pedersen, Antti-Pekka Jauho, Niels Asger Mortensen, and Kjeld Pedersen, “Optical properties of graphene antidot lattices,” Phys. Rev. B 77, 245431 (2008).