# Magic Angle Butterfly in Twisted Trilayer Graphene

Fedor K. Popov and Grigory Tarnopolsky  Department of Physics, New York
University, New York, NY 10003, USA  Department of Physics, Carnegie Mellon
University, Pittsburgh, PA 15213, USA

###### Abstract

We consider a configuration of three stacked graphene monolayers with
commensurate twist angles $\theta_{12}/\theta_{23}=p/q$, where $p$ and $q$ are
coprime integers with $0<p<|q|$ and $q$ can be positive or negative. We study
this system using the continuum model in the chiral limit when interlayer
coupling terms between $\textrm{AA}_{12}$ and $\textrm{AA}_{23}$ sites of the
moiré patterns $12$ and $23$ are neglected. There are only three inequivalent
displacements between the moiré patterns $12$ and $23$, at which the three
monolayers’ Dirac zero modes are protected. Remarkably, for these
displacements and an arbitrary $p/q$ we discover exactly flat bands at an
infinite set of twist angles (magic angles). We provide theoretical
explanation and classification of all possible configurations and topologies
of the flat bands.

##  I Introduction and Summary

Two graphene layers placed one on top of the other with a relative small twist
angle form a periodic moiré pattern, which alters significantly the low energy
electronic spectrum. At the special twist ”magic” angle $\theta_{*}\approx
1.05^{\circ}$ a dramatic flattening of the lowest energy bands was observed in
[1, 2]. These almost dispersionless electronic bands hinted to a possibility
of interesting strongly interacting phenomena, which was subsequently
confirmed through a series of experimental studies [3, 4, 5]. These
tantalizing experiments have inspired numerous theoretical and experimental
investigations [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34], and the field continues
to advance with fascinating new proposals and insights.

Multiple layers of graphene stacked on top of each other with small relative
twist angles [35, 36, 37, 38, 36, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48]
provide even greater versatility compared to Twisted Bilayer Graphene (TBG)
due to the increased number of parameters. Initial theoretical investigations
of alternating-twist trilayer graphene (aTTG) [35] demonstrated a similar
flattening of electronic bands at ”magic” angles, ultimately leading to the
experimental discovery of various correlated phenomena [49, 50, 51, 52, 53,
54] at the angles predicted in [35]. The study of the interaction effects in
such systems remains an active area of ongoing theoretical investigations [55,
56, 57].

In [36] a general configuration of twisted trilayer graphene was proposed,
where the three layers are consecutively twisted by small angles
$\theta_{12}=p\theta$ and $\theta_{23}=q\theta$, where $p$ and $q$ are coprime
integers $0<p<q$, as shown in Fig. 1.

Figure 1: Picture of a pq-Twist Trilayer Graphene (pqTTG) with $(p,q)=(1,2)$.

The analysis of electronic bands in this system, was initially performed using
the continuum model with equal interlayer coupling terms between AA and AB
sites of the moiré patterns. However, this model did not exhibit a distinct
phenomenon of band flattening.

The existence of magic angles and perfectly flat bands was recently discovered
in the case of equal-twist trilayer graphene (eTTG) [58], where the twist
angles between layers 1 and 2 ($\theta_{12}$) and between layers 2 and 3
($\theta_{23}$) are equal. The chiral limit of the twisted graphene continuum
model [29, 30] was crucial in revealing these magic angles and flat bands. In
this limit, the interlayer coupling terms between AA sites in the moiré
pattern are disregarded, resulting in the Hamiltonian exhibiting exactly flat
bands at an infinite series of magic angles.

In addition to the twist angles, the spectrum of a general twisted trilayer
graphene (TTG) system is influenced by the relative displacement d between the
two moiré patterns formed by layers $12$ and $23$. The investigation of equal-
twist trilayer graphene (eTTG) without any displacement between the moiré
patterns ($\textbf{d}=0$, AAA stacking) revealed an intriguing connection
between the flat bands of twisted bilayer graphene and eTTG, establishing a
relation between the magic angles in these distinct systems. A subsequent
study by Guerci et al. [59] (see also [60, 61]) demonstrated that eTTG system
exhibits a ”moiré of moiré” pattern, resulting in large triangular regions of
ABA and BAB stacking
($\textbf{d}=\pm\frac{1}{3}(\textbf{a}_{1}-\textbf{a}_{2})$, where
$\textbf{a}_{1,2}$ are the single moiré lattice unit vectors) separated by
smaller AAA regions. By employing the chiral model with ABA or BAB stacking of
the graphene layers, another set of magic angles and perfectly flat bands were
unveiled. Furthermore, it was shown that the flat bands in eTTG can possess
the Chern number $C=2$, and the mathematical origin of such flat bands was
explained.

In this letter, we analyze the general twisted trilayer graphene configuration
in which the consecutive layers are twisted at small but commensurate angles
$\theta_{12}=p\theta$ and $\theta_{23}=q\theta$, where $p$ and $q$ are coprime
integer numbers $0<p<|q|$ and $q$ can be positive or negative. We refer to
this system as $pq$-Twist Trilayer Graphene (pqTTG), and it is schematically
depicted in Fig. 1. A configuration of this type has been recently achieved in
an experiment [54]. We analyze the continuum model for this trilayer graphene
system and take the chiral limit, where interlayer coupling terms are
$w_{\textrm{AB}}\approx 110\textrm{meV}$ and $w_{\textrm{AA}}=0$.

Figure 2: The Magic Angle Butterfly. Plot of the leading magic angles $\alpha_{12}$ for all possible ratios $p/q$ of coprime integers $p$ and $q$ with $0<p<|q|$ and up to $|q|=20$ and displacements $\textbf{d}=0,\pm\textbf{d}_{0}$. Near the center ($p/q=0$) the magic angles approach those of TBG. $(p,q)_{\textbf{d}}$ | $\alpha_{12}^{(1)}$ | $\alpha_{12}^{(2)}$ | $\alpha_{12}^{(3)}$ | $\alpha_{12}^{(4)}$
---|---|---|---|---
$(1,4)_{0}$ | $0.5673$ | $1.4042$ | $1.8253$ | $1.8295$
$(1,4)_{\pm\textbf{d}_{0}}$ | $0.5672$ | $1.6416_{4}$ | $1.7963$ | $2.0260_{4}$
$(1,3)_{0}$ | $0.5542$ | $1.2025$ | $2.8487$ | $4.0074$
$(1,3)_{\pm\textbf{d}_{0}}$ | $0.5527$ | $1.2464$ | $1.3919$ | $2.0304$
$(1,2)_{0}$ | $0.5194$ | $0.9589$ | $1.0255_{4}$ | $1.5309_{4}$
$(1,2)_{\pm\textbf{d}_{0}}$ | $0.5094$ | $1.0698$ | $1.6226$ | $2.9886$
$(2,3)_{0}$ | $0.4403$ | $0.4549$ | $0.4954_{4}$ | $0.7662_{4}$
$(2,3)_{\pm\textbf{d}_{0}}$ | $0.4527$ | $0.4548$ | $0.4668$ | $0.5010$
$(3,4)_{0}$ | $0.4356$ | $0.4710$ | $0.4941_{4}$ | $0.5012$
$(3,4)_{\pm\textbf{d}_{0}}$ | $0.4269$ | $0.4302$ | $0.4612$ | $0.4646$
$(1,1)_{0}$ | $0.8283_{4}$ | $3.1412_{4}$ | $5.3053_{4}$ | $7.4621_{4}$
$(1,1)_{\pm\textbf{d}_{0}}$ | $0.3771$ | $1.1967_{4}$ | $1.7549$ | $2.4136_{4}$
Table 1: The table of magic angles for various $p,q>0$

For small twist angles we assume that the moiré patterns $12$ and $23$ are
formed by vectors $p\textbf{a}_{1,2}$ and $q\textbf{a}_{1,2}$, where
$\mathbf{a}_{1,2}=(4\pi/3k_{\theta})(\pm\sqrt{3}/2,1/2)$ with
$k_{\theta}=2k_{D}\sin(\theta/2)\approx k_{D}\theta$ and
$k_{D}=4\pi/3\sqrt{3}a$ is the Dirac momentum of the monolayer graphene with
lattice constant $a\approx 1.42$Å. Similarly to the previously discussed eTTG
configuration, pqTTG also has a moiré of moiré pattern, resulting in local
variations of the displacement vector d between the moiré patterns $12$ and
$23$. There are only three inequivalent displacement vectors, at which the
three monolayers’ Dirac zero modes are protected. These are $\textbf{d}=0$ and
$\textbf{d}=\pm\textbf{d}_{0}$, where
$\textbf{d}_{0}=\frac{1}{3p|q|}(\textbf{a}_{1}-\textbf{a}_{2})$ and the
spectrum is invariant under the shifts
$\textbf{d}\to\textbf{d}+\frac{1}{pq}\textbf{a}_{1,2}$. The values of the
magic angles are identical for the displacements $\pm\textbf{d}_{0}$.
Remarkably, in the chiral limit and these displacements the electronic energy
spectrum exhibits perfectly flat bands at an infinite sequence of magic angles
for any combination of coprime integers $p$ and $q$. We introduce
dimensionless twist parameters $\alpha_{12}=\alpha/p$ and
$\alpha_{23}=\alpha/q$ where $\alpha=w_{\textrm{AB}}/(v_{F}k_{D}\theta)$, and
$v_{F}\approx 10^{6}$m/s is the monolayer graphene Fermi velocity (for brevity
we also refer to $\alpha_{12}$, $\alpha_{23}$ and $\alpha$ as twist angles).
We plot magic angles $\alpha_{12}$ as a function of $p/q$ in Fig. 2. Plot of
the magic angles $(\alpha_{12}$, $\alpha_{23})$ for different $p/q$ is shown
in Fig. 3. Finally we present Tables 1 and 2 of the first four magic angles
$\alpha_{12}$ for various values of $p$ and $q$ and displacements
$\textbf{d}=0$ and $\textbf{d}=\pm\textbf{d}_{0}$.

Figure 3: Plot of the leading magic angles $(\alpha_{12},\alpha_{23})$ for all possible ratios $p/q$ of coprime integers $p$ and $q$ with $0<p<|q|$ and up to $|q|=20$. The black dashed line is a semicirlce $\alpha_{12}^{2}+\alpha_{23}^{2}=1/3$ derived in [38] using perturbation theory. $(p,q)_{\textbf{d}}$ | $\alpha_{12}^{(1)}$ | $\alpha_{12}^{(2)}$ | $\alpha_{12}^{(3)}$ | $\alpha_{12}^{(4)}$
---|---|---|---|---
$(1,-1)_{0}$ | $0.4141$ | $1.5706$ | $2.6526$ | $3.7310$
$(1,-1)_{\pm\textbf{d}_{0}}$ | $0.5072$ | $2.2295$ | $2.9043$ | $4.4080$
$(3,-4)_{0}$ | $0.4821$ | $0.4887_{4}$ | $0.6493$ | $0.6626_{4}$
$(3,-4)_{\pm\textbf{d}_{0}}$ | $0.5121$ | $0.6167$ | $0.6520$ | $0.6954$
$(2,-3)_{0}$ | $0.4912_{4}$ | $0.5054$ | $1.2890$ | $1.3259_{4}$
$(2,-3)_{\pm\textbf{d}_{0}}$ | $0.9458$ | $1.0945$ | $1.7939$ | $2.2357$
$(1,-2)_{0}$ | $0.5131$ | $1.5014$ | $2.0490$ | $3.2157$
$(1,-2)_{\pm\textbf{d}_{0}}$ | $0.5857$ | $0.7614$ | $1.3797$ | $1.8178_{4}$
$(1,-3)_{0}$ | $0.5574$ | $1.4373$ | $2.6738_{4}$ | $3.9906_{4}$
$(1,-3)_{\pm\textbf{d}_{0}}$ | $0.5587$ | $1.2915$ | $1.3920$ | $2.1730$
$(1,-4)_{0}$ | $0.5685$ | $1.7225_{4}$ | $1.8548_{4}$ | $3.2869$
$(1,-4)_{\pm\textbf{d}_{0}}$ | $0.5685$ | $1.7130$ | $1.7298$ | $1.9467$
Table 2: The table of magic angles for various $p,-q>0$

Below we formulate the continuum model for twisted trilayer graphene, present
the criteria for the emergence of perfectly flat bands and provide a complete
classification of the structures of the flat bands at magic angles.

##  II Continuum model for Twisted Trilayer Graphene

We consider a system of three stacked graphene monolayers, where each layer
$\ell=1,2,3$ is rotated counterclockwise by an angle $\theta_{\ell}$ around an
atom site and then shifted by a vector $\textbf{d}_{\ell}$, so atoms in each
layer are parametrized by
$\textbf{r}=R_{\theta_{\ell}}(\textbf{R}+\tau_{\alpha})+\textbf{d}_{\ell}$,
where $R_{\theta}=e^{-i\theta\sigma_{y}}$ is the rotation matrix and R and
$\tau_{\alpha}$ run over the lattice and sub-lattice sites. The continuum
model Hamiltonian for twisted trilayer graphene can be written as [35]:

$\displaystyle
H=\left(\begin{array}[]{ccc}-iv_{F}\bm{\sigma}_{\theta_{1}}\bm{\nabla}&T^{12}(\textbf{r}-\textbf{d}_{12})&0\\\
T^{12{\dagger}}(\textbf{r}-\textbf{d}_{12})&-iv_{F}\bm{\sigma}_{\theta_{2}}\bm{\nabla}&T^{23}(\textbf{r}-\textbf{d}_{23})\\\
0&T^{23{\dagger}}(\textbf{r}-\textbf{d}_{23})&-iv_{F}\bm{\sigma}_{\theta_{3}}\bm{\nabla}\\\
\end{array}\right)\,,$

where $\bm{\sigma}_{\theta}\equiv
e^{i\frac{\theta}{2}\sigma_{z}}\bm{\sigma}e^{-i\frac{\theta}{2}\sigma_{z}}$,
$\bm{\sigma}=(\sigma_{x},\sigma_{y})$ and
$\textbf{d}_{\ell\ell^{\prime}}=\frac{1}{2}(\textbf{d}_{\ell}+\textbf{d}_{\ell^{\prime}}+i\cot(\theta_{\ell^{\prime}\ell}/2)\sigma_{y}(\textbf{d}_{\ell}-\textbf{d}_{\ell^{\prime}}))$
is the moiré pattern displacement vector. The moiré potential between adjacent
layers $\ell$ and $\ell^{\prime}$ is

$\displaystyle
T^{\ell\ell^{\prime}}(\textbf{r})=\sum_{n=1}^{3}T^{\ell\ell^{\prime}}_{n}e^{-i\textbf{q}_{n}^{\ell\ell^{\prime}}\textbf{r}}\,,$
(1)

where
$T^{\ell\ell^{\prime}}_{n+1}=w_{\textrm{AA}}^{\ell\ell^{\prime}}\sigma_{0}+w_{\textrm{AB}}^{\ell\ell^{\prime}}(\sigma_{x}\cos
n\phi+\sigma_{y}\sin n\phi)$ and

$\displaystyle\textbf{q}_{1}^{\ell\ell^{\prime}}=2k_{D}\sin(\theta_{\ell^{\prime}\ell}/2)R_{\phi_{\ell\ell^{\prime}}}(0,-1),\;\;\textbf{q}^{\ell\ell^{\prime}}_{2,3}=R_{\pm\phi}\textbf{q}^{\ell\ell^{\prime}}_{1}$
(2)

with $\theta_{\ell\ell^{\prime}}=\theta_{\ell}-\theta_{\ell^{\prime}}$,
$\phi_{\ell\ell^{\prime}}=(\theta_{\ell}+\theta_{\ell^{\prime}})/2$,
$\phi=2\pi/3$. The coupling between adjacent layers $\ell$ and $\ell^{\prime}$
is characterized by two parameters $w_{\textrm{AA}}^{\ell\ell^{\prime}}$ and
$w_{\textrm{AB}}^{\ell\ell^{\prime}}$ representing intra- and intersub-lattice
couplings. The chiral limit corresponds to
$w_{\textrm{AA}}^{\ell\ell^{\prime}}=0$. Also using translation invariance we
can make a replacement $\textbf{r}\to\textbf{r}+\mathbf{d}_{12}$, and
therefore the Hamiltonian depends only on a single displacement vector
$\mathbf{d}=\mathbf{d}_{23}-\mathbf{d}_{12}$.

We consider a trilayer configuration where the relative twists are
commensurate $\theta_{21}=p\theta$ and $\theta_{32}=q\theta$ and $p$ and $q$
are two coprime integers, which satisfy $0<p<|q|$. For a small angle $\theta$
we can set $\phi_{ll^{\prime}}=0$ leading to
$\textbf{q}_{1}^{12}=p\textbf{q}_{1}$ and
$\textbf{q}_{1}^{23}=q\textbf{q}_{1}$ and $\textbf{q}_{1}=k_{\theta}(0,-1)$
with $k_{\theta}=2k_{D}\sin(\theta/2)\approx k_{D}\theta$. Thus we obtain the
following Hamiltonian

$\displaystyle
H_{\textrm{pqTTG}}=\left(\begin{array}[]{ccc}-iv_{F}\bm{\sigma}_{-p\theta}\bm{\nabla}&T(p\textbf{r})&0\\\
T^{{\dagger}}(p\textbf{r})&-iv_{F}\bm{\sigma}\bm{\nabla}&T(q(\textbf{r}-\textbf{d}))\\\
0&T^{{\dagger}}(q(\textbf{r}-\textbf{d}))&-iv_{F}\bm{\sigma}_{q\theta}\bm{\nabla}\\\
\end{array}\right)\,,$ (6)

where we assume equal couplings between layers. For a small twist angle
$\theta$ we can neglect the phase factors in the Pauli matrices
$\bm{\sigma}_{-p\theta}\to\bm{\sigma}$ and
$\bm{\sigma}_{q\theta}\to\bm{\sigma}$.

The moiré Brillouin zone (mBZ) for this Hamiltonian is depicted in Fig. 4. The
reciprocal moiré unit vectors are
$\mathbf{b}_{1,2}=\mathbf{q}_{2,3}-\mathbf{q}_{1}$. In the coordinate space
the unit vectors are $\mathbf{a}_{1,2}=(4\pi/3k_{\theta})(\pm\sqrt{3}/2,1/2)$.
It is useful to introduce complex coordinates $z,\bar{z}=\mathbf{r}_{x}\pm
i\mathbf{r}_{y}$ in real space and $k,\bar{k}=\mathbf{k}_{1}\pm
i\mathbf{k}_{2}$ in momentum space.

Figure 4: Original Brillouin zones of three graphene layers with their Dirac
points $K_{1}$, $K_{2}$ and $K_{3}$ and the moiré Brillouin zones for pqTTG.
The layers are twisted by the angles $p\theta$ and $q\theta$, where $p$ and
$q$ are coprime integers such that $0<p<|q|$. The wave vector k is zero at the
Dirac point $K_{2}$. We neglect a relative rotation between vectors
$\textbf{q}_{1}^{12}$ and $\textbf{q}_{1}^{23}$ and assume that
$\textbf{q}_{1}^{12}=p\textbf{q}_{1}$ and
$\textbf{q}_{1}^{23}=q\textbf{q}_{1}$. (a) Case $(p,q)=(1,-3)$. (b) Case
$(p,q)=(2,3)$.

The Hamiltonian (6) acts on a spinor
$\Psi(\textbf{r})=(\psi_{1},\chi_{1},\psi_{2},\chi_{2},\psi_{3},\chi_{3})$,
where the indices $1,2,3$ represent the graphene layer. Introducing the
dimensionless twist parameter $\alpha=w_{\textrm{AB}}/(v_{F}k_{\theta})$, and
writing the Hamiltonian (6) in the sublattice basis
$\Psi(\textbf{r})=(\psi_{1},\psi_{2},\psi_{3},\chi_{1},\chi_{2},\chi_{3})$ we
obtain

$\displaystyle\mathcal{H}_{\textrm{pqTTG}}=\begin{pmatrix}\mathcal{M}(\textbf{r})&\mathcal{D}^{{\dagger}}(\textbf{r})\\\
\mathcal{D}(\textbf{r})&\mathcal{M}(\textbf{r})\end{pmatrix}\,,$ (7)

where we have rescaled coordinates $\textbf{r}\to k_{\theta}\textbf{r}$ and
the Hamiltonian, so the energies of (7) are measured in units of
$v_{F}k_{\theta}$. The operators $\mathcal{D}$ and $\mathcal{M}$ are

$\displaystyle\mathcal{D}(\textbf{r})=\left(\begin{array}[]{ccc}-2i\bar{\partial}&\alpha
U_{1}(p\textbf{r})&0\\\ \alpha U_{1}(-p\textbf{r})&-2i\bar{\partial}&\alpha
U_{1}(q(\textbf{r}-\textbf{d}))\\\ 0&\alpha
U_{1}(q(\textbf{d}-\textbf{r}))&-2i\bar{\partial}\\\ \end{array}\right)\,,$
(11)
$\displaystyle\mathcal{M}(\mathbf{r})=\frac{w_{\textrm{AA}}}{w_{\textrm{AB}}}\left(\begin{array}[]{ccc}0&U_{0}(p\textbf{r})&0\\\
U_{0}(-p\textbf{r})&0&U_{0}(q(\textbf{r}-\textbf{d}))\\\
0&U_{0}(q(\textbf{d}-\textbf{r}))&0\\\ \end{array}\right)\,,$ (15)

where $\partial,\bar{\partial}=\frac{1}{2}(\partial_{x}\mp i\partial_{y})$ are
complex derivatives and we defined
$U_{m}(\mathbf{r})=\sum^{3}_{n=1}\omega^{m(n-1)}e^{-i\textbf{q}_{n}\textbf{r}}$
with $\omega=e^{i\phi}$. The Bloch states
$\Psi_{\textbf{k}}(\textbf{r})=(\psi_{\textbf{k}}(\textbf{r}),\chi_{\textbf{k}}(\textbf{r}))$
of (7) are parametrized by the wave vector k from mBZ and satisfy the
following boundary conditions

$\displaystyle\Psi_{\mathbf{k}}(\mathbf{r}+\mathbf{a}_{1,2})=e^{i\mathbf{k}\,\mathbf{a}_{1,2}}U_{\phi}\Psi_{\mathbf{k}}(\mathbf{r})\,,$
(16)

where the matrix
$U_{\phi}=\mathds{1}_{\textrm{AB}}\otimes\operatorname{diag}\left(\omega^{p},1,\omega^{-q}\right)$.

Below we focus on the chiral limit $w_{\textrm{AA}}=0$. In this case the
Hamiltonian (7) is particle-hole symmetric
$\\{\mathcal{H}_{\textrm{pqTTG}},\sigma_{z}\otimes\mathds{1}\\}=0$. The chiral
and anti-chiral zero energy equations are

$\displaystyle\mathcal{D}(\textbf{r})\psi_{\textbf{k}}(\textbf{r})=0,\quad\mathcal{D}^{{\dagger}}(\textbf{r})\chi_{\textbf{k}}(\textbf{r})=0\,,$
(17)

and they have solutions for every wave vector k of the mBZ (flat bands) only
at the special values (magic angles) of the twist angle $\alpha$. We can
represent the chiral operator $\mathcal{D}(\mathbf{r})$ in (15) in the form

$\displaystyle\mathcal{D}(\mathbf{r})=-2i\bar{\partial}+\alpha\bar{A}(\textbf{r})\,,$
(18)

where we introduced the matrix vector-potential $\bar{A}(\textbf{r})$. The
magic angles $\alpha$ can be found numerically as a generalized eigenvalues of
the operators $-2i\bar{\partial}$ and $\bar{A}(\textbf{r})$ [62, 63].
Generically the magic angles $\alpha$ form an infinite set of isolated points
in the entire complex plane. In order to have magic angles, it is necessary
for the chiral operator $\mathcal{D}(\textbf{r})$ to adhere to specific
symmetries. Namely we would like the chiral operator to have three Dirac zero
modes $\psi_{K_{1}}$, $\psi_{K_{2}}$ and $\psi_{K_{3}}$ at the Dirac points
$K_{1}=p\textbf{q}_{1}$, $K_{2}=0$ and $K_{3}=-q\textbf{q}_{1}$ for an
arbitrary twist angle $\alpha$. This is guaranteed by $C_{3z}$ symmetry
together with the particle-hole symmetry [29, 36]. A non-zero displacement d
can break the $C_{3z}$ symmetry of the Hamiltonian. One can show that spectrum
of the Hamiltonian (7) is invariant under the following shifts of the
displacement

$\displaystyle\mathbf{d}\to\mathbf{d}+\frac{1}{pq}\textbf{a}_{1,2}\,$ (19)

and the inequivalent displacements that preserve the three Dirac zero modes
are $\textbf{d}=0$ and $\textbf{d}=\pm\textbf{d}_{0}$, where

$\displaystyle\textbf{d}_{0}=\frac{1}{3p|q|}\left(\textbf{a}_{1}-\textbf{a}_{2}\right)\,.$
(20)

The magic angles of the the Hamiltonian are identical for the displacements
$\textbf{d}=\pm\textbf{d}_{0}$, because they are related to each other by the
$C_{2z}\mathcal{T}$ symmetry, that exchanges the chiral and anti-chiral
components of the Bloch states.

##  III Origin of the flat bands

In this section we derive general properties of the zero-modes of the
Hamiltonian (7) in the chiral limit $w_{\textrm{AA}}=0$. As we discussed
above, the chiral operator admits three Dirac zero modes $\psi_{K_{1}}$,
$\psi_{K_{2}}$ and $\psi_{K_{3}}$ at the Dirac points $K_{1}=p\textbf{q}_{1}$,
$K_{2}=0$ and $K_{3}=-q\textbf{q}_{1}$

$\displaystyle\mathcal{D}(\mathbf{r})\psi_{K_{i}}(\mathbf{r})=0,\quad i=1,2,3$
(21)

for the displacements $\mathbf{d}=0,\pm\mathbf{d}_{0}$ and an arbitrary twist
angle $\alpha$. Using that $\textrm{Tr}\bar{A}(\textbf{r})=0$ and equations
(21) one can derive that the Wronskian of the Dirac spinors

$\displaystyle
W(\mathbf{r})\equiv\det\left(\psi_{K_{1}},\psi_{K_{2}},\psi_{K_{3}}\right)$
(22)

satisfies $\bar{\partial}W(\mathbf{r})=0$ [64],[59]. (Note that the Wronskian
in this case is simply a triple product:
$W(\textbf{r})=\psi_{K_{1}}\cdot(\psi_{K_{2}}\times\psi_{K_{3}})$). Therefore
we must conclude that the Wronskian is a constant $W(\textbf{r})=W$, since we
can not have a non-constant holomorphic function on a compact manifold, which
is a torus in our case.

First we prove that if $W\neq 0$ there is no flat band. $W\neq 0$ implies that
$\psi_{K_{i}}(\textbf{r})$, $i=1,2,3$ as $3$-dimensional vectors are linearly
independent at each point $\mathbf{r}$ of the moiré unit cell, and form a
basis. Therefore any other chiral zero mode $\psi_{\mathbf{k}}(\mathbf{r})$
($\mathcal{D}\psi_{\textbf{k}}=0$) can be expanded as

$\displaystyle\psi_{\mathbf{k}}(\mathbf{r})=\mathcal{C}_{1}(\mathbf{r})\psi_{K_{1}}(\mathbf{r})+\mathcal{C}_{2}(\mathbf{r})\psi_{K_{2}}(\mathbf{r})+\mathcal{C}_{3}(\mathbf{r})\psi_{K_{3}}(\mathbf{r})\,,$

where the scalar functions $\mathcal{C}_{1,2,3}$ must depend on the wave
vector k in order for $\psi_{\mathbf{k}}$ to satisfy the the Bloch boundary
conditions (16). Applying the operator $\mathcal{D}(\mathbf{r})$ to both parts
of the above equality we conclude that
$\bar{\partial}\mathcal{C}_{1,2,3}(\mathbf{r})=0$. And since
$\psi_{\mathbf{k}}$ and $\psi_{K_{i}}$ are finite everywhere
$\mathcal{C}_{i}(\mathbf{r})$ can only be constants. But then
$\psi_{\textbf{k}}$ would violate the boundary conditions (16). Therefore if
$W\neq 0$ we can not have other zero energy solutions (and thus flat band)
except the Dirac zero modes.

Now, we prove that if $W=0$, we necessarily have a flat band. $W=0$ implies
that the Dirac spinors $\psi_{K_{i}}(\mathbf{r})$ as vectors are linearly
dependent at every point $\mathbf{r}$ of the moiré unit cell. We stress that
it does not mean that they are linearly dependent functions in the Hilbert
space. There are two possible scenarios.

Scenario 1. All Bloch spinors of the zero energy flat bands are collinear at
every point $\mathbf{r}$ of the moiré unit cell. We refer to such case as rank
1 flat bands. In this situation any flat band’s wave spinor can be described
by the following equation:

$\displaystyle\psi_{\mathbf{k}}(\mathbf{r})=\mathcal{C}(\mathbf{r})\psi_{K_{1}}(\mathbf{r})\,,$
(23)

where $\mathcal{C}(\mathbf{r})$ is some scalar function. If we apply the
operator $\mathcal{D}(\textbf{r})$ to both sides of this equation we obtain
$\bar{\partial}\mathcal{C}(\mathbf{r})=0$ and thus $\mathcal{C}$ is a
constant, which violates the boundary conditions (16) for the spinor
$\psi_{\textbf{k}}$. Therefore the eq. (23) is impossible unless
$\psi_{K_{1}}(\mathbf{r})$ has zeros. If $\psi_{K_{1}}(\mathbf{r})$ has $n$
zeros at the points $\textbf{r}_{\lambda}$, $\lambda=1,\dots,n$, we can
construct $n$ flat bands, using the following wave functions:

$\displaystyle\psi^{(\lambda)}_{\mathbf{k}}(\mathbf{r})=f_{\textbf{k}-K_{1}}(z-z_{\lambda})\psi_{K_{1}}(\mathbf{r})\,,$
(24)

where
$z_{\lambda}=\left(\mathbf{r}_{\lambda}\right)_{x}+i\left(\mathbf{r}_{\lambda}\right)_{y}$
and we introduced the following meromorphic function

$\displaystyle
f_{\textbf{k}}(z)=e^{i\frac{\textbf{k}\textbf{a}_{1}}{a_{1}}z}\frac{\vartheta_{1}(z/a_{1}-k/b_{2}|\omega)}{\vartheta_{1}(-k/b_{2}|\omega)\vartheta_{1}(z/a_{1}|\omega)}\,.$
(25)

Normalization of the function $f_{\textbf{k}}(z)$ is chosen such that
$f_{\textbf{k}+\textbf{b}_{1,2}}(z)=f_{\textbf{k}}(z)$ and one can compute
that the Chern number of the flat bands in (24) is $C=1$.

Scenario 2. The other case is that all wave spinors of the flat bands form a
two-dimensional vector space at every point $\mathbf{r}$ of the moiré unit
cell. We refer to such a case as rank 2 flat bands. In this case we have

$\displaystyle\psi_{K_{3}}(\mathbf{r})=\mathcal{C}_{1}(\mathbf{r})\psi_{K_{1}}(\mathbf{r})+\mathcal{C}_{2}(\mathbf{r})\psi_{K_{2}}(\mathbf{r})\,,$
(26)

where $\mathcal{C}_{1}(\mathbf{r})$ and $\mathcal{C}_{2}(\mathbf{r})$ are some
scalar functions. If we assume that the spinors $\psi_{K_{1}}(\textbf{r})$ and
$\psi_{K_{2}}(\textbf{r})$ are linearly independent at every point r of the
moiré unit cell and apply the operator $\mathcal{D}(\textbf{r})$ to both sides
of the equation (26) we obtain $\bar{\partial}\mathcal{C}_{1,2}(\mathbf{r})=0$
and thus $\mathcal{C}_{1,2}$ are constants, which violates the boundary
conditions (16) for the spinor $\psi_{K_{3}}$. Therefore we conclude that
there exists at least one point $\textbf{r}_{0}$ where the vectors
$\psi_{K_{1}}(\mathbf{r}_{0})$ and $\psi_{K_{2}}(\mathbf{r}_{0})$ are linearly
dependent. More generally, there could be $n$ such points
$\mathbf{r}_{\lambda}$, $\lambda=1,\dots,n$, where

$\displaystyle
c_{1\lambda}\psi_{K_{1}}(\mathbf{r}_{\lambda})+c_{2\lambda}\psi_{K_{2}}(\mathbf{r}_{\lambda})=0\,,$
(27)

for $\lambda=1,\dots,n$. (Notice that the coefficients $c_{1,2\lambda}$ are
not equal to $\mathcal{C}_{1,2}(\textbf{r}_{\lambda})$). Now it is easy to see
that for each zero point $\mathbf{r}_{\lambda}$ we can construct a chiral zero
mode at every point k of the mBZ:

$\displaystyle\psi^{(\lambda)}_{\mathbf{k}}(\mathbf{r})=\big{(}c_{1\lambda}$
$\displaystyle f_{\mathbf{k}-K_{1}}(z-z_{\lambda})\psi_{K_{1}}(\mathbf{r})$
(28)
$\displaystyle\qquad\quad+c_{2\lambda}f_{\mathbf{k}-K_{2}}(z-z_{\lambda})\psi_{K_{2}}(\mathbf{r})\big{)}\,.$

The function $f_{\textbf{k}}$ defined in (25) behaves as

$\displaystyle
f_{\textbf{k}}(z-z_{\lambda})\to\frac{a_{1}}{\vartheta_{1}^{\prime}(0|\omega)}\frac{1}{z-z_{\lambda}},\quad
z\to z_{\lambda}$ (29)

and the poles of the functions $f_{\mathbf{k}-K_{1}}(z)$ and
$f_{\mathbf{k}-K_{2}}(z)$ are cancelled at the point $\mathbf{r}_{\lambda}$,
thus $\psi_{\textbf{k}}^{(\lambda)}$ is finite at every point of the moiré
unit cell. A similar construction within the context of eTTG was proposed by
Guerci et al. [59].

Since the function $f_{\textbf{k}}$ in (25) is periodic in the mBZ, but has a
pole at $\mathbf{k}=0$ the function $\psi^{(\lambda)}_{\mathbf{k}}$ can have
two poles in mBZ if both coefficients $c_{1\lambda}$ and $c_{1\lambda}$ are
not zero. In this case the flat band has the Chern number $C=2$. If either of
the coefficients $c_{1\lambda}$ or $c_{2\lambda}$ is zero, the function
$\psi^{(\lambda)}_{\mathbf{k}}$ has one simple pole within the mBZ, we have a
rank $1$ flat sub-band with the Chern number $C=1$.

Thus, the rank 2 flat band is characterized by two wave functions
$\psi_{K_{1}}(\mathbf{r})$ and $\psi_{K_{2}}(\mathbf{r})$. These wave spinors
are linearly independent vectors at every point in the moiré unit cell besides
$n$ isolated points $\textbf{r}_{\lambda}$. Among these points, there may
exist $n_{1}$ points where either of the spinors $\psi_{K_{1}}$ or
$\psi_{K_{2}}$ is zero. These points give rise to $n_{1}$ flat bands with the
Chern number $C=1$. At the rest of $n_{2}=n-n_{1}$ points $\psi_{K_{1}}$ and
$\psi_{K_{2}}$ are non-zero and linearly dependent. These points give rise to
$n_{2}$ flat bands with the Chern number $C=2$.

We note that the above construction can be easliy generalized for the case of
twisted $\ell$-layer graphene in the chiral limit and shows that one can have
flat bands with the Chern number up to $C=\ell-1$.

Finally we note that we can have only a single Dirac cone on top of the flat
bands. To demonstrate this, we observe that if there are two Dirac cones on
top of the flat bands, their zero energy wave spinors must be linearly
dependent at some point. Consequently, these wave functions would generate a
rank 2 flat band, contradicting the assumption that they constitute distinct
zero-energy states.

Below we elaborate more on the possible structures of the flat bands in pqTTG.
For this we separately consider cases of $\textbf{d}=0$ and
$\textbf{d}=\pm\mathbf{d}_{0}$ displacements. We summarize all possible
configurations of the flat bands in the Fig. 5 and 6.

##  IV Structures of the flat bands

In this section we discuss configurations of the flat bands and their Chern
numbers for two different cases of the displacement: $\textbf{d}=0$ and
$\textbf{d}=\pm\mathbf{d}_{0}$.

###  IV.1 Displacement $\textbf{d}=0$

Absence of the $C=2$ flat bands. In the case of zero displacement one can
construct anti-chiral zero modes $\chi_{\textbf{k}}$
($\mathcal{D}^{{\dagger}}\chi_{\textbf{k}}=0$) from chiral zero modes
$\psi_{\mathbf{k}}$ ($\mathcal{D}\psi_{\textbf{k}}=0$) as

$\displaystyle\chi_{\mathbf{k}_{1,2}}(\mathbf{r})=Q\bar{\psi}_{\mathbf{k}_{1,2}}(-\mathbf{r})\,$
$\displaystyle\chi_{\mathbf{k}}(\mathbf{r})=Q\left(\bar{\psi}_{\mathbf{k}_{1}}(\mathbf{r})\times\bar{\psi}_{\mathbf{k}_{2}}(\mathbf{r})\right),$
(30)

where $Q=\textrm{diag}(1,-1,1)$,
$\mathbf{k}=\mathbf{k}_{*}-\mathbf{k}_{1}-\mathbf{k}_{2}$ and
$\mathbf{k}_{*}=(p-q)\mathbf{q}_{1}$. By combining the above formulas we can
also construct an additional chiral zero mode $\psi_{\mathbf{k}}$ at the point
$\mathbf{k}=\mathbf{k}_{*}-\mathbf{k}_{1}-\mathbf{k}_{2}$:

$\displaystyle\psi_{\mathbf{k}}(\mathbf{r})=\psi_{\mathbf{k}_{1}}(-\mathbf{r})\times\psi_{\mathbf{k}_{2}}(-\mathbf{r})\,.$
(31)

Moreover one can check that that the following function

$\displaystyle
v(\mathbf{r})=\psi_{\mathbf{k}_{1}}(\mathbf{r})\cdot\psi_{\mathbf{k}_{2}}(-\mathbf{r})$
(32)

satisfies $\bar{\partial}v(\mathbf{r})=0$ [29], therefore to obey the Bloch
boundary conditions $v$ must be zero if $\mathbf{k}_{1}\neq\mathbf{k}_{2}$.

Now let us assume that we have a rank $2$ flat band discussed in the Scenario
2 in the previous section. This means that we have two wave spinors
$\psi_{\mathbf{k}_{1}}$ and $\psi_{\mathbf{k}_{2}}$ that are linearly
independent as vectors at every point in the moiré unit cell besides $n$
isolated points. By virtue of (31) we have an additional solution
$\psi_{\mathbf{k}}$. Let us take a point $\textbf{k}_{3}$ in the mBZ such that
$\mathbf{k}_{3}\neq\mathbf{k}_{1,2}$, then it is easy to check that

$\displaystyle\psi_{\mathbf{k}_{3}}\times\psi_{\mathbf{k}}=\psi_{\mathbf{k}_{1}}\left(\psi_{\mathbf{k}_{2}}\cdot\psi_{\mathbf{k}_{3}}\right)-\psi_{\mathbf{k}_{2}}\left(\psi_{\mathbf{k}_{1}}\cdot\psi_{\mathbf{k}_{3}}\right)=0\,,$
(33)

where we have used that $\psi_{\mathbf{k}_{1,2}}\cdot\psi_{\mathbf{k}_{3}}=0$.
But it means that
$\psi_{\mathbf{k}_{3}}=\mathcal{C}(\mathbf{r})\psi_{\mathbf{k}}$ for some
function $\mathcal{C}(\mathbf{r})$ (two spinors are collinear). Since
$\mathbf{k}_{3}$ can be arbitrary we conclude that each wave function in the
rank $2$ flat band is linearly dependent to each other. That is in
contradiction with our initial assumption that the flat band has rank $2$.
Therefore we conclude that for $\textbf{d}=0$ the flat bands can not realize
Scenario 2 and thus have the Chern number $C=2$.

Dirac cone on top of the flat bands. Now we prove that if we have rank $1$
flat bands then we must have a Dirac cone on top of them. Namely we prove that
we have an additional zero energy mode of the chiral (and thus also zero
energy mode of the anti-chiral) operator at an isolated point of the mBZ. We
consider a case of a single chiral flat band, when the corresponding wave
function $\psi_{\mathbf{k}}(\mathbf{r})$ has only one zero. We can pick such
$\mathbf{k}_{0}$ for which $\psi_{\textbf{k}_{0}}(0)=0$. Using discussion
below (32) one can notice that
$\psi_{\mathbf{k}_{0}}(\textbf{r})\cdot\psi_{\mathbf{k}_{0}}(-\mathbf{r})=0$.
Due to this fact the following equation

$\displaystyle\psi_{\mathbf{k}_{0}}(\textbf{r})=\hat{\phi}_{\mathbf{k}_{*}-2\mathbf{k}_{0}}(-\textbf{r})\times\psi_{\mathbf{k}_{0}}(-\textbf{r})$
(34)

has a finite solution $\hat{\phi}_{\mathbf{k}_{*}-2\mathbf{k}_{0}}$ at the
point $\mathbf{k}_{*}-2\mathbf{k}_{0}$ of the mBZ. Then if we apply the
operator $\mathcal{D}(\mathbf{r})$ to both sides of this equation we get

$\displaystyle
0=\left(\mathcal{D}(-\mathbf{r})\hat{\phi}_{\mathbf{k}_{*}-2\mathbf{k}_{0}}(-\mathbf{r})\right)\times\psi_{\mathbf{k}_{0}}(-\mathbf{r})\,,$
(35)

which leads to

$\displaystyle\mathcal{D}(\mathbf{r})\hat{\phi}_{\mathbf{k}_{*}-2\mathbf{k}_{0}}(\mathbf{r})=\mathcal{C}(\mathbf{r})\psi_{\mathbf{k}_{*}-2\mathbf{k}_{0}}(\mathbf{r})\,.$
(36)

Since $\phi_{\mathbf{k}_{*}-2\mathbf{k}_{0}}$ has a zero at some point we can
always find a function $g(\mathbf{r})$ [64] such that

$\displaystyle\left(\mathcal{C}(\mathbf{r})+\bar{\partial}g(\mathbf{r})\right)\psi_{\mathbf{k}_{*}-2\mathbf{k}_{0}}(\mathbf{r})=0$
(37)

and therefore we can construct a solution
$\phi_{\mathbf{k}_{*}-2\mathbf{k}_{0}}(\mathbf{r})=\hat{\phi}_{\mathbf{k}_{*}-2\mathbf{k}_{0}}(\mathbf{r})+g(\mathbf{r})\psi_{\mathbf{k}_{*}-2\mathbf{k}_{0}}(\mathbf{r})$
that belongs to a specific point of the mBZ and satisfies the equation

$\displaystyle\mathcal{D}(\mathbf{r})\phi_{\mathbf{k}_{*}-2\mathbf{k}_{0}}(\mathbf{r})=0\,.$
(38)

This concludes the classification of the flat bands structures for
$\mathbf{d}=0$. In summary, we can have only chiral and anti-chiral flat bands
of rank $1$ with the Chern numbers $C=1$ and $C=-1$ with a single Dirac cone
on top of them. This is schematically depicted in Fig. 5.

Figure 5: A structure of the zero energy flat bands in the case of the
displacement $\mathbf{d}=0$. Figure 6: Possible structures of the zero energy
flat bands in the case of the displacement $\mathbf{d}=\pm\textbf{d}_{0}$.

###  IV.2 Displacement $\textbf{d}=\pm\textbf{d}_{0}$

Possible topologies of the flat bands. In the case of non-zero displacement
$\textbf{d}=\pm\textbf{d}_{0}$ we can construct an anti-chiral zero mode by
taking cross product of two chiral zero modes:

$\displaystyle\chi_{\mathbf{k}_{*}-\mathbf{k}_{1}-\mathbf{k}_{2}}(\mathbf{r})=Q\left(\bar{\psi}_{\mathbf{k}_{1}}(\mathbf{r})\times\bar{\psi}_{\mathbf{k}_{2}}(\mathbf{r})\right)\,.$
(39)

Even though we can not construct an anti-chiral mode from a chiral one using
complex conjugation and inversion as in (30), the number of anti-chiral flat
bands must be equal to the number of the chiral ones, because of the particle-
hole symmetry.

Also similarly to the $\textbf{d}=0$ case one can show that the following
function

$\displaystyle
v(\mathbf{r})=\bar{\chi}_{\mathbf{k}_{1}}(-\mathbf{r})Q\psi_{\mathbf{k}_{2}}(\mathbf{r})$
(40)

satisfies $\bar{\partial}v(\mathbf{r})=0$ and thus is equal to zero if
$\mathbf{k}_{1}\neq\mathbf{k}_{2}$.

Unlike the case $\textbf{d}=0$, when $\textbf{d}=\pm\textbf{d}_{0}$ there are
no obstacles preventing the existence of rank $2$ flat bands. Let us assume
that we have a rank $2$ chiral flat band, thus we have two linearly
independent (besides $n$ points) wave spinors $\psi_{\mathbf{k}_{1}}$ and
$\psi_{\mathbf{k}_{2}}$. Using (39) we construct an anti-chiral zero mode
$\chi_{\mathbf{k}}$. The number of zeros of $\chi_{\mathbf{k}}$ is equal to
the number of points where $\psi_{\mathbf{k}_{1}}$ and $\psi_{\mathbf{k}_{2}}$
are linearly dependent. If we assume that there is some other anti-chiral zero
mode $\chi_{\mathbf{k}_{3}}$ with
$\mathbf{k}_{3}\neq\mathbf{k}_{1},\mathbf{k}_{2}$ we have

$\displaystyle\chi_{\mathbf{k}}(\textbf{r})\times\chi_{\mathbf{k}_{3}}(\textbf{r})=Q\left(Q\left(\psi_{\mathbf{k}_{1}}\times\psi_{\mathbf{k}_{1}}\right)\times\bar{\chi}_{\mathbf{k}_{3}}\right)=0\,,$
(41)

where we used (40). Therefore the anti-chiral flat band must have rank $1$.
From that it follows that we can not have at the same time rank $2$ chiral and
anti-chiral flat bands. Also one can prove that the Dirac cone on top of the
flat bands is prohibited in this case.

Finally, if we assume that we have a rank $1$ chiral flat band
$\psi_{\textbf{k}}$ and one Dirac cone on top of it with its zero mode
$\phi_{\textbf{k}_{0}}$. Then using (39) for $\phi_{\textbf{k}_{0}}$ and
$\psi_{\textbf{k}}$ we can construct a rank $1$ anti-chiral flat band. And
vice versa, using the same argument as in the previous subsection we can prove
that if we have rank $1$ chiral and anti-chiral flat bands, we must have a
single Dirac cone on top of them. We show schematically in Fig. 6 all the
possible cases for $\textbf{d}=\pm\textbf{d}_{0}$.

## Acknowledgments

We are grateful to Simon Becker, Tristan Humbert, Igor R. Klebanov, Jens
Wittsten, Mengxuan Yang for useful discussions. We would like to thank Igor R.
Klebanov for valuable comments on the draft. F.K.P. is currently a Simons
Junior Fellow at New York University and supported by a grant 855325FP from
the Simons Foundation.

## References

  * [1] R. Bistritzer and A. H. MacDonald, “Moiré bands in twisted double-layer graphene,” Proceedings of the National Academy of Sciences, vol. 108, no. 30, pp. 12233–12237, 2011.
  * [2] 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, vol. 82, p. 121407, Sep 2010.
  * [3] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature, vol. 556, no. 7699, pp. 43–50, 2018.
  * [4] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, “Correlated insulator behaviour at half-filling in magic-angle graphene superlattices,” Nature, vol. 556, no. 7699, pp. 80–84, 2018.
  * [5] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, “Tuning superconductivity in twisted bilayer graphene,” Science, vol. 363, no. 6431, pp. 1059–1064, 2019.
  * [6] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, “Origin of mott insulating behavior and superconductivity in twisted bilayer graphene,” Phys. Rev. X, vol. 8, p. 031089, Sep 2018.
  * [7] A. Thomson, S. Chatterjee, S. Sachdev, and M. S. Scheurer, “Triangular antiferromagnetism on the honeycomb lattice of twisted bilayer graphene,” Phys. Rev. B, vol. 98, p. 075109, Aug 2018.
  * [8] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, “Band structure of twisted bilayer graphene: Emergent symmetries, commensurate approximants, and wannier obstructions,” Phys. Rev. B, vol. 98, p. 085435, Aug 2018.
  * [9] F. Guinea and N. R. Walet, “Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers,” Proceedings of the National Academy of Sciences, vol. 115, no. 52, pp. 13174–13179, 2018.
  * [10] S. Carr, S. Fang, P. Jarillo-Herrero, and E. Kaxiras, “Pressure dependence of the magic twist angle in graphene superlattices,” Phys. Rev. B, vol. 98, p. 085144, Aug 2018.
  * [11] Y. Su and S.-Z. Lin, “Pairing symmetry and spontaneous vortex-antivortex lattice in superconducting twisted-bilayer graphene: Bogoliubov-de gennes approach,” Phys. Rev. B, vol. 98, p. 195101, Nov 2018.
  * [12] J. González and T. Stauber, “Kohn-luttinger superconductivity in twisted bilayer graphene,” Phys. Rev. Lett., vol. 122, p. 026801, Jan 2019.
  * [13] F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, “Topological insulators in twisted transition metal dichalcogenide homobilayers,” Phys. Rev. Lett., vol. 122, p. 086402, Feb 2019.
  * [14] D. K. Efimkin and A. H. MacDonald, “Helical network model for twisted bilayer graphene,” Phys. Rev. B, vol. 98, p. 035404, Jul 2018.
  * [15] N. F. Q. Yuan and L. Fu, “Model for the metal-insulator transition in graphene superlattices and beyond,” Phys. Rev. B, vol. 98, p. 045103, Jul 2018.
  * [16] C. Xu and L. Balents, “Topological superconductivity in twisted multilayer graphene,” Phys. Rev. Lett., vol. 121, p. 087001, Aug 2018.
  * [17] M. Ochi, M. Koshino, and K. Kuroki, “Possible correlated insulating states in magic-angle twisted bilayer graphene under strongly competing interactions,” Phys. Rev. B, vol. 98, p. 081102, Aug 2018.
  * [18] F. Wu, A. H. MacDonald, and I. Martin, “Theory of phonon-mediated superconductivity in twisted bilayer graphene,” Phys. Rev. Lett., vol. 121, p. 257001, Dec 2018.
  * [19] Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T. Senthil, “Nearly flat chern bands in moiré superlattices,” Phys. Rev. B, vol. 99, p. 075127, Feb 2019.
  * [20] J. Kang and O. Vafek, “Symmetry, maximally localized wannier states, and a low-energy model for twisted bilayer graphene narrow bands,” Phys. Rev. X, vol. 8, p. 031088, Sep 2018.
  * [21] J. M. Pizarro, M. J. Calderón, and E. Bascones, “The nature of correlations in the insulating states of twisted bilayer graphene,” Journal of Physics Communications, vol. 3, p. 035024, mar 2019.
  * [22] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, “Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene,” Phys. Rev. X, vol. 8, p. 031087, Sep 2018.
  * [23] D. M. Kennes, J. Lischner, and C. Karrasch, “Strong correlations and $d+\mathit{id}$ superconductivity in twisted bilayer graphene,” Phys. Rev. B, vol. 98, p. 241407, Dec 2018.
  * [24] H. Isobe, N. F. Q. Yuan, and L. Fu, “Unconventional superconductivity and density waves in twisted bilayer graphene,” Phys. Rev. X, vol. 8, p. 041041, Dec 2018.
  * [25] L. Rademaker and P. Mellado, “Charge-transfer insulation in twisted bilayer graphene,” Phys. Rev. B, vol. 98, p. 235158, Dec 2018.
  * [26] T. J. Peltonen, R. Ojajärvi, and T. T. Heikkilä, “Mean-field theory for superconductivity in twisted bilayer graphene,” Phys. Rev. B, vol. 98, p. 220504, Dec 2018.
  * [27] V. Kozii, H. Isobe, J. W. F. Venderbos, and L. Fu, “Nematic superconductivity stabilized by density wave fluctuations: Possible application to twisted bilayer graphene,” Phys. Rev. B, vol. 99, p. 144507, Apr 2019.
  * [28] V. Kozii, M. P. Zaletel, and N. Bultinck, “Spin-triplet superconductivity from intervalley goldstone modes in magic-angle graphene,” Phys. Rev. B, vol. 106, p. 235157, Dec 2022.
  * [29] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, “Origin of magic angles in twisted bilayer graphene,” Phys. Rev. Lett., vol. 122, p. 106405, Mar 2019\.
  * [30] P. San-Jose, J. González, and F. Guinea, “Non-abelian gauge potentials in graphene bilayers,” Phys. Rev. Lett., vol. 108, p. 216802, May 2012.
  * [31] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, “All magic angles in twisted bilayer graphene are topological,” Phys. Rev. Lett., vol. 123, p. 036401, Jul 2019.
  * [32] K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Balents, “Multiple topological transitions in twisted bilayer graphene near the first magic angle,” Phys. Rev. B, vol. 99, p. 035111, Jan 2019.
  * [33] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, “Faithful tight-binding models and fragile topology of magic-angle bilayer graphene,” Phys. Rev. B, vol. 99, p. 195455, May 2019.
  * [34] Y. Sheffer, R. Queiroz, and A. Stern, “Symmetries as the guiding principle for flattening bands of dirac fermions,” Phys. Rev. X, vol. 13, p. 021012, Apr 2023.
  * [35] E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, “Magic angle hierarchy in twisted graphene multilayers,” Phys. Rev. B, vol. 100, p. 085109, Aug 2019.
  * [36] C. Mora, N. Regnault, and B. A. Bernevig, “Flatbands and perfect metal in trilayer moiré graphene,” Phys. Rev. Lett., vol. 123, p. 026402, Jul 2019\.
  * [37] T. Cea, N. R. Walet, and F. Guinea, “Twists and the electronic structure of graphitic materials,” Nano Letters, vol. 19, pp. 8683–8689, 12 2019.
  * [38] Z. Zhu, S. Carr, D. Massatt, M. Luskin, and E. Kaxiras, “Twisted trilayer graphene: A precisely tunable platform for correlated electrons,” Phys. Rev. Lett., vol. 125, p. 116404, Sep 2020.
  * [39] Y. Mao, D. Guerci, and C. Mora, “Supermoiré low-energy effective theory of twisted trilayer graphene,” 2023.
  * [40] X. Lin, C. Li, K. Su, and J. Ni, “Energetic stability and spatial inhomogeneity in the local electronic structure of relaxed twisted trilayer graphene,” Physical Review B, vol. 106, no. 7, p. 075423, 2022.
  * [41] Z. Ma, S. Li, M. Lu, D.-H. Xu, J.-H. Gao, and X. Xie, “Doubled moiré flat bands in double-twisted few-layer graphite,” Science China Physics, Mechanics & Astronomy, vol. 66, no. 2, p. 227211, 2023.
  * [42] M. Liang, M.-M. Xiao, Z. Ma, and J.-H. Gao, “Moiré band structures of the double twisted few-layer graphene,” Phys. Rev. B, vol. 105, p. 195422, May 2022.
  * [43] J. Y. Lee, E. Khalaf, S. Liu, X. Liu, Z. Hao, P. Kim, and A. Vishwanath, “Theory of correlated insulating behaviour and spin-triplet superconductivity in twisted double bilayer graphene,” Nature Communications, vol. 10, no. 1, p. 5333, 2019.
  * [44] P. J. Ledwith, A. Vishwanath, and E. Khalaf, “Family of ideal chern flatbands with arbitrary chern number in chiral twisted graphene multilayers,” Phys. Rev. Lett., vol. 128, p. 176404, Apr 2022.
  * [45] P. J. Ledwith, E. Khalaf, Z. Zhu, S. Carr, E. Kaxiras, and A. Vishwanath, “Tb or not tb? contrasting properties of twisted bilayer graphene and the alternating twist $n$-layer structures ($n=3,4,5,\backslash$),” arXiv preprint arXiv:2111.11060, 2021.
  * [46] J. Wang and Z. Liu, “Hierarchy of ideal flatbands in chiral twisted multilayer graphene models,” Phys. Rev. Lett., vol. 128, p. 176403, Apr 2022.
  * [47] S. Zhang, B. Xie, Q. Wu, J. Liu, and O. V. Yazyev, “Chiral decomposition of twisted graphene multilayers with arbitrary stacking,” Nano Letters, 03 2023.
  * [48] M. Yang, “Flat bands and high chern numbers in twisted multilayer graphene,” arXiv:2303.00103, 2023.
  * [49] X. Zhang, K.-T. Tsai, Z. Zhu, W. Ren, Y. Luo, S. Carr, M. Luskin, E. Kaxiras, and K. Wang, “Correlated insulating states and transport signature of superconductivity in twisted trilayer graphene superlattices,” Phys. Rev. Lett., vol. 127, p. 166802, Oct 2021.
  * [50] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, “Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene,” Nature, vol. 590, no. 7845, pp. 249–255, 2021.
  * [51] Z. Hao, A. M. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, “Electric field-tunable superconductivity in alternating-twist magic-angle trilayer graphene,” Science, vol. 371, no. 6534, pp. 1133–1138, 2021.
  * [52] X. Liu, N. J. Zhang, K. Watanabe, T. Taniguchi, and J. I. A. Li, “Isospin order in superconducting magic-angle twisted trilayer graphene,” Nature Physics, vol. 18, no. 5, pp. 522–527, 2022.
  * [53] S. Turkel, J. Swann, Z. Zhu, M. Christos, K. Watanabe, T. Taniguchi, S. Sachdev, M. S. Scheurer, E. Kaxiras, C. R. Dean, and A. N. Pasupathy, “Orderly disorder in magic-angle twisted trilayer graphene,” Science, vol. 376, no. 6589, pp. 193–199, 2022.
  * [54] A. Uri, S. C. de la Barrera, M. T. Randeria, D. Rodan-Legrain, T. Devakul, P. J. Crowley, N. Paul, K. Watanabe, T. Taniguchi, R. Lifshitz, et al., “Superconductivity and strong interactions in a tunable moire quasiperiodic crystal,” arXiv:2302.00686, 2023.
  * [55] D. Călugăru, F. Xie, Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, “Twisted symmetric trilayer graphene: Single-particle and many-body hamiltonians and hidden nonlocal symmetries of trilayer moiré systems with and without displacement field,” Phys. Rev. B, vol. 103, p. 195411, May 2021.
  * [56] F. Xie, N. Regnault, D. Călugăru, B. A. Bernevig, and B. Lian, “Twisted symmetric trilayer graphene. ii. projected hartree-fock study,” Phys. Rev. B, vol. 104, p. 115167, Sep 2021.
  * [57] M. Christos, S. Sachdev, and M. S. Scheurer, “Correlated insulators, semimetals, and superconductivity in twisted trilayer graphene,” Phys. Rev. X, vol. 12, p. 021018, Apr 2022.
  * [58] F. K. Popov and G. Tarnopolsky, “Magic angles in equal-twist trilayer graphene,” 2023.
  * [59] D. Guerci, Y. Mao, and C. Mora, “Chern mosaic and ideal flat bands in equal-twist trilayer graphene,” arXiv preprint arXiv:2305.03702, 2023.
  * [60] T. Devakul, P. J. Ledwith, L.-Q. Xia, A. Uri, S. de la Barrera, P. Jarillo-Herrero, and L. Fu, “Magic-angle helical trilayer graphene,” arXiv preprint arXiv:2305.03031, 2023.
  * [61] N. Nakatsuji, T. Kawakami, and M. Koshino, “Multi-scale lattice relaxation in chiral twisted trilayer graphenes,” arXiv preprint arXiv:2305.13155, 2023\.
  * [62] S. Becker, M. Embree, J. Wittsten, and M. Zworski, “Spectral characterization of magic angles in twisted bilayer graphene,” Physical Review B, vol. 103, no. 16, p. 165113, 2021.
  * [63] S. Becker, M. Embree, J. Wittsten, and M. Zworski, “Mathematics of magic angles in a model of twisted bilayer graphene,” Probability and Mathematical Physics, vol. 3, no. 1, pp. 69–103, 2022.
  * [64] F. K. Popov and A. Milekhin, “Hidden wave function of twisted bilayer graphene: The flat band as a landau level,” Phys. Rev. B, vol. 103, p. 155150, Apr 2021.


This essay received an Honorable Mention in the 2018 Essay Competition of the
Gravity Research Foundation

# Fermat’s principle in black-hole spacetimes

Shahar Hod  The Ruppin Academic Center, Emeq Hefer 40250, Israel  The Hadassah
Institute, Jerusalem 91010, Israel

###### Abstract

Black-hole spacetimes are known to possess closed light rings. We here present
a remarkably compact theorem which reveals the physically intriguing fact that
these unique null circular geodesics provide the fastest way, as measured by
asymptotic observers, to circle around spinning Kerr black holes.

Email<EMAIL_ADDRESS>

Fermat’s principle, also known as the principle of least time Lip ; Notewik ,
asserts that among all possible null trajectories, the path taken by a ray of
light between two given points A and B in a flat spacetime geometry is the
path that minimizes the traveling time $T_{A\to B}$. This remarkably elegant
principle implies, in particular, that the unique null trajectory taken by a
ray of light between two given points is generally distinct from the straight
line trajectory which minimizes the spatial distance $d_{AB}$ between these
points.

In the present Essay we would like to highlight an intriguing and closely
related physical phenomenon which characterizes curved spacetime geometries.
In particular, we here raise the physically interesting question: Among all
possible closed paths that circle around a black hole in a curved spacetime,
which path provides the fastest way, as measured by asymptotic observers, to
circle the central black hole?

We first note that, in flat spacetimes, the characteristic orbital period
$T_{\odot}$ of a test particle that moves around a spatially compact object of
radius $R$ is trivially bounded from below by the compact relation Noteunit

$T_{\odot}\geq T^{\text{flat}}_{\odot\text{min}}=2\pi R\ ,$ (1)

where the equality sign in (1) is attained by massless particles that circle
the compact object on the shortest possible (tangential) trajectory with
$r_{\text{fast}}=R$.

It should be emphasized, however, that the simple lower bound (1) is not valid
in realistic curved spacetimes. In particular, it does not take into account
the important time dilation (red-shift) effect which is caused by the
gravitational field of the central compact object Chan . In addition, the
flat-space relation (1) does not take into account the well-known phenomenon
of dragging of inertial frames by spinning compact objects in curved
spacetimes Chan .

As we shall explicitly show below, due to the influences of these two
interesting physical effects, the shortest possible orbital period $T_{\odot}$
of a test particle around a central compact object, as measured by asymptotic
observers, is larger than the naive flat-space estimate (1). In particular, we
shall prove that, in generic curved spacetimes, the unique circular trajectory
$r=r_{\text{fast}}$ that minimizes the traveling time $T_{\odot}$ around a
central Kerr black hole is distinct from the tangential trajectory with
$r=r_{\text{short}}$ which could minimize the traveling distance around the
spinning black hole.

The fastest circular orbit around a spinning Kerr black hole.— We shall
analyze the physical and mathematical properties of equatorial circular
trajectories around spinning Kerr black holes. In Boyer-Lindquist coordinates
$(t,r,\theta,\phi)$, the asymptotically flat black-hole spacetime can be
described by the curved line element Chan

$\displaystyle ds^{2}=-{{\Delta}\over{\rho^{2}}}(dt-a\sin^{2}\theta
d\phi)^{2}+{{\rho^{2}}\over{\Delta}}dr^{2}+\rho^{2}d\theta^{2}+{{\sin^{2}\theta}\over{\rho^{2}}}\big{[}adt-(r^{2}+a^{2})d\phi\big{]}^{2}\
,$ (2)

where $M$ is the black-hole mass, $J\equiv Ma$ is its angular momentum,
$\Delta\equiv r^{2}-2Mr+a^{2}$, and $\rho^{2}\equiv
r^{2}+a^{2}\cos^{2}\theta$. The black-hole (event and inner) horizons are
determined by the spatial zeros of the metric function $\Delta$:

$r_{\pm}=M\pm(M^{2}-a^{2})^{1/2}\ .$ (3)

We would like to identify the unique circular trajectory which minimizes the
orbital period $T_{\odot}$, as measured by asymptotic observers, around the
central black hole. We shall therefore assume the relation $v/c\to 1^{-}$ for
the velocity of the orbiting test particle Notegng . The corresponding radius-
dependent orbital periods $T_{\odot}(r)$ of the test particles can easily be
obtained from the characteristic black-hole curved line element (2) with
$ds=dr=d\theta=0$ and $\Delta\phi=\pm 2\pi$ Notecocoun ; Noteco ; Hod1 . This
yields the compact functional relation

$T_{\odot}(r)=2\pi\cdot{{\sqrt{r^{2}-2Mr+a^{2}}-{{2Ma}\over{r}}}\over{1-{{2M}\over{r}}}}\
$ (4)

for the orbital periods of co-rotating test particles around the central
spinning black hole.

The physically interesting co-rotating circular orbit with
$r=r_{\text{fast}}$, which is characterized by the shortest possible orbital
period $T_{\odot\text{min}}=\min_{r}\\{T_{\odot}(r)\\}$ around the central
spinning black hole, is determined by the functional relation
$dT(r=r_{\text{fast}})/dr=0$. This yields the characteristic algebraic
equation

$r^{2}-3Mr+2a^{2}+2a\sqrt{r^{2}-2Mr+a^{2}}=0\ \ \ \ \text{for}\ \ \ \
r=r_{\text{fast}}\ .$ (5)

Remarkably, this equation can be solved analytically to yield the simple
functional relation

$r_{\text{fast}}=2M\cdot\big{\\{}1+\cos[{2\over 3}\cos^{-1}(-{a/M})]\big{\\}}\
.$ (6)

for the unique orbital radius $r_{\text{fast}}(M,a)$ which characterizes the
fastest co-rotating circular trajectory (the closed circular path with the
shortest possible orbital period) around the central spinning Kerr black hole.

What we find most intriguing is the fact that the spin-dependent radii
$r_{\text{fast}}(M,a)$ of the fastest circular trajectories, as given by the
functional expression (6), exactly coincide with the corresponding radii
$r_{\gamma}(M,a)$ of the null circular geodesics Bar  which characterize the
spinning Kerr black-hole spacetimes. One therefore concludes that co-rotating
null circular geodesics (closed light rings) provide the fastest way, as
measured by asymptotic observers, to circle around generic Kerr black holes.

It is physically interesting to define the dimensionless ratio [see Eqs. (1)
and (3)]

$\Theta({\bar{a}})\equiv{{T_{\odot\text{min}}}\over{2\pi r_{+}}}\ \ \ \ ;\ \ \
\ {\bar{a}}\equiv a/M\ ,$ (7)

which characterizes the unique closed circular trajectories [with
$r=r_{\text{fast}}({\bar{a}})$] that minimize the orbital periods around the
central spinning black holes. As emphasized above, a naive flat-space
calculation predicts the relation $\Theta^{\text{flat}}_{\text{min}}\equiv
T^{\text{flat}}_{\odot\text{min}}/2\pi R=1$ [see Eq. (1)]. However,
substituting Eqs. (4) and (6) into (7), one finds the characteristic
inequality Notetheta

$\Theta^{\text{Kerr}}({\bar{a}})>1\ $ (8)

for all Kerr black-hole spacetimes in the physically allowed regime
${\bar{a}}\in[0,1]$. In particular, the dimensionless function
$\Theta^{\text{Kerr}}({\bar{a}})$ exhibits a non-trivial (non-monotonic)
functional dependence on the dimensionless black-hole rotation parameter
${\bar{a}}$ with the property Noteexpan

$\text{min}_{\bar{a}}\big{\\{}\Theta^{\text{Kerr}}({\bar{a}})\big{\\}}\simeq
2-{{3(13-7\sqrt{3})}\over{88}}\ \ \ \ \text{at}\ \ \ \
{\bar{a}}^{\text{Kerr}}_{\text{min}}\simeq 1-{{126-45\sqrt{3}}\over{1936}}\ .$
(9)

Summary.— Fermat’s principle asserts that, in a flat spacetime geometry, the
path taken by a ray of light is unique in the sense that it represents the
spatial trajectory with the shortest possible traveling time between two given
points Lip ; Notewik . This intriguing principle implies, in particular, that
the paths taken by light rays are generally distinct from the straight line
trajectories which could minimize the traveling distances between two given
points.

In the present short Essay we have highlighted an intriguing and closely
related phenomenon in curved black-hole spacetimes. In particular, we have
raised the physically interesting question: Among all possible trajectories
that circle around a spinning Kerr black hole, which closed trajectory
provides the fastest way, as measured by asymptotic observers, to circle the
central black hole?

Our compact theorem has revealed the physically intriguing fact that the
equatorial null circular geodesics (closed light rings), which characterize
the curved black-hole spacetimes, provide the fastest way to circle around
spinning Kerr black holes. In particular, we have explicitly proved that, in
analogy with the Fermat principle in flat spacetime geometries, the unique
curved trajectories $r=r_{\text{fast}}(M,a)$ [see Eq. (6)] which minimize the
traveling times $T_{\odot}$ of test particles around central black holes are
distinct from the tangential trajectories $r=r_{+}(M,a)$ [see Eq. (3)] which
could minimize the traveling distances around the black holes.

ACKNOWLEDGMENTS

This research is supported by the Carmel Science Foundation. I would like to
thank Yael Oren, Arbel M. Ongo, Ayelet B. Lata, and Alona B. Tea for helpful
discussions.

## References

  * (1) A. Lipson, S. G. Lipson, H. Lipson, Optical Physics 4th Edition, Cambridge University Press (2010).
  * (2) See also https://en.wikipedia.org/wiki/Fermat$\%$27s_principle .
  * (3) We use natural units in which $G=c=1$.
  * (4) S. Chandrasekhar, The Mathematical Theory of Black Holes, (Oxford University Press, New York, 1983).
  * (5) In the present essay we consider geodesic as well as non-geodesic trajectories of test particles around the central Kerr black holes. It is well known that the null circular geodesic of the black-hole spacetime, on which massless particles that move exactly at the speed of light can circle the central black hole, is characterized by a well defined radius $r_{\gamma}=r_{\gamma}(M,a)$. It is important to stress the fact that, by using man-made rockets which are based on non-gravitational forces, massive particles can also circle the central black hole on non-geodesic trajectories with orbital velocities that, in principle, may approach arbitrarily close to the speed of light and with orbital radii that may differ from the unique radius $r_{\gamma}(M,a)$ of the black-hole null circular geodesic.
  * (6) Here the upper/lower signs correspond respectively to co-rotating/counter-rotating trajectories of the test particles around the central spinning black holes.
  * (7) In the present essay we are interested in circular trajectories that minimize the orbital periods $T_{\odot}$ around the central spinning black holes as measured by asymptotic observers. We shall therefore focus on co-rotating circular orbits.
  * (8) S. Hod, Phys. Rev. D 84, 104024 (2011).
  * (9) See equation (2.18) of J. M. Bardeen, W. H. Press and S. A. Teukolsky, Astrophys. Jour. 178, 347 (1972).
  * (10) In particular, one finds $\Theta^{\text{Kerr}}({\bar{a}}\to 0)\to 3\sqrt{3}/2$ in the limit of non-rotating Schwarzschild black holes and $\Theta^{\text{Kerr}}({\bar{a}}\to 1)\to 2$ in the limit of maximally spinning (extremal) Kerr black holes.
  * (11) Here we have used the series expansion $\Theta^{\text{Kerr}}({\bar{a}})=2+\sqrt{2}(\sqrt{3}-2)\cdot\sqrt{\epsilon}+(14/3-2\sqrt{3})\cdot\epsilon+O(\epsilon^{3/2})$ in the near-extremal $0\leq\epsilon\equiv 1-{\bar{a}}\ll 1$ regime.