# Nearly flat bands in twisted triple bilayer graphene

Jiseon Shin  Department of Physics, University of Seoul, Seoul 02504, Korea
Bheema Lingam Chittari  Department of Physical Sciences, Indian Institute of
Science Education and Research Kolkata, Mohanpur 741246, West Bengal, India
Yunsu Jang  Department of Physics and Astronomy, Seoul National University,
Seoul 08826, Republic of Korea  Hongki Min  Department of Physics and
Astronomy, Seoul National University, Seoul 08826, Republic of Korea  Jeil
Jung  Department of Physics, University of Seoul, Seoul 02504, Korea
Department of Smart Cities, University of Seoul, Seoul 02504, Korea

###### Abstract

We investigate the electronic structure of alternating-twist triple Bernal-
stacked bilayer graphene (t3BG) as a function of interlayer coupling $\omega$,
twist angle $\theta$, interlayer potential difference $\Delta$, and top-bottom
bilayers sliding vector $\bm{\tau}$ for three possible configurations
AB/AB/AB, AB/BA/AB, and AB/AB/BA. The parabolic low-energy band dispersions in
a Bernal-stacked bilayer and gap-opening through a finite interlayer potential
difference $\Delta$ allows the flattening of bands in t3BG down to $\sim 20$
meV for twist angles $\theta\lesssim 2^{\circ}$ regardless of the stacking
types. The easier isolation of the flat bands and associated reduction of
Coulomb screening thanks to the intrinsic gaps of bilayer graphene for finite
$\Delta$ facilitate the formation of correlation-driven gaps when it is
compared to the metallic phases of twisted trilayer graphene under electric
fields. We obtain the stacking dependent Coulomb energy versus bandwidth
$U/W\gtrsim 1$ ratios in the $\theta$ and $\Delta$ parameter space. We also
present the expected $K$-valley Chern numbers for the lowest-energy nearly
flat bands.

##  I Introduction

Figure 1:  (Color online) Schematic diagrams for the starting stacking
configurations before a twist of (a) type I(-,0,+), (b) type II(-,0,+), and
(c) type III(-,0,+). (d) The moire patterns at $\theta=1.5^{\circ}$ in real
space for type I0 and the location of the corresponding three local
commensurate stackings a, b, and c are indicated. (e) Schematic diagrams for
four different cases where the couplings between bilayers $\omega$ are
considered as a perturbation. The three parabolic bands have the intrinsic
bilayer gap ($\Delta$) and the two parabolic bands at $\tilde{K}^{\prime}$
have the band offset (2$\Delta^{\prime}$+$\Delta$) due to an external electric
field in the $z$-direction. The three parabolic dispersive bands for the
second ($\omega\neq 0$, $\Delta^{(\prime)}=0$) and the fourth ($\omega\neq 0$,
$\Delta^{(\prime)}\neq 0$) cases are plotted with dotted lines as a reference
in Fig. 2. When the tunneling is considered perturbatively, the interband
hybridizations have $\Delta_{c}\approx 7.9$ meV and
$\Delta_{c}^{\prime}\approx 5.0$ meV at $\theta=1.5^{\circ}$ for $\omega=0.01$
eV and $\Delta=\Delta^{\prime}=0.01$ eV.

Intense research activity on the physics of twisted structures has unfolded
since experimental observations of correlated phases and flat band
superconductivity near its magic-angle in twisted bilayer graphene (t2G) Kim
_et al._ (2017); Cao  _et al._ (2018a, b); Yankowitz  _et al._ (2019); Cao
_et al._ (2020a) and twisted trilayer graphene (t3G) Park  _et al._ (2021a);
Hao  _et al._ (2021). Many conceivable combinations of twisted graphene
systems beyond t2G have been studied, including the twisted double bilayer
graphene (t2BG) Lee  _et al._ (2019); Chebrolu  _et al._ (2019); Choi and Choi
(2019); Koshino (2019); Burg  _et al._ (2019); Cao  _et al._ (2020b); Shen
_et al._ (2020); Liu  _et al._ (2020) that shows electric-field-tunable
lowest-energy bands up to two times narrower than that of t2G for similar
system parameters Lee  _et al._ (2019); Chebrolu  _et al._ (2019); Choi and
Choi (2019); Koshino (2019); Burg  _et al._ (2019); Cao  _et al._ (2020b);
Shen  _et al._ (2020); Liu  _et al._ (2020), and the ABC-stacked trilayer
graphene on hexagonal boron nitride Chittari  _et al._ (2019) that shows
electric field tunable bandwidths and band isolation by opening a gap at the
charge neutrality point. Other trilayer systems that have been studied
recently include the twisted mono-bilayers Suárez Morell  _et al._ (2013); Li
_et al._ (2019); Szendrő  _et al._ (2020); Carr  _et al._ (2020); Ma  _et al._
(2021); Xu  _et al._ (2021); Chen  _et al._ (2021); Polshyn  _et al._ (2020);
Park  _et al._ (2020); Lei  _et al._ (2020); Wu  _et al._ (2021), the
consecutive-twist Mora  _et al._ (2019) and alternating twisted trilayer
graphene (t3G) Park  _et al._ (2021a); Zuo  _et al._ (2018); Khalaf  _et al._
(2019); Carr  _et al._ (2020); Li  _et al._ (2019); Lopez-Bezanilla and Lado
(2020); Lei  _et al._ (2020); Hao  _et al._ (2021); Călugăru  _et al._ (2021);
Shin  _et al._ (2021); Mora  _et al._ (2019).

This work expands the ongoing effort in search of graphene-based flat band
systems by proposing that the electronic structure of twisted triple bilayer
graphene (t3BG) formed by three Bernal-stacked bilayer graphene with
alternating twists can offer advantages thanks to the electric field tunable
intrinsic gaps of each bilayer graphene Min  _et al._ (2007) and the
possibility of achieving narrow bandwidths at larger twist angles due to the
effective increase of the interlayer tunneling strength expected in
alternating twist triple layer systems Lopez-Bezanilla and Lado (2020). We
perform numerical calculations for sufficiently small twist angles $\theta\leq
2.5^{\circ}$ in search of nearly flat low-energy bands by using continuum
moire bands models used in earlier works Bistritzer and MacDonald (2011); Jung
_et al._ (2014). A total of 9 different starting stacking combinations are
considered, to account for the three different possible stacking orientations
and three sliding positions of the top bilayer with respect to the bottom
bilayer, and we consider the remote hopping terms in the bilayer graphene Jung
and MacDonald (2014) for a more accurate description of the electronic
structure. Our phase diagram calculations in the parameter space of twist
angle $\theta$ and perpendicular electric fields induced by interlayer
potential differences $\Delta$ indicate that narrow bandwidths below $\sim 20$
meV can be achieved for twist angles below $\theta\lesssim 2^{\circ}$
regardless of stacking types. A perpendicular electric field can be used as a
control knob of the bandwidths and band isolation that can facilitate the
onset of correlation gaps and valley Chern bands.

Our manuscript is organized as follows. In Sec II, we give a detailed
description of the stacking configurations at each moire interface and the
associated model Hamiltonian. We devote Sec. III to describe how the
electronic properties of t3BG change as a function of the twist angles,
interlayer couplings, starting stackings, sliding of the top bilayer, and
electric fields. In Sec. IV, we explore a wide parameter space of the
interlayer potential differences and twist angles in search of large Coulomb
energy versus bandwidth ($U/W>1$) regime. In Sec. V we describe the valley
Chern number phase diagrams that characterizes the topological properties of
the lowest-energy electron and hole bands in t3BG. We close the paper in Sec.
VI with the conclusions.

##  II Model

Figure 2:  (Color online) Energy dispersions and density of states (DOS) of
t3BG for type I0 at $\theta=2^{\circ}$ and $\theta=1.5^{\circ}$ for (a) rigid
($\omega=\omega^{\prime}$) and (b) relaxed ($\omega\neq\omega^{\prime}$)
lattice where $\Delta~{}(\Delta^{\prime})$ represents the displacement field
within (between) bilayers. The parabolic bands due to decoupled bilayers in
Fig. 1 (e) which are indicated by the purple and yellow dotted lines in this
figure are superposed on the energy dispersions. For zero interlayer
displacement field, the effective mass of one of the two parabolic dispersions
at $\tilde{K}^{\prime}$ remains constant regardless of the twist angle, while
the other one at $\tilde{K}^{\prime}$ and the one at $\tilde{K}$ hybridize
with each other and generate the nearly flat bands at the magic-angle. For
finite interlayer displacement field, the three parabolic dispersions are
depicted with pink, purple, and yellow dotted lines and note that the
effective mass of the two parabolic dispersions at $\tilde{K}^{\prime}$ varies
upon turning on the interlayer displacement field. Bandwidths of the lowest-
energy valence (red) and conduction (blue) bands in type I0 as a function of
the twist angle $\theta$ for three different interlayer tunnelings
$\omega=0.12,0.15,0.18$ eV when (c) $\omega=\omega^{\prime}$ and (d)
$\omega\neq\omega^{\prime}$, and as a function of $\omega$ for three different
twist angles $\theta=1.0^{\circ},1.5^{\circ},2^{\circ}$ when (e)
$\omega=\omega^{\prime}$, (f) $\omega\neq\omega^{\prime}$ such that
$\omega^{\prime}=C_{2}\omega^{2}+C_{1}\omega+C_{0}$ where $C_{2}=-0.5506$
eV-1, $C_{1}=1.036$, and $C_{0}=-0.02245$ eV [see Sec. II].  Figure 3:  (Color
online) Bandwidth, secondary bandgap, and primary bandgap of the lowest-energy
valence (V) and conduction (C) bands of type I0 for
$\Delta^{\prime}=\Delta=0.0$ eV for the (a) minimal and (b) remote hopping
models, and $\Delta^{\prime}=\Delta=0.02$ eV for the (c) minimal and (d)
remote hopping model as a function of the pressure and the twist angle where
$\Delta~{}(\Delta^{\prime})$ represents the gate potential difference within
(between) bilayers. The solid lines in the first and the third columns
represent contours of the bandwidth of the lowest-energy valence and the
conduction bands. The solid lines in the second, fourth and fifth columns
represent the contour at the zero gaps. The cross-section at
$\theta=1.5^{\circ}$ is presented on the rightmost column. The dashed lines
indicate zero pressure.

We consider three vertically stacked graphene bilayers with alternating twists
such that the top and bottom bilayers are rotated by $-\theta/2$ (clockwise),
while the middle bilayer is rotated by $+\theta/2$ (counterclockwise) as shown
in Fig. 1. We note that this twist configuration leading to commensurate moire
patterns is favored in experimental devices Park  _et al._ (2021a); Hao  _et
al._ (2021); Park  _et al._ (2021b). We distinguish three starting stacking
configurations AB/AB/AB, AB/BA/AB, and AB/AB/BA and call them type I0, II0,
and III0 respectively. Subsequently, we take into account slidings of the top
bilayer with respect to the bottom bilayer by $0$, $\pm a/\sqrt{3}$ along the
$y$-axis parallel to the carbon dimers in the unit cell which we express
compactly as $\bm{\tau}_{0,\pm}$. We note that $\pm 3\bm{\tau}_{\pm}$
translation brings the system back to the original configuration. Also, a
sliding by $\bm{\tau}_{+}$ ($\bm{\tau}_{-}$) is equivalent to
$-2\bm{\tau}_{-}$ ($+2\bm{\tau}_{+}$). Therefore, we consider 9 in total
possible starting stacking geometries that we name as I(0,+,-), II(0,+,-), and
III(0,+,-) as shown in Figs. 1(a)-(c) where a translation of the top bilayer
by $\bm{\tau}_{0}$ or $\bm{\tau}_{\pm}$ is denoted by $0$ or $\pm$ in the
subscript of the stacking types, I, II or III. The moire patterns for type I0
are illustrated in Fig. 1(d) where we indicate the three representative local
commensurate stacking configurations by the letters a, b, c. All 9 starting
stacking geometries in Figs. 1(a)-(c) correspond to the local commensurate
stacking at location a in the moire patterns where we put the rotation axis.
We will show in Sec. III how the electronic band structures and local density
of states are modified depending on the initial stacking geometries.

Figure 4:  (Color online) Schematic diagrams for the atomic configurations of
local commensurate stackings at a, b, and c on the upper row, the band
structures and the associated DOS at $\theta=1.5^{\circ}$ for the minimal
model (red dotted lines for bands and DOS) and for the remote model (black
solid lines for bands, and blue filled area for DOS) with LDOS in the
hexagonal frames at the van Hove singularities denoted by the red (minimal
model) and the blue (remote model) arrows, when $\Delta^{\prime}=\Delta=0.0$
eV in the middle row, and when $\Delta^{\prime}=\Delta=0.02$ eV in the lower
row for (a) type I-, (b) type I0, (c) type I+.

The Hamiltonian of our model can be expressed concisely as

$H_{\textrm{t3BG}}(\theta)=\left(\begin{array}[]{cccccc}h_{l,b}^{-}+\bar{V}_{1}&t_{s}^{-}&0&0&0&0\\\
t_{s}^{-\dagger}&h_{u,b}^{-}+\bar{V}_{2}&T_{1}(\bm{r})&0&0&0\\\
0&T_{1}^{\dagger}(\bm{r})&h_{l,m}^{+}+\bar{V}_{3}&t_{s}^{+}&0&0\\\
0&0&t_{s}^{+\dagger}&h_{u,m}^{+}+\bar{V}_{4}&T_{2}(\bm{r})&0\\\
0&0&0&T_{2}^{\dagger}(\bm{r})&h_{l,t}^{-}+\bar{V}_{5}&t_{s}^{-}\\\
0&0&0&0&t_{s}^{-\dagger}&h_{u,t}^{-}+\bar{V}_{6}\\\ \end{array}\right),$ (1)

where $h^{\pm}_{r,(b,m,t)}=h(\pm\theta/2)+(\mathbb{1}+rm_{s}\sigma_{z})\hskip
1.42271pt\delta/2$ for the upper ($r=+1$) or the lower ($r=-1$) layer within
an AB-stacked ($m_{s={\rm AB}}=1$) or BA-stacked ($m_{s={\rm BA}}=-1$) bilayer
located at the bottom (b), middle (m), or top (t). The first term represents a
$2\times 2$ Dirac Hamiltonian
$h(\theta=0)={\upsilon}_{F}\bm{p}\cdot\bm{\sigma}$ rotated by $\pm\theta/2$
which is given by

$h(\pm\theta/2)=D^{\dagger}(\pm\theta/2)~{}h(\theta=0)~{}D(\pm\theta/2),$ (2)

where $D(\phi)=e^{-i\phi\sigma_{z}/2}$ is the general form of the spin $1/2$
rotation operator with the $z$-component Pauli matrix $\sigma_{z}$, and
$\upsilon_{F}=\sqrt{3}|t_{0}|a/2\hbar$, where we choose $t_{0}=-3.1$ eV. The
second term describes the intrinsic sublattice potential $\delta=0.015$ eV due
to the Bernal-stacked bilayer emerging at higher-energy dimers at the upper or
lower layers Jung and MacDonald (2014).

Figure 5:  (Color online) Same as Fig. 4 for type II.

The $2\times 2$ matrices $\bar{V}_{i}=V_{i}\mathbb{1}$ ($i=1,\cdots,6$) on the
diagonal are the potential energies on the $i^{th}$ layer which can be
expressed in terms of the interlayer potential differences within a bilayer
$\Delta$ and the potential differences between bilayers $\Delta^{\prime}$.
When we assume that the potential $V_{i}$ drops monotonically at the same rate
neglecting screening effects, the $V_{i}$ values satisfy
$V_{1}=3\Delta/2+\Delta^{\prime}$, $V_{2}=\Delta/2+\Delta^{\prime}$,
$V_{3}=\Delta/2$, $V_{4}=-V_{3}$, $V_{5}=-V_{2}$, and $V_{6}=-V_{1}$.

The interlayer tunneling terms within each bilayer are given by 2 $\times$ 2
matrix $t^{\pm}_{s}$ (where $s=$ AB or BA) consisting of (1) the vertical
hopping term in the Bernal-stacked bilayer $t_{1}=0.361$ eV originating from
$t_{1}(\sigma_{x}-im_{s}\sigma_{y})/2$ for either AB- or BA-stacked bilayer,
and (2) the two additional hopping terms, $v_{3}$ and $v_{4}$, responsible for
the trigonal warping and electron-hole asymmetry, respectively, as proposed in
Ref. Jung and MacDonald (2014), where $v_{i}=\sqrt{3}|t_{i}|a/2\hbar$,
$t_{3}=0.283$ eV and $t_{4}=0.183$ eV. The explicit form of the matrix for an
AB-stacked bilayer $t^{\pm}_{AB}$, as an example, is given as follows,

$t^{\pm}_{AB}=\left(\begin{array}[]{cc}-v_{4}\pi^{\pm\dagger}&-v_{3}\pi^{\pm}\\\
t_{1}&-v_{4}\pi^{\pm\dagger}\\\ \end{array}\right),$ (3)

resulting in $t^{\pm}_{BA}=t^{\pm\dagger}_{AB}$. Here, $\pi^{\pm}$ represents
$p_{x}+ip_{y}$ with an implicit phase of $\pm\theta/2$ due to the twist. In
the minimal model the remote hopping terms are set to zero, i.e.
$\delta=v_{3}=v_{4}=0$.

Figure 6:  (Color online) Same as Fig. 4 for type III.

The interlayer tunneling at the twisted interface is captured by Bistritzer
and MacDonald (2011); Jung  _et al._ (2014)

$T_{k}(\bm{r})=\sum_{j=0,\pm}e^{-im_{k}\bm{q}_{j}\cdot\bm{r}}T^{j}_{ll^{\prime}}.$
(4)

where $m_{k}=(-1)^{k}$, and $\bm{q}_{0}$, $\bm{q}_{\pm}$ are given as
$\bm{q}_{0}=\theta k_{D}(0,-1)$, $\bm{q}_{\pm}=\theta
k_{D}(\pm\sqrt{3}/2,1/2)$ in the small angle approximation. Here,
$k_{D}=4\pi/3a$ is equal to the length of a side of the Brillouin zone of a
monolayer graphene where $a=2.461\textrm{\AA}$.

The interlayer tunneling at moire interfaces can be expressed in terms of the
horizontal displacement vector $\bm{d}$ which measures the horizontal
direction and distance between two atoms on the same sublattice – one atom
from the upper layer and the other from the lower layer. The coplanar
displacement between two atoms of the same sublattice from the different
layers, for example, is given as $\bm{d}=(d_{x},d_{y})=(0,0)$ for AA′ and BB′
stackings, $(0,+a/\sqrt{3})$ for AB′ stacking, and $(0,-a/\sqrt{3})$ for BA′
stacking at the moire interfaces. Then, the interlayer tunneling is explicitly
expressed as $T_{ll’}=\omega_{ll^{\prime}}\exp(ij\bm{G}_{j}\cdot\bm{d})$ for
$j=0,\pm 1$ corresponding to $\bm{G}_{0}=(0,0)$ and
$\bm{G}_{\pm}=k_{D}(-3/2,\pm\sqrt{3}/2)$ with sublattice indices $l$ and
$l^{\prime}$, yielding

$T^{0}=\left(\begin{array}[]{cc}\omega^{\prime}&\omega\\\
\omega&\omega^{\prime}\\\ \end{array}\right),\hskip
8.5359ptT^{\pm}=\left(\begin{array}[]{cc}\omega^{\prime}&\omega e^{\mp
i2\pi/3}\\\ \omega e^{\pm i2\pi/3}&\omega^{\prime}\\\ \end{array}\right),$ (5)

for AA stacking, for example. Here we use
$\omega^{\prime}=\omega_{AA^{\prime}}=\omega_{BB^{\prime}}=0.0939$ eV for the
intra-sublattice tunneling energy and
$\omega=\omega_{AB^{\prime}}=\omega_{BA^{\prime}}=0.12$ eV for the inter-
sublattice tunneling at zero pressure when we capture the out-of-plane lattice
relaxation effects parametrized through the polynomial

$\omega^{\prime}=C_{2}\omega^{2}+C_{1}\omega+C_{0}$ (6)

where $C_{2}=-0.5506$ eV-1, $C_{1}=1.036$, and $C_{0}=-0.02245$ eV as reported
in Ref. Chebrolu  _et al._ (2019). Otherwise, we can use equal
$\omega^{\prime}=\omega=0.12$ eV when we do not have the lattice corrugations
in the $z$-direction.

Figure 7:  (Color online) Bandwidth, secondary bandgap, and primary bandgap of
the lowest-energy valence (V) and conduction (C) bands of type I0 [see Fig. 4
(b) for the atomic configurations of the local commensurate stackings] for (a)
the minimal and (b) the remote hopping model as a function of the displacement
field and the twist angle. Here we set $\Delta^{\prime}=\Delta$ where
$\Delta~{}(\Delta^{\prime})$ represents the potential difference within
(between) bilayers. The solid lines in the first and the third columns
represent contours of the bandwidth of the lowest-energy valence and the
conduction bands. The solid lines in the second, fourth and fifth columns
represent the contour at the zero gaps. The cross-section at
$\theta=1.5^{\circ}$ is presented on the rightmost column with (purple) stars
for the V bandwidth, (red) triangles for the C bandwidth, (orange) squares for
the V secondary gap, (blue) diamonds for the C secondary gap, and (green)
circles for the primary gap.  Figure 8:  (Color online) The ratio of the bare
Coulomb interaction to the bandwidth $U/W$ and the ratio of the effective
Coulomb interaction to the bandwidth $U_{\textrm{eff}}/W$ [see Eq. (8)] for
the minimal and the remote hopping model as a function of the interlayer
potential difference $\Delta$ and twist angle $\theta$ of the lowest-energy
valence and conduction bands of type I(-,0,+) in (a), (b), (c), type II(-,0,+)
in (d), (e), (f), and type III(-,0,+). We set $\Delta^{\prime}=\Delta$ where
$\Delta~{}(\Delta^{\prime})$ represents the potential difference within
(between) bilayers.

In Fig. 1(e) we show schematic diagrams for four cases when the interlayer
coupling $\omega$ at twisted interfaces is considered as a perturbation. We
have three parabolic bands stemming from each bilayer as shown in the first
panel where there is no interlayer coupling at moire interfaces. If we have a
small interlayer coupling, for example, $\omega=0.01$ eV, the one parabolic
bands at $\tilde{K}^{\prime}$ and the one at $\tilde{K}$ start hybridizing as
denoted by $\Delta_{c}$ as seen in the second panel. $\Delta_{c}$ increases as
the twist angle grows such that $\Delta_{c}\approx 7.3$ meV, $7.9$ meV, $8.9$
meV at $\theta=1.0^{\circ}$, $1.5^{\circ}$, $2.5^{\circ}$, respectively for
$\omega=0.01$ eV. As shown in the third case, an external electric field
causes a band offset $2\Delta^{\prime}+\Delta$ between the two parabolic bands
at $\tilde{K}^{\prime}$ as well as the intrinsic gaps $\Delta$ in all three
parabolic bands without the interlayer coupling at moire interfaces. In the
fourth panel where we turn on interlayer coupling under the electric field
$\Delta$, the strength of hybridization $\Delta_{c}^{\prime}$ does not
increase monotonically with respect to the interlayer coupling unlike
$\Delta_{c}$ because of the another variable $\Delta$ in this case, and has
values $\approx 3.3$ meV, $5.0$ meV, $3.9$ meV at $\theta=1.0^{\circ}$,
$1.5^{\circ}$, $2.5^{\circ}$, respectively.

##  III Electronic properties

###  III.1 Lattice relaxation and electric field

The electronic band structures for the 9 possible cases as a starting stacking
due to (1) AB and BA alignments of the Bernal-stacked bilayer, and (2)
translations of the top bilayer with respect to the bottom bilayer by
$\bm{\tau}_{0}$, or $\bm{\tau}_{\pm}$ denoted as I(0,+,-), II(0,+,-), and
III(0,+,-) are quite distinct from each other. Here we first consider the type
I0 and investigate the effects of the lattice relaxation and the electric
field on the electronic band structures. See the Supplementary information for
the analog of Fig. 2 for the other types. We present in Fig. 2 the energy
dispersions and density of states (DOS) of the rigid lattice case
($\omega=\omega^{\prime}$) and the case of the out-of-plane relaxed lattice
($\omega\neq\omega^{\prime}$) for the minimal model in Figs. 2 (a) and 2 (b),
respectively, for the two twist angles, $\theta=2^{\circ}$ and $1.5^{\circ}$,
and for the two cases of interlayer potential differences,
$\Delta^{\prime}=\Delta=0.0$ eV and $\Delta^{\prime}=\Delta=0.02$ eV. In the
rigid lattice for zero electric field, at $\theta=2^{\circ}$, the energy
dispersions are similar in appearance such that the two parabolic bands
located at $\tilde{K}$ and $\tilde{K}^{\prime}$ are hybridized in common in
all the 9 cases. It is noteworthy that the other parabolic band dispersion at
$\tilde{K}^{\prime}$ occurs differently according to the starting stackings
and sliding types (not shown). It opens a large gap $|E|\gtrsim$ 0.05 eV at
$\tilde{K}^{\prime}$ in type I-, II-, and III0, while it remains parabolic as
seen in type I0 (Fig. 2), II0, and III+, or has a form of slightly linear
bands in types I+, II+, and III-. At $\theta=1.5^{\circ}$, the band
dispersions are quite different for each type, however, they have a common
feature that they have a gap at the $\tilde{\Gamma}$ point of the size of
$E_{gap,\tilde{\Gamma}}\gtrsim$ 0.05 eV regardless of the type. To be
specific, for type I0 as shown in Fig. 2, the two parabolic bands make nearly
flat bands with a gap-opening at $\tilde{\Gamma}$.

It has been reported that a corrugation effect induced by out-of-plane lattice
relaxation, which was observed in both experiments Cao  _et al._ (2016, 2018a,
2018b) and theoretical works Uchida  _et al._ (2014); van Wijk  _et al._
(2015) causes flat bands to be isolated from adjacent dispersing bands Koshino
_et al._ (2018). The changes in the electronic properties due to the lattice
relaxation are captured through unequal interlayer tunneling strengths at the
moire interface, $\omega$ and $\omega^{\prime}$ such that it follows the
relation of Eq. (6). The effects appear remarkably in type I0, II-, III-
lifting the degeneracies between the lowest-energy and the next higher-energy
bands whereas there are no qualitative changes in the electronic properties in
the other types.

On the other hand, the finite electric field $\Delta=0.02$ eV opens a primary
gap $\delta_{p}$ which is defined as a gap between the lowest-energy electron
and the hole bands as indicated by the red shaded region on the right panel of
lower row in Fig. 2 (b) for all the 9 cases for the relaxed lattice. The
secondary gap $\delta_{s}$ which is defined as a gap between the lowest-energy
electron (hole) band and the next higher-energy electron (hole) band as
indicated by the blue shaded regions on the right panel of lower row in Fig. 2
(b) barely opens only in type I0, I+, II+, III-, and III+ for the relaxed
lattice.

In the first row of Figs. 2(c) and 2(d), we show the bandwidths of the lowest-
energy electron and hole bands of the type I0 as a function of twist angle for
the rigid and relaxed lattice cases, respectively. Due to the gap-opening at
$\tilde{\Gamma}$, the bandwidth increases near $\theta=1.5^{\circ}$ as the
twist angle increases. The bandwidths evolve with twist angle similarly for
different interlayer coupling strengths $\omega=0.15$, $0.18$ eV but are
shifted approximately linearly with $\omega$Chittari  _et al._ (2018);
Chebrolu  _et al._ (2019). The different interlayer tunneling can be achieved
by applying pressure in the $z$-direction. We defer the detailed discussions
on the effect of pressure to the next subsection. Interestingly, the t3BG
system does not have a local minimum in bandwidth at the twist angle
$\theta=1.5^{\circ}$ where t3G has its first magic-angleKhalaf  _et al._
(2019); Carr  _et al._ (2020); Li  _et al._ (2019); Lopez-Bezanilla and Lado
(2020); Lei  _et al._ (2020); Hao  _et al._ (2021); Călugăru  _et al._ (2021);
Shin  _et al._ (2021), while t2BG has the local minimum in bandwidth at the
magic-angle $\theta=1.06^{\circ}$ which is also the magic-angle of t2G
Chittari  _et al._ (2018). Figs. 2 (e) and (f) give us the bandwidth as a
function of $\omega$. The bandwidth appears to be smaller than 20 meV at
$\theta=1.5^{\circ}$ in the range of $\omega\gtrsim$ 0.12 eV. This behavior
repeats at $\theta=1.0^{\circ}$ but in the range of $\omega\gtrsim$ 0.09 eV
and at $\theta=2.0^{\circ}$ in the range of $\omega\gtrsim$ 0.18 eV. The out-
of-plane lattice relaxation makes the graph of bandwidth smoother but the
overall physics does not change.

Similar to t3G where the electronic structure can be decomposed into a pair of
t2G at different angles with an additional decoupled monolayer graphene Khalaf
_et al._ (2019), we find that the electronic structure of t3BG can be
approximately described by a pair of t2BG at different angles with an
additional decoupled bilayer graphene at low energies near $\tilde{K}$ for I0
without an external field. Importantly, there are two parabolic bands at
$\tilde{K}^{\prime}$ which are exactly overlaid indicated by the dotted purple
lines, and one parabolic band at $\tilde{K}$ indicated by the dotted yellow
lines in the first row in Figs. 2 (a) and 2(b). One of the two parabolic bands
at $\tilde{K}^{\prime}$ corresponds to the decoupled bilayer graphene, and its
effective mass remains the same irrespective of the twist angle under the zero
electric field. For a finite electric field, the two parabolic bands at
$\tilde{K}^{\prime}$ start hybridizing with each other and the effective
masses of them change as indicated by the dotted purple, pink lines in the
second row in Figs. 2 (a) and 2(b).

###  III.2 Under pressure

In Fig. 3 we present an extended perspective of the electronic band structure
of I0 under pressure for the relaxed lattice case such as the bandwidth of the
lowest-energy electron and hole bands, the primary gap $\delta_{p}$, and the
secondary gaps $\delta_{s}$. See the Supplementary information for a more
complete survey. Based on the earlier studies Jung  _et al._ (2014); Leconte
_et al._ (2017); Yankowitz  _et al._ (2018); Chebrolu  _et al._ (2019), here
we introduce the same strategy used in Ref. Chebrolu  _et al._ (2019) given by
the two relations, between the external pressure and the interlayer distance,
and between the interlayer tunneling $\omega$ and the interlayer distance,
leading to a second-order polynomial as follows.

$P=A_{1}\omega^{2}+A_{2}\omega+A_{3},$ (7)

where we use the fitting parameters $A_{1}=324.7$ GPa/eV2, $A_{2}=-35.47$
GPa/eV and $A_{3}=-0.4671$ GPa obtained from EXX+RPA. The zero pressure (P = 0
GPa) is indicated by a black dashed line. The contour for the positive
(negative) level is drawn by the solid (dashed) lines. We juxtapose two
$y$-axes together – one for the pressure in the unit of GPa on the left side,
and the other for the interlayer tunneling strength $\omega$ on the right
side. We get electronic properties under a certain condition of external
pressure such that there appears a linear contour with an approximate relation
between the pressure and the twist angle such as $P\sim 3.3\theta$ for both
the minimal and remote hopping models for all the 9 cases as shown in Figs.
3(a) and 3(b). Type II- and III0 have finite secondary gaps even under the
zero electric field, while others have finite secondary gaps only when the
system is gated by an electric field. The result with the a finite interlayer
potential difference $\Delta^{\prime}=\Delta=0.02$ eV for type I0 is shown in
Figs. 3(c) and 3(d) as an example, and it has the positive secondary and
primary gaps in the range of $\theta\gtrsim 0.5^{\circ}$ along the linear
contour, while the primary and secondary gaps are closed in the entire
parameter space under the zero interlayer potential difference. Furthermore,
the finite interlayer potential difference tends to enlarge the size of the
primary gap in both minimal and remote cases. The cross-section at
$\theta=1.5^{\circ}$ is illustrated on the rightmost column. In the minimal
model without the displacement field the bandwidth (secondary gap) of the
conduction is symmetric to that of the valence band regardless of the
pressure. Upon turning on the remote hopping terms, the symmetry in bandwidth
(secondary gap) between the conduction band and the valence band starts
breaking as seen in the rightmost column in Figs. 3 (b) and 3(d). In contrast,
the displacement field preserves the analogous features between the conduction
and valence bands but turns the monotonically decreasing (increasing)
bandwidth (secondary gap) into a non-monotonic function as the pressure
increases as shown in the rightmost column in Fig. 3 (c).

Figure 9:  (Color online) Topological phase diagrams for the $K$-valley Chern
numbers of the lowest-energy valence and conduction bands of t3BG of type
I(-,0,+) in (a), (b), (c), type II(-,0,+) in (d), (e), (f), and type
III(-,0,+) in (g), (h), (i) in the minimal (upper row) and the remote (lower
row) model as a function of the interlayer potential difference $\Delta$ and
the twist angle. Here we use the condition $\Delta^{\prime}=\Delta$ where
$\Delta~{}(\Delta^{\prime})$ represents the gated potential difference within
(between) bilayers.

###  III.3 Local density of states

Each of the 9 stacking geometries shown in Figs. 1 (a)-(c) gives moire
patterns consisting of three different local commensurate stackings, a, b, and
c as shown in Fig. 1 (d) with the three local commensurate stackings in the
moire patterns for I0 as an example. In principle, one can get the same band
structure no matter which stackings among these three local commensurate
stackings we take as the starting stacking. We listed the atomic
configurations of the local commensurate stackings for type I(-,0,+) in the
first row in Fig. 4, for type II(-,0,+) in Fig. 5, and for type III(-,0,+) in
Fig. 6. Throughout Figs. 4-6, the twist angle is set to $\theta=1.5^{\circ}$,
the energy dispersions and DOS for the minimal model are indicated with red
dashed lines, and the ones with the remote hopping terms are in black solid
lines for the energy band and blue-filled areas for the DOS.

The electronic features are quite distinct for the 9 different geometries I,
II, III(-,0,+) as shown in Figs. 4-6. It is common across all types of
stackings that the minimal model has particle-hole symmetric band, and the
introduction of the remote hopping terms breaks the electron-hole symmetry
especially of the lowest-energy bands resulting in DOS broadening. We show the
band structures under a finite displacement field
$\Delta^{\prime}=\Delta=0.02$ eV and it is noteworthy that all the 9 cases
have a finite primary gap $\delta_{p}$ and some of cases have finite secondary
gaps $\delta_{s}$ as well, leading to the isolated lowest-energy bands.

We show the local density of states (LDOS) at the van Hove singularities
denoted by the red arrow for the minimal and the blue arrow for the remote
hopping terms model. We found in common for all the 9 cases that the local-AA
stacking regions at the moire interface have relatively high LDOS
concentrations. For instance, among the local commensurate stackings of type
I-, the stacking c has the local-AA configurations at both moire interfaces,
yielding the highly concentrated LDOS, unlike the stackings a and b. In type
I0, on the other hand, the local commensurate stackings b and c have the
local-AA configurations at one of the two moire interfaces but a does not have
the local-AA stackings at both interfaces. Thus, the intensity of LDOS at the
stackings b and c surpasses that of the stacking a for a moire length of
$l_{M}=9.4$ nm.

###  III.4 Bandwidths and bandgaps

In the following, we explore the evolution of the bandwidths, secondary, and
primary bandgaps in the parameter space of the displacement field, and the
twist angle for I0 taking Fig. 7 as an example. See the supplementary
information for other stackings. All 9 starting stacking cases have in common
that the valence and conduction bandwidths are not suppressed around
$\theta=1.5^{\circ}$ but the secondary gaps $\delta_{s}$ become open at
several large or small islands around this angle for finite displacement
fields. Besides, the primary gap $\delta_{p}$ is open in a large range of the
parameter space, a fact that will facilitate the band isolation. The rightmost
column represents cross-sections at $\theta=1.5^{\circ}$ for the bandwidths
and bandgaps. In these plots the difference between both minimal and remote
hopping models becomes more apparent that the electron-hole symmetry of the
bands in the minimal model is preserved except for type III+ where the
electron-hole asymmetry with respect to electric field is already broken in
the minimal model. More specifically, the type I0 in Fig. 7 has finite
secondary gaps at $\theta\lesssim 0.2^{\circ}$ for $\Delta>0$, at several
small islands $0.2^{\circ}\lesssim\theta\lesssim 0.7^{\circ}$ in the range of
$\Delta\lesssim 0.04$ eV, and at larger islands
$1^{\circ}\lesssim\theta\lesssim 2^{\circ}$ in the range of $\Delta\lesssim
0.06$ eV.

Figure 10:  (Color online) Band structures at $\theta=1.5^{\circ}$ and the
corresponding Berry curvatures for the lowest-energy valence and conduction
bands with the $K$-valley Chern number [see Fig. 9] on top of them for type
I(-,0,+) in (a), (b), (c), type II(-,0,+) in (d), (e), (f), and type
III(-,0,+) in (g), (h), (i). Here we set $\Delta^{\prime}=\Delta=0.02$ eV for
both the minimal model and the remote hopping model where
$\Delta~{}(\Delta^{\prime})$ represents the potential differences within
(between) bilayers.

##  IV Effective Coulomb interactions

We evaluate the relative strength of the Coulomb interactions with partial
screening due to the metallicity of the bands as in Ref. Chebrolu  _et al._
(2019):

$U_{\textrm{eff}}=U\exp{(-l_{M}/\lambda_{D})},$ (8)

where $U$ represents the bare Coulomb interaction which is defined as
$\frac{e^{2}}{4\pi\epsilon_{r}\epsilon_{0}l_{M}}$, the moire length $l_{M}\sim
a/\theta$ when the twist angle $\theta$ is small enough, and $\lambda_{D}$ is
the Debye length which can be expressed as
$\lambda_{D}=2\epsilon_{0}/[e^{2}D(\delta_{p},\delta_{s})]$ with the two-
dimensional DOS $D(\delta_{p},\delta_{s})$. Here, we define
$D(\delta_{p},\delta_{s})$ as

$D(\delta_{p},\delta_{s})=4\frac{|\delta_{p}|u(-\delta_{p})+|\delta_{s}|u(-\delta_{s})}{W^{2}A_{M}},$
(9)

where $u(x)$ is the Heaviside step function so that $D(\delta_{p},\delta_{s})$
is proportional to the ratio of overlapping bands to the bandwidth
$|\delta_{p(s)}|/W$ if $\delta_{p}<0$ ($\delta_{s}<0$), $A_{M}$ is the area of
a moire unit cell in real space given as $\sqrt{3}~{}l^{2}_{M}/2$. In the
current work, we set the dielectric constant of graphene as $\epsilon_{r}=4$
Jung and MacDonald (2013).

We discuss the ratio of the bare $U$ and the effective Coulomb interaction
$U_{\textrm{eff}}$ to the bandwidth $W$ as a function of the displacement
field within bilayers $\Delta$ and the twist angle for both minimal and remote
hopping models for the all 9 cases I(-,0,+), II(-,0,+), and III(-,0,+) shown
in Fig. 8. For all cases, the bare Coulomb interaction over the bandwidth has
the value of $U/W\gtrsim 1$ in almost the entire parameter space and has the
maximum values at around $0.2^{\circ}\lesssim\theta\lesssim 1.5^{\circ}$ for
weak displacement fields $|\Delta|\lesssim 25$ meV. Apart from type I-, II-,
III0, the strongly correlated electron states extended like a long tail above
$|\Delta|\sim 25$ meV for both the minimal and remote models. On the other
hand, if we take into account the screening effect due to the conducting
electrons, namely the effective Coulomb interactions $U_{\textrm{eff}}$, the
territory where the strongly correlated states appear is suppressed for both
the minimal and remote hopping cases and remain finite only at small islands
in the range of $1^{\circ}\lesssim\theta\lesssim 2^{\circ}$ when
$|\Delta|\lesssim 0.05$ eV in all the cases. To be more specific, in type I0,
I+, and III- relatively larger islands survive under the inclusion of the
screening effects by the itinerant electrons in the valence band of the remote
hopping model at $\theta=1.5^{\circ}$ under a week displacement field
$|\Delta|\lesssim 0.05$ eV. On the other hand, the strong effective Coulomb
regions are barely found at around $\theta=1.5^{\circ}$ in type I-, I+, II0,
II+, and III+ since it is cut off by the electrons in the first and second
lowest-energy bands overlapping with each other. It is also noteworthy in Fig.
8 that both $U/W$ and $U_{\textrm{eff}}/W$ are symmetric against $\Delta=0$ eV
line except the three types III(0,-,+) where the stacking configurations of
the bilayers on top (BA) and bottom (AB) are different to each other.

##  V Valley Chern numbers

We now turn our attention to the topological nature of the bands in t3BG by
evaluating the Chern number of the lowest-energy valence and conduction bands
only for one valley. Here we calculate the Chern number of $n^{\rm th}$ energy
band $C_{n}$ for the $+K$ valley which is defined as

$C_{n}=\frac{1}{2\pi}\int_{\textrm{mBZ}}d^{2}\bm{k}~{}\Omega_{n}(\bm{k}),$
(10)

by integrating the Berry curvature $\Omega_{n}(\bm{k})$ within the moire
Brillouin zone (mBZ) which is given as Xiao  _et al._ (2010)

$\Omega_{n}(\bm{k})=-2\sum_{n^{\prime}\neq n}\textrm{Im}\Bigg{[}\frac{\langle
n|\frac{\partial H}{\partial k_{x}}|n^{\prime}\rangle\langle
n^{\prime}|\frac{\partial H}{\partial
k_{y}}|n\rangle}{(E_{n^{\prime}}-E_{n})^{2}}\Bigg{]}.$ (11)

Here, $E_{n}$ represents an energy eigenvalue of the moire superlattice Bloch
state $|n\rangle$.

In Fig. 9 we show the $K$-valley Chern numbers of the lowest-energy valence
and conduction bands for all 9 starting stackings without/with remote hopping
terms. In all cases, the valence and conduction bands have well-defined Chern
numbers in almost the entire parameter space of the displacement field within
bilayers $|\Delta|\leq 0.1$ eV and the twist angle $1^{\circ}\leq\theta\leq
2.5^{\circ}$. The inclusion of the remote hopping terms barely changes the
phase diagram but the borders between topologically distinct phases are
slightly shifted. Apart from the three types III(-,0,+) the $K$-valley Chern
numbers of the valence and conduction bands tend to be the opposite to each
other such that the sum of them is zero.

For the I0, I+, and III- cases where the regions of strongly correlated states
survive in a slightly larger area under the screening effects in the remote
hopping model as discussed in the previous section, type I0 has the
topologically trivial valence and conduction bands near $\theta=1.5^{\circ}$,
while type I+ has a topologically nontrivial phase with C$=+1~{}(-1)$ in the
valence (conduction) band for $\Delta>0$ [C$=-1~{}(+1)$ in the valence
(conduction) band for $\Delta<0$], and type III- with C$=+1~{}(0)$ in the
valence (conduction) band for $\Delta>0$ [C$=0~{}(+1)$ in the valence
(conduction) band for $\Delta<0$] in the areas of the correlated electrons.

In Fig. 10 we present the band structures and the corresponding Berry
curvatures $\Omega_{n}(\bm{k})$ for the lowest-energy electron and hole bands
at $\theta=1.5^{\circ}$ for all 9 cases I(-,0,+), II(-,0,+), and III(-,0,+)
under a certain positive finite value of interlayer difference field,
$\Delta^{\prime}=\Delta=0.02$ eV for both the minimal and remote cases. We
highlight the valence (conduction) band with the red (blue) line in the band
structure figure and indicate the individual band topological invariant on top
of each Berry curvature figure. The mBZ is indicated by a black solid hexagon.
Thus one can see the effect of the remote hopping terms on the Berry
curvatures. We see that the topological properties generally remain intact
upon inclusion of the remote hopping terms except for the valence band of type
II- and the conduction band of type II0.

##  VI SUMMARY

We have presented the electronic structure analysis of alternatingly twisted
triple Bernal-stacked bilayer graphene in search of new multilayer graphene
devices that give rise to nearly flat lowest-energy Chern bands. Twisted
systems made up of bilayer graphene are interesting because each layer can
open a bandgap upon application of an electric field. Since they have distinct
AB or BA stacking orientations we have considered 9 starting stackings by
considering three different sliding vectors $\bm{\tau}_{0,\pm}$ for the top
bilayer with respect to the bottom bilayer.

More specifically, we have first described the detailed electronic properties,
such as bandwidths and bandgaps of the lowest-energy bands in a wide parameter
space of the interlayer couplings strength $\omega$ varied by the pressure,
the twist angle $\theta$, and the interlayer potential difference $\Delta$ in
the perpendicular direction that can be introduced by an electric field. We
have found that in all 9 starting stackings, the lowest-energy bandwidths can
be easily suppressed to values of the order of $\sim$20 meV near
$\theta\approx 1.5^{\circ}$ and even narrower bands can be obtained for
smaller twist angles. We have compared the band structures for the minimal
model with the remote hopping model, and we have shown the DOS and LDOS map
giving rise to relatively high densities at the local-AA stacking sites of the
moire interfaces.

Our numerical analysis shows the electric field and twist angles where strong
electron-electron interactions can be expected from the ratio between the
bandwidth and the screened Coulomb potentials. We have found that type I0, I+,
and III- are more favorable than other stackings to find correlated phases in
the vicinity of $\theta=1.5^{\circ}$ under a finite interlayer potential
difference of $\Delta\sim 20$ meV.

The valley Chern numbers phase diagrams for the lowest-energy electron and
hole bands for the minimal and remote hopping models have been obtained for
all 9 stacking configurations for a range of twist angles below 2.5∘ and
$\left|\Delta\right|<0.1$ eV and the associated Berry curvatures in the
momentum space at the twist angle of $\theta=1.5^{\circ}$ and
$\Delta~{}(\Delta^{\prime})\approx 20$ meV. The well-defined valley Chern
numbers in the minimal model are mostly robust upon inclusion of the remote
hopping terms except for two cases, the valence band in type II- and the
conduction band in type II0, where the remote hopping terms lead to band-
touching with higher-energy bands. We expect that our work will provide
guidance in the search of nearly flat bands in twisted triple bilayer devices.

###### Acknowledgements.

This work was supported by Samsung Science and Technology Foundation Grant No.
SSTF-BA1802-06 (J.S.), Korean NRF through the Grants No. 2021R1A6A3A01087281
(J.S.), No. 2020R1A5A1016518 (B.L.C.), No. 2018R1A2B6007837 (Y.J. and H.M.),
No. 2020R1A2C3009142 (J.J.), and Creative-Pioneering Researchers Program
through Seoul National University (SNU) (Y.J. and H.M.). We acknowledge
computational support from KISTI Grant No. KSC-2021-CRE-0389 and by the
resources of Urban Big data and AI Institute (UBAI) at UOS. J.J also
acknowledges support by the Korean Ministry of Land, Infrastructure and
Transport (MOLIT) from the Innovative Talent Education Program for Smart
Cities and the KREONET network infrastructure.

## References

  * Kim  _et al._ (2017) K. Kim, A. DaSilva, S. Huang, B. Fallahazad, S. Larentis, T. Taniguchi, K. Watanabe, B. J. LeRoy, A. H. MacDonald, and E. Tutuc, PNAS 114, 3364 (2017).
  * Cao  _et al._ (2018a) 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, Nature 556, 80 (2018a).
  * Cao  _et al._ (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Yankowitz  _et al._ (2019) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Cao  _et al._ (2020a) Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020a).
  * Park  _et al._ (2021a) J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 590, 249 (2021a).
  * Hao  _et al._ (2021) Z. Hao, A. M. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Science 371, 1133 (2021).
  * Lee  _et al._ (2019) J. Y. Lee, E. Khalaf, S. Liu, X. Liu, Z. Hao, P. Kim, and A. Vishwanath, Nat. Commun. 10, 5333 (2019).
  * Chebrolu  _et al._ (2019) N. R. Chebrolu, B. L. Chittari, and J. Jung, Phys. Rev. B 99, 235417 (2019).
  * Choi and Choi (2019) Y. W. Choi and H. J. Choi, Phys. Rev. B 100, 201402(R) (2019).
  * Koshino (2019) M. Koshino, Phys. Rev. B 99, 235406 (2019).
  * Burg  _et al._ (2019) G. W. Burg, J. Zhu, T. Taniguchi, K. Watanabe, A. H. MacDonald, and E. Tutuc, Phys. Rev. Lett. 123, 197702 (2019).
  * Cao  _et al._ (2020b) Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 583, 215 (2020b).
  * Shen  _et al._ (2020) C. Shen, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, Y. Chu, K. Watanabe, T. Taniguchi, R. Yang, Z. Y. Meng, D. Shi, and G. Zhang, Nat. Phys. 16, 520 (2020).
  * Liu  _et al._ (2020) X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, Y. Ronen, H. Yoo, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Nature 583, 221 (2020).
  * Chittari  _et al._ (2019) B. L. Chittari, G. Chen, Y. Zhang, F. Wang, and J. Jung, Phys. Rev. Lett. 122, 016401 (2019).
  * Suárez Morell  _et al._ (2013) E. Suárez Morell, M. Pacheco, L. Chico, and L. Brey, Phys. Rev. B 87, 125414 (2013).
  * Li  _et al._ (2019) X. Li, F. Wu, and A. H. MacDonald,  (2019), arXiv:1907.12338 .
  * Szendrő  _et al._ (2020) M. Szendrő, P. Süle, G. Dobrik, and L. Tapasztó, (2020), arXiv:2001.11462 .
  * Carr  _et al._ (2020) S. Carr, C. Li, Z. Zhu, E. Kaxiras, S. Sachdev, and A. Kruchkov, Nano Lett. 20, 3030 (2020).
  * Ma  _et al._ (2021) Z. Ma, S. Li, Y.-W. Zheng, M.-M. Xiao, H. Jiang, J.-H. Gao, and X. Xie, Science Bulletin 66, 18 (2021).
  * Xu  _et al._ (2021) S. Xu, M. M. A. Ezzi, N. Balakrishnan, A. Garcia-Ruiz, B. Tsim, C. Mullan, J. Barrier, N. Xin, B. A. Piot, T. Taniguchi, K. Watanabe, A. Carvalho, A. Mishchenko, A. K. Geim, V. I. Fal’ko, S. Adam, A. H. C. Neto, K. S. Novoselov, and Y. Shi, Nat. Phys. 17, 619 (2021).
  * Chen  _et al._ (2021) S. Chen, M. He, Y.-H. Zhang, V. Hsieh, Z. Fei, K. Watanabe, T. Taniguchi, D. H. Cobden, X. Xu, C. R. Dean, and M. Yankowitz, Nat. Phys. 17, 374 (2021).
  * 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, Nature 588, 66 (2020).
  * Park  _et al._ (2020) Y. Park, B. L. Chittari, and J. Jung, Phys. Rev. B 102, 035411 (2020).
  * Lei  _et al._ (2020) C. Lei, L. Linhart, W. Qin, F. Libisch, and A. H. MacDonald, (2020), arXiv:2010.05787 .
  * Wu  _et al._ (2021) Z. Wu, Z. Zhan, and S. Yuan, Sci. China Phys. Mech. Astron. 64, 267811 (2021).
  * Mora  _et al._ (2019) C. Mora, N. Regnault, and B. A. Bernevig, Phys. Rev. Lett. 123, 026402 (2019).
  * Zuo  _et al._ (2018) W.-J. Zuo, J.-B. Qiao, D.-L. Ma, L.-J. Yin, G. Sun, J.-Y. Zhang, L.-Y. Guan, and L. He, Phys. Rev. B 97, 035440 (2018).
  * Khalaf  _et al._ (2019) E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Phys. Rev. B 100, 085109 (2019).
  * Lopez-Bezanilla and Lado (2020) A. Lopez-Bezanilla and J. L. Lado, Phys. Rev. Research 2, 033357 (2020).
  * Călugăru  _et al._ (2021) D. Călugăru, F. Xie, Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Phys. Rev. B 103, 195411 (2021).
  * Shin  _et al._ (2021) J. Shin, B. L. Chittari, and J. Jung, Phys. Rev. B 104, 045413 (2021).
  * Min  _et al._ (2007) H. Min, B. Sahu, S. K. Banerjee, and A. H. MacDonald, Physical Review B 75, 155115 (2007).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, PNAS 108, 12233 (2011).
  * Jung  _et al._ (2014) J. Jung, A. Raoux, Z. Qiao, and A. H. MacDonald, Phys. Rev. B 89, 205414 (2014).
  * Jung and MacDonald (2014) J. Jung and A. H. MacDonald, Phys. Rev. B 89, 035405 (2014).
  * Park  _et al._ (2021b) J. M. Park, Y. Cao, L. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero,  (2021b), 10.48550/ARXIV.2112.10760.
  * Cao  _et al._ (2016) Y. Cao, J. Y. Luo, V. Fatemi, S. Fang, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Phys. Rev. Lett. 117, 116804 (2016).
  * Uchida  _et al._ (2014) K. Uchida, S. Furuya, J.-I. Iwata, and A. Oshiyama, Phys. Rev. B 90, 155451 (2014).
  * van Wijk  _et al._ (2015) M. M. van Wijk, A. Schuring, M. I. Katsnelson, and A. Fasolino, 2D Mater. 2, 034010 (2015).
  * Koshino  _et al._ (2018) M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys. Rev. X 8, 031087 (2018).
  * Chittari  _et al._ (2018) B. L. Chittari, N. Leconte, S. Javvaji, and J. Jung, Electron. Struct. 1, 015001 (2018).
  * Leconte  _et al._ (2017) N. Leconte, J. Jung, S. Lebègue, and T. Gould, Phys. Rev. B 96, 195431 (2017).
  * Yankowitz  _et al._ (2018) M. Yankowitz, J. Jung, E. Laksono, N. Leconte, B. L. Chittari, K. Watanabe, T. Taniguchi, S. Adam, D. Graf, and C. R. Dean, Nature 557, 404 (2018).
  * Jung and MacDonald (2013) J. Jung and A. H. MacDonald, Phys. Rev. B 88, 075408 (2013).
  * Xiao  _et al._ (2010) D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).