# Extended magic phase in twisted graphene multilayers

D. C. W. Foo<EMAIL_ADDRESS>Centre for Advanced 2D Materials, National
University of Singapore, 6 Science Drive 2, Singapore 117546  Z. Zhan  IMDEA
Nanociencia, C/ Faraday 9, 28049 Madrid, Spain  Mohammed M. Al Ezzi  Centre
for Advanced 2D Materials, National University of Singapore, 6 Science Drive
2, Singapore 117546  Department of Physics, Faculty of Science, National
University of Singapore, 2 Science Drive 3, Singapore 117542  L. Peng
Department of Physics, Faculty of Science, National University of Singapore, 2
Science Drive 3, Singapore 117542  S. Adam  Centre for Advanced 2D Materials,
National University of Singapore, 6 Science Drive 2, Singapore 117546
Department of Physics, Faculty of Science, National University of Singapore, 2
Science Drive 3, Singapore 117542  Department of Materials Science and
Engineering, National University of Singapore, 9 Engineering Drive 1,
Singapore 117575  Yale-NUS College, 16 College Ave West, Singapore 138527  F.
Guinea  IMDEA Nanociencia, C/ Faraday 9, 28049 Madrid, Spain  Donostia
International Physics Center, Paseo Manuel de Lardizábal 4, 20018, San
Sebastián, Spain

###### Abstract

Theoretical and experimental studies have verified the existence of “magic
angles” in twisted bilayer graphene, where the twist between layers gives rise
to flat bands and consequently highly correlated phases. Narrow bands can also
exist in multilayers with alternating twist angles, and recent theoretical
work suggests that they can also be found in trilayers with twist angles
between neighboring layers in the same direction. We show here that flat bands
exist in a variety of multilayers where the ratio between twist angles is
close to coprime integers. We generalize previous analyses, and, using the
chiral limit for interlayer coupling, give examples of many combinations of
twist angles in stacks made up of three and four layers which lead to flat
bands. The technique we use can be extended to systems with many layers. Our
results suggest that flat bands can exist in graphene multilayers with angle
disorder, that is, narrow samples of turbostatic graphite.

Introduction  – Twisted bilayer graphene (tBG) has proven to be a fertile
platform for condensed matter physics studies owing to the remarkable
tunability of the Fermi velocity and electronic structure through the twist
angle [1]. In particular, at the so-called “magic angles” [2] for tBG,
$\theta^{\mathrm{M}}_{\mathrm{tBG}}\approx 1.1^{\circ}$, the Fermi velocity
vanishes, allowing for strongly correlated phases where interactions dominate
over the kinetic energy [3, 4].

This has led to a resurgence of interest in moiré systems, including but not
limited to Van der Waals heterostructures [5, 6, 7, 8], transition metal
dichalcogenides [9, 10, 11, 12], and of course graphene multilayers [13, 14,
15, 16]. One particularly successful avenue of investigation has been the
extension of the arguments demonstrating the existence of magic angles in tBG
[17, 18] to graphene multilayers [19], which found simple geometric relations
between the bilayer and multilayer magic angles, predictions that were later
verified experimentally [20, 21, 22]. Later numerical studies have built upon
these foundations, demonstrating for example that the existence of flat bands
is robust against perpendicular magnetic fields [23] and asymmetry of the
layer parameters [24], and that peaks in the density of states (DOS) appear,
and are resistant to relaxation, for trilayer configurations beyond the
canonical alternating-twist case [25]. Recent work has predicted magic angles
in the equal twist case [26], and a magic line, with magic angle a function of
pressure, has been predicted in twisted WSe2 [27]. In this work we demonstrate
the existence of flat bands in multilayers with generic, commensurate twist
angles, indicating an extended $N-2$ dimensional magic surfaces in the space
of $N-1$ twist angles characterising $N$ layers, with particular emphasis on
generic twisted trilayer graphene (tTG).

A prior study [19] noted that in the particular case of unshifted, alternating
twist graphene multilayers, the Hamiltonian may be exactly decomposed into a
direct sum of bilayer graphene Hamiltonians (with an additional monolayer for
odd number of layers), with the moiré potential identical up to a scalar
multiple. This allows one to directly read off the magic angles of multilayer
stacks. In contrast, other trilayer configurations have more complex effective
moiré potentials, though magic angle physics may still be possible [28]. In
particular, a strongly correlated superconducting state has been observed in a
trilayer system with incommensurate twist angles [29].

Moiré lattices in twisted graphene trilayers and in other twisted stacks. –
The electronic structure of a twisted bilayer can be reasonably approximated
by a periodic moiré structure [1, 2], irrespective of whether the two layers
are commensurate at the atomic scale. The situation is more complex in twisted
stacks with more than two layers [30, 31, 32, 33, 34, 35, 36, 25]. In a
trilayer, for instance, two twist angles between nearest neighbor layers can
be defined, $\theta_{12}$ and $\theta_{23}$, shown in Fig. 1(a). Except for
the case mentioned above of multilayers where the twist angle is of equal
magnitude and opposite sign between nearest neighbor pairs of layers [19], no
simple global moiré structure can be defined.

An approximate moiré lattice can be defined if the relative angle between the
Brillouin Zones defined by the twist angles between different pairs of layers
is neglected [32, 34, 35]. For instance, in a trilayer with two different
twist angles, $\theta_{12},\theta_{23}$, this approximation becomes exact when
$\theta_{12}=\theta_{23}$ in the limit $\theta_{12}\rightarrow 0$, although
corrections [36] are required for physically relevant angles,
$\theta_{12},\theta_{23}\approx 1^{\circ}$. It is worth noting that this
approximation, for a trilayer, has a dependence not only on the angles
$\theta_{12},\theta_{23}$, but also on the relative displacement of the moiré
potentials, $\mathbf{D}$ illustrated in Fig. 1(b) [33]. The dependence is
periodic on $\mathbf{D}$, with wavelength of the order of the moiré length
scales.

Figure 1: Schematic illustrating the trilayer parameters (a) twist angle pairs
and (b) the moiré displacement $\mathbf{D}$. The moiré unit cell of tBG with
$\theta_{12}$ is outlined by a rhombus with alternating colors red and blue,
and $\theta_{23}$ with colors red and green. The rhombus vertices lie at $AA$
stacking regions of the respective layers. (c) Sketch of connected momentum
states between layers with misaligned moiré patterns (orange and purple),
resulting in a moiré-of-moiré modulation. Large angles
$\theta_{12}=4^{\circ}$, $\theta_{23}=3^{\circ}$ are used for clarity, and
less misalignment is expected for real magic angle samples. The lengthscale of
this modulation is determined by the closest approach of second layer lattice
sites formed from the two moiré patterns, indicated with a dotted line. The
local expansion used in numerical calculation neglects this misalignment and,
for commensurate angles, collapses this distance to zero. (d) Log-scale plot
of the moiré-of-moiré length scale, $l_{m2}$, in microns, as a function of
twist angles. The left diagonal in white marks the line
$\theta_{12}=-\theta_{23}$, where there is no misalignment and $l_{m2}$
diverges. $l_{m2}$ also diverges when either angle is zero and there is no
moiré-of-moiré pattern to speak of.

It can be shown [37] that the approximation mentioned above defined by a moiré
lattice commensurate with the two twist angles and a displacement between the
top and bottom layers describes locally the full problem of the trilayer. For
the case of approximately commensurate twist angles, $\theta_{12}\approx
m\theta_{0},\theta_{23}\approx n\theta_{0}$, where $m,n$ are coprime integers,
the interlayer tunneling terms in the trilayer continuum Hamiltonian can be
written as containing two periodicities: i) the periodicity, $\ell_{m}$, of
the moiré problem obtained by neglecting the angle between the two Brillouin
zones defined by $\theta_{12},\theta_{23}$, of order $d/\theta_{0}$, where $d$
is the lattice constant of graphene, and ii) a periodicity, $\ell_{m2}$
obtained by modulating the displacement $\mathbf{D}$, over a larger unit cell
whose dimensions scale as $d/\theta_{0}^{2}$. Each tBG pair in tTG defines a
momentum-space moiré unit cell, purple and orange in Fig. 1(c), whose
misalignment leads to misalignment of the moiré lattices, shown with the black
lines and the finite $\delta k$. Neglecting this moiré misalignment forces the
black lines to coincide, and for commensurate angles forces $\delta k=0$.
$\delta k$ therefore is conjugate to $\ell_{m2}$. For small values of
$\theta_{0}$, shown in Fig. 1(d), this scaling behavior implies that
$\ell_{m}\ll\ell_{m2}$, as $\ell_{m2}$ can reach values of microns for
$\theta_{0}\lesssim 1^{\circ}$. Further details may be found in [38].

It is worth noting that small and layer-dependent biaxial strains can exactly
align the Brillouin zones associated to different twist angles, $\theta_{12}$,
$\theta_{23}$, so that $\ell_{m2}\rightarrow\infty$ and the approximation
which makes use of a single moiré lattice defined by $\ell_{m}$ becomes exact,
see [39]. These strains imply a small relaxation of the lattice [39]. This
analysis has been extended to some commensurate cases of trilayers with
unequal twist angles [40].

Results. Flat bands in the chiral limit. – We first analyze numerically the
existence of flat bands in chiral model trilayers [17, 18] for different
combinations of twist angles, $\theta_{12},\theta_{23}$ and relative
displacements between the top and bottom layers. The chiral limit for graphene
trilayers with the same twist angle between nearest layers,
$\theta_{12}=\theta_{23}$ has been recently explored in [26, 41, 39]. The
parameters that we use [42] in the following are $\gamma_{0}=2.6$ eV
(intralayer nearest neighbor hopping), and $t_{AB}=0.0975$ eV (interlayer
hopping). The choice of $\gamma_{0}$ leads to a Fermi velocity $v_{F}\approx
840$ km/s.

The analytical arguments for the existence of infinitely flat central bands in
twisted bilayer graphene [18] is based on the existence of zeros at the
corners of the unit cell of a 2 component spinor which describes the periodic
part of the wavefunction at the $K$ point of the Brillouin Zone. As one of the
two components, by symmetry, vanishes at the relevant points, the problem
reduces to the tuning of the zeros of a single component. The argument can be
generalized to a trilayer, where a three component spinor is involved. When
the Hamiltonian has ${\cal C}_{3}$ symmetry [18] two of the three components
of the spinor vanish at the relevant point in the real space unit cell, so
that the existence of infinitely flat bands reduces to the existence of a zero
in a single component of the spinor. The ${\cal C}_{3}$ symmetry exists at
least for a displacement $\mathbf{D}=0$, and for the displacements which take
the trilayer from the $AAA$ to the $ABA$ configuration. An alternative
analysis, for the case of a trilayer with $\theta_{12}=\theta_{23}$, can be
found in [26]. We numerically find flat bands, whose width is continuously
lowered as the angle is determined with better precision, for all combinations
of integers, $m,n$ and displacements, $\mathbf{D}$, tried. Results are shown
in Fig. 2. Interestingly, in some cases we find a quadruplet of bands nearly
degenerate at zero energy. We show in Fig. 2 the local DOS at the magic angles
corresponding to $\\{m,n\\}=\\{1,1\\}$ and $\\{m,n\\}=\\{3,-4\\}$ obtained for
a displacement between top and bottom layers $\mathbf{D}=0$. For each magic
angle, we plot the DOS for the selected displacement, $\mathbf{D}=0$, and for
two other displacements, shown in the inset in Fig. 2 (a). There is a
significant dependence of the local DOS on $\mathbf{D}$ in the $\\{1,1\\}$
case, and a less pronounced dependence in the $\\{3,-4\\}$ case. This trend
towards a weaker dependence on $\mathbf{D}$ with increasing values of the
integers $\\{m,n\\}$ is observed for all combinations explored, see [38].

Figure 2: (a) Local DOS of a graphene trilayer in the chiral limit with
$\theta_{12}=\theta_{23}=\theta_{0}=0.72^{\circ}$ (magic angle for zero
displacement) for three different displacements $\mathbf{D}$ between the top
and bottom layers. The insert on the top panel are the three displacements,
which are marked by a dot in the hexagon with the color corresponding to that
of the local DOS. (b) As in (a), but for a combination of angles
$\theta_{12}=3\theta_{0}$, $\theta_{23}=-4\theta_{0}$, and
$\theta_{0}=0.404^{\circ}$. Note the difference in the scales between the
plots in (a) and (b). (c) Calculated magic angle pairs ($\theta_{12}$,
$\theta_{23}$) for twisted trilayer graphene. The red and green colors are for
magic angles with red and green displacements illustrated in insert of (a),
respectively. The magic angle pair of (a) is marked by a blue star and that of
(b) by a blue triangle. The blue square is the magic angle in twisted trilayer
graphene with mirror symmetry. (d) Magic angle triplets ($\theta_{12}$,
$\theta_{23}$, $\theta_{34}$) for stacks of four layers. Different colors are
used for clarity.

Alternatively, the $3\times 3$ matrix equation which describes the spinor at
the magic angles of the chiral bilayer [18, 26] can be approximately reduced
to an effective bilayer, making significantly easier the calculation of magic
angles as function of commensuration and displacement, see [38].

Figure 3: Local DOS for a trilayer in the chiral limit with angles
$\theta_{12}=3\theta_{0}$, $\theta_{23}=-4\theta_{0}$,
$\theta_{0}=0.47^{\circ}$ and different displacements between the top and
bottom layers. The displacements are described by the phases $\\{0,0,0\\}$
(a), $\\{0,-\pi/18,\pi/18\\}$ (b), and $\\{\pi/36,-\pi/18,\pi/36\\}$ (c). (d)
An average of the three densities of states with the weights $\\{1:2:3\\}$, is
shown in (d).

The scheme discussed here for a trilayer [32, 34, 35, 37], based on two twist
angles, two integers, and one interlayer displacement can be extended to an
$N$ layer model, with $N-1$ twist angles, and $N-2$ independent displacements.
For the case of $N-$layer stacks where the twist angles between neighboring
layers are equal, an analytical expression for the first magic angle, which
generalizes the results in [26], gives:

$\displaystyle\theta_{N}^{M}$
$\displaystyle=\theta_{2}^{M}\left[\frac{1}{(N-1)!}\right]^{\frac{1}{N-1}}$
(1)

This expression gives, for $N\rightarrow\infty$
$\lim_{N\rightarrow\infty}\theta_{N}^{M}\approx\theta_{2}^{M}\times[]e(N-1)^{-1}]$
where $e$ is Euler’s number. Figure 2(d) shows results for magic angle
combinations in stacks of four layers, defined by three twist angles, and
three integers, $l,m,n$. Calculations for the chiral helical case,
$\theta_{1}=\theta_{23}=\theta_{34}=\theta_{0}$, can be found in [38].

The magic angles where flat bands appear depend on the integers $m,n$, and on
the displacement $\mathbf{D}$. To a first approximation, we can approximate
the electronic density of states of the system by an average over the values
of $\mathbf{D}$. This approximation can be expected to become exact in the
limit $\ell_{m2}/\ell_{m}\sim 1/\theta_{0}\rightarrow\infty$, and also when
the dependence of the band structure on $\mathbf{D}$ is small.

The flat bands at the magic angles reported here show a greater degree of
plasmon screening compared to tBG and thus, a greater resilience against
distortion from interaction effects. They also have a finite Berry curvature
and, typically, non zero Chern numbers, see [38].

Results. Analysis of the trilayer studied in [29]. We now analyze the
experimental situation reported in [29]. The trilayer studied there shows
angles $\theta_{12}\approx 1.42^{\circ},\theta_{23}\approx-1.88^{\circ}$,
which are well approximated by $\theta_{12}=3\theta_{0}$ and
$\theta_{23}=-4\theta_{0}$ with $\theta_{0}\approx 0.47^{\circ}$. Results for
local and average densities of states are shown in Fig. 3 (we use
$\\{t_{AA},t_{AB}\\}=\\{0.0797,0.0975\\}$ eV, and $\gamma_{0}=3.06$ eV, which
implies $v_{F}\approx 987$ Km/s, in order to compare with [29]). As shown in
the figure, the dependence of the density of states on the displacement
$\mathbf{D}$ is very weak, which suggests that the weighted average shown in
Fig. 3 (d) gives a good approximation to the exact value.

Figure 4: The tight-binding results. (a) The density of states for a trilayer
with angles $\theta_{12}=1.41^{\circ}$, $\theta_{23}=-1.88^{\circ}$. The gray
line is the result in the chiral limit. (b) The density of states for a
trilayer with angles $\theta_{12}=2^{\circ}$, $\theta_{23}=2^{\circ}$. (c) The
local density of states at energy of 7 meV in the real space for trilayer in
(a), with a zoom in of the states with displacements $\mathbf{D}$ marked by
dashed square with colors. Red and gray correspond to the color in Fig. 2(a),
and blue are in the middle region.

Figure 4(a) is the tight-binding (TB) density of states for trilayer with
$\theta_{12}\approx 1.42^{\circ},\theta_{23}\approx-1.88^{\circ}$, which is
obtained via the parameters $\gamma_{0}=3.1$ eV and $\gamma_{1}=0.43$ eV [43,
44, 45], the same as the ones in [29]. The discrepancy between the TB and
chiral limit results is due to the difference of the Fermi velocity. Moreover,
if we reduce the Fermi velocity, a flat band appears in the charge neutrality
point (not shown here). The states of the peak at the energy of 7 meV in Fig.
4(a) is highly localized with an independence of the displacement
$\mathbf{D}$, shown in Fig. 4(c), which is in agreement with the results in
Fig. 3 in the chiral limit.

We extend our discussion to another configuration, the helical trilayer
graphene with $\theta_{12}=\theta_{23}=\theta_{0}$. One magic angle is
$\theta_{0}=2^{\circ}$ [25], which has peaks located at the charge neutrality
point, shown in Fig. 4. If we consider the lattice relaxation [46, 47, 48],
the magic angle reduces to $1.57^{\circ}$ [39]. For $AAA$ helical trilayer
graphene with identical twist angles, the moiré-of-moiré structure locally
relaxes into large areas that contains only one moiré pattern. That is, a
large region of the periodic $ABA$ could be reconstructed if a strong
intralayer and interlayer potential is utilized. The states of the flat band
in the relaxed system is highly dependent on the displacement $\mathbf{D}$
(see [38] Fig. 10), which is consistent with the results in Fig. 2(a).
Furthermore, the helical trilayer graphene with $\theta_{12}=\theta_{23}$
could be an ideal platform for realization of the chiral limit due to
significant corrugation, or via heterostrain engineering [39].

Conclusions. We have considered the existence of flat bands in twisted stacks
of graphene layers where the twist angles do not allow for the definition of a
simple moiré structure, generalizing previous works [32, 34, 35, 37, 26, 41,
39, 40] which have focused on special configurations yielding exact or
approximate moiré structures in trilayers. By analyzing the chiral limit [17,
18] we conclusively demonstrate that flat bands appear for all combinations of
twist angles where the angles are close to multiples of coprime integers. Our
further extension to generic, incommensurate angles [38] connects the
previously disparate magic points of specific trilayer configurations into an
extended magic line. Building on the local approximation as developed in [37]
for the specific case of helical trilayers with nearly equal twist angles, we
accurately describe the electronic structure of generic twisted multilayers
even away from the magic angle, with our results consistent with both large
scale tight binding calculations [25] and contemporary experiments [29]. The
analysis presented here opens the way to the study of topological features in
many types of twisted multilayers, and they suggest that the rich phase
diagram found in twisted bilayer graphene may exist in many other graphene
stacks.

Note added. We have noticed ref.([54], posted a day before our paper. This
reference discusses topics related to our work. Insofar the two manuscript
overlap, the results are in agreement.

Acknowledgements. We are thankful to B. Amorim, B. A. Bernevig, and E. Castro
for helpful discussions. IMDEA Nanociencia acknowledges support from the
“Severo Ochoa” Programme for Centres of Excellence in R&D (CEX2020-001039-S /
AEI / 10.13039/501100011033). F.G. acknowledges funding from the European
Commission, within the Graphene Flagship, Core 3, grant number 881603 and from
grants NMAT2D (Comunidad de Madrid, Spain), SprQuMat (Ministerio de Ciencia e
Innovación, Spain) and financial support through the (MAD2D-CM)-MRR MATERIALES
AVANZADOS-IMDEA-NC. DCWF, MA, LP and SA acknowledge support from the Singapore
National Research Foundation Investigator Award (NRF-NRFI06-2020-0003). Z.Z.
acknowledges support funding from the European Union’s Horizon 2020 research
and innovation programme under the Marie Skłodowska-Curie grant agreement No
101034431. Numerical calculations presented in this paper have been performed
on the supercomputing system in the Supercomputing Center of Wuhan University.

Supplemental Material for:
Extended magic phase in generalized trilayer graphene

##  I Definition of periodicities.

For completeness, we specify the two relevant periodicities in a twisted
graphene trilayer with almost commensurate twist angles, $\theta_{12}$ and
$\theta_{23}$ following the steps described in [37]. We write the two angles
as:

$\displaystyle\theta_{12}$ $\displaystyle=m\theta_{0}+m\delta\theta$
$\displaystyle\theta_{23}$ $\displaystyle=n\theta_{0}-n\delta\theta$ (S1)

where $m,n$ are coprime integers, and we assume that $\theta_{0}\ll 1$ and
$\delta\theta\ll\theta_{0}$. The angles $\theta_{12}$ and $\theta_{23}$ define
moiré structures. The vectors connecting the points $\vec{K}$ and
$\vec{K}^{\prime}$ of the respective Brillouin zones are:

$\displaystyle\Delta\vec{K}_{12}$
$\displaystyle\approx\frac{4\pi}{3d}\left\\{-\frac{m^{2}(\theta_{0}^{2}+2\theta_{0}\delta\theta)}{2},m\theta_{0}+m\delta\theta\right\\}\approx
m\vec{g}_{0}+m\delta\vec{g}$ $\displaystyle\Delta\vec{K}_{23}$
$\displaystyle\approx\frac{4\pi}{3d}\left\\{-\frac{n^{2}(\theta_{0}^{2}-2\theta_{0}\delta\theta)}{2},n\theta_{0}-n\theta_{0}\right\\}\approx
n\vec{g}_{0}-n\delta\vec{g}$ (S2)

where we have assumed that $\theta_{0}\ll 1$. We obtain:

$\displaystyle\vec{g}_{0}$
$\displaystyle\approx\frac{4\pi}{3d}\left\\{\frac{m+n}{4}\theta_{0}^{2}-\frac{m-n}{2}\theta_{0}\delta\theta,\theta_{0}\right\\}\approx\frac{4\pi}{3d}\left\\{0,\theta_{0}\right\\}$
$\displaystyle\delta\vec{g}$
$\displaystyle\approx\frac{4\pi}{3d}\left\\{-\frac{m-n}{4}\theta_{0}^{2}-\frac{m+n}{2}\theta_{0}\delta\theta,\delta\theta\right\\}$

Vectors $\vec{g}_{0}$ and $\delta\vec{g}$ and three equivalent vectors rotated
by $\pm 120^{\circ}$ define the relevant periodicities of the trilayer. The
two vectors satisfy $|\delta\vec{g}|\ll|\vec{g_{0}}|$, so that, in real space,
the periodicity associated with $\delta\vec{g}$ changes much more slowly than
that asociated to $\vec{g}_{0}$. Locally, $\vec{g}_{0}$ defines a moiré
lattice with unit vector $\ell_{m}\approx d/\theta_{0}$, where $d$ is the
lattice constant of graphene.

In the following, for simplicity, we assume that
$\delta\theta\ll\theta_{0}^{2}$, although this assumption is not crucial.
Then, $\delta\vec{g}$ defines a new periodicity, with unit cell rotated
$90^{\circ}$ with respect to the unit cell defined by $\vec{g}_{0}$ and a
lattice unit of length
$\ell_{m2}\approx(4d)/(|m-n|\theta_{0}^{2})\gg\ell_{m}$. Note, finally, that
our definition of $\delta\vec{g}$, for $m\neq 1$ or $n\neq 1$, differs from
the one used in [37], although both definitions scale as
$|\delta\vec{g}|\propto\theta_{0}^{2}$.

The analysis described above can be extended to an arbitrary number, $N$, of
layers, and $N_{1}$ twist angles between nearest layers, provided that these
angles are almost commensurate, defined by $N-1$ coprime integers. Such a
calculation leads to a periodicity, $\vec{g}_{0}$, defined by an angle
$\theta_{0}$, and $N-2$ additional periodicities defined by vectors
$\delta\vec{g}_{i=1,\cdots,N-2}$ such that
$|\delta\vec{g}_{i}|\ll|\vec{g}_{0}|$. Insofar as the role of $\delta\vec{g}$
can be reduced to an overall displacement between the layers (see below) the
two definitions should give the same results.

##  II Slow periodicity as a displacement between layers.

We now elaborate on the argument that the slow periodicity, $\delta\vec{g}$,
identified in the previous section can be interpreted as a local displacement
of the top layer with respect to the bottom layer [37].

Figure S1: Sketch of the momenta of the Bloch waves included in the continuum
model, see eq.(S4), for $m=2,n=1$. The magenta arrows represent the vectors
$m\vec{g}_{i}$, and the orange arrows represent the vectors $n\vec{g}_{i}$,
see eq.(S8).

The interlayer tunneling part of the continuum Hamiltonian [1, 2], in real
space, is:

$\displaystyle{\cal H}_{T}$ $\displaystyle={\cal H}_{T}^{12}+{\cal
H}_{T}^{23}$ $\displaystyle{\cal H}_{T}^{12}$
$\displaystyle=\sum_{j=1,2,3}T_{j}e^{im\vec{g}_{j}\cdot\vec{r}}$
$\displaystyle{\cal H}_{T}^{23}$
$\displaystyle=\sum_{j=1,2,3}T_{j}e^{in\vec{g}_{j}\cdot\vec{r}}$ (S4)

where $T_{j}$ define inter sublattice $2\times 2$ matrices:

$\displaystyle T_{j}$
$\displaystyle\equiv\left(\begin{array}[]{cc}t_{AA}&t_{AB}e^{i\frac{2\pi(j-1)}{3}}\\\
t_{AB}e^{-i\frac{2\pi(j-1)}{3}}&t_{AA}\end{array}\right)$ (S7)

and:

$\displaystyle\vec{g}_{1}$
$\displaystyle=\frac{\ell_{m}}{\sqrt{3}}\left\\{0,1\right\\}$
$\displaystyle\vec{g}_{2}$
$\displaystyle=\frac{\ell_{m}}{\sqrt{3}}\left\\{\frac{\sqrt{3}}{2},-\frac{1}{2}\right\\}$
$\displaystyle\vec{g}_{3}$
$\displaystyle=\frac{\ell_{m}}{\sqrt{3}}\left\\{-\frac{\sqrt{3}}{2},-\frac{1}{2}\right\\}$
(S8)

The momenta of the Bloch waves used in the calculation are shown in Fig. S1.

We can generalize the tunneling Hamiltonian to the case where the moiré
pattern of layers 1 and 2 is displaced by a shift $\vec{\mathbf{D}}_{1}$ and
the moiré pattern between layers 2 and 3 is displaced by a shift
$\vec{\mathbf{D}}_{3}$ [33]. In the initial Hamiltonian, the rotation is
around a point with $AA$ stacking between layers 1 and 2, and $AA$ stacking
between layers 2 and 3, that is, $AAA$ stacking. The shifts lead to the
Hamiltonian:

$\displaystyle{\cal H}_{T}$ $\displaystyle={\cal H}_{T}^{12}+{\cal
H}_{T}^{23}$ $\displaystyle{\cal H}_{T}^{12}$
$\displaystyle=\sum_{j=1,2,3}T_{j}e^{im\vec{g}_{j}\cdot\left(\vec{r}-\vec{\mathbf{D}}_{1}\right)}$
$\displaystyle{\cal H}_{T}^{23}$
$\displaystyle=\sum_{j=1,2,3}T_{j}e^{in\vec{g}_{j}\cdot\left(\vec{r}-\vec{\mathbf{D}}_{3}\right)}$
(S9)

Each Bloch wave in the calculation, see Fig. S1, can be modified by a gauge
factor, $e^{i\vec{G}\cdot\vec{\mathbf{D}^{\prime}}}$, where $\vec{G}$ is the
momentum of the wave. This gauge transformation induces a phase in the
interlayer tunnelings, which is equivalent to a redefinition of the shifts:

$\displaystyle\vec{\mathbf{D}}_{1}^{\prime}$
$\displaystyle=\vec{\mathbf{D}}_{1}-\vec{\mathbf{D}}^{\prime}$
$\displaystyle\vec{\mathbf{D}}_{3}^{\prime}$
$\displaystyle=\vec{\mathbf{D}}_{3}-\vec{\mathbf{D}}^{\prime}$ (S10)

This result is independent of the integers $m$ and $n$, and of the angle
$\theta_{0}$. An equal shift, $\vec{\mathbf{D}}_{1}=\vec{\mathbf{D}}_{3}$, of
the two moiré patterns can be gauged away. This implies that the calculated
electronic structure is only dependent on the value of
$\vec{\mathbf{D}}_{1}-\vec{\mathbf{D}}_{3}=\mathbf{D}$.

Figure S2: Sketch of the real space unit cell structure for (a) trilayer with
$\\{m,n\\}=\\{1,2\\}$ and (c) trilayer with $\\{m,n\\}=\\{3,-4\\}$. The black
hexagon defines the large unit cell, associated with moiré length $\ell_{m2}$.
The red hexagons show the moiré structure associated with angle $\theta_{0}$
and moiré length $\ell_{m}$. Blue and green hexagons show moirés associated
with angles $\theta_{12}$ and $\theta_{23}$. Note that the moiré lattices
associated to $\theta_{0}$ and $\theta_{12}$ coincide. Special points in the
momentum space unit cell with the $AAA$ electronic structure for (b)
$\\{m,n\\}=\\{1,2\\}$ and (d) $\\{m,n\\}=\\{3,-4\\}$ cases.

The two periodicities defined by $\vec{g}_{0}$ and $\delta\vec{g}$ in the
previous section modify the continuum Hamiltonian, Eq. (S4), leading to:

$\displaystyle{\cal H}_{T}^{12}$
$\displaystyle=\sum_{j=1,2,3}T_{j}e^{im\vec{g}_{0j}\cdot\vec{r}}e^{im\delta\vec{g}_{j}\cdot\vec{r}}$
$\displaystyle{\cal H}_{T}^{23}$
$\displaystyle=\sum_{j=1,2,3}T_{j}e^{in\vec{g}_{0j}\cdot\vec{r}}e^{-in\delta\vec{g}_{j}\cdot\vec{r}}$
(S11)

The vectors $\delta\vec{g}_{j}$ define a periodicity with unit vectors in real
space:

$\displaystyle\\{\vec{D}_{1},\vec{D}_{3}\\}$
$\displaystyle=\frac{2\pi}{3d|\delta\vec{g}_{i}|}\\{\vec{R}_{1},\vec{R}_{2}\\}\approx\frac{1}{\theta_{0}^{2}}\\{\vec{R}_{1},\vec{R}_{2}\\}$
(S12)

where $\vec{R}_{1}$ and $\vec{R}_{2}$ are the lattice vectors of an untwisted
graphene layer.

We now assume that the arguments $\delta\vec{g}_{j}\cdot\vec{r}$ vary slowly
with respect the functions $\vec{g}_{0j}\vec{r}$, and assume that the value of
$\vec{r}=\vec{\mathbf{D}}$ in $\delta\vec{g}_{j}\cdot\vec{r}$ is constant.
Then, using Eq. (S9) we obtain that the local Hamiltonian is equivalent to
displacements of the top and bottom layers by
$\pm(d/\ell_{m2})R_{\pi/2}(\vec{\mathbf{D}})$, where $R_{\theta}$ is a
rotation operator.

For $m\neq 1$ or $n\neq 1$ the loops in momentum space which connect a
momentum in the central layer to itself through a path which visits both the
top and bottom layers require at least $m\times n$ steps. The increased
periodicity of these loops for $m\neq 1$ or $n\neq 1$ imply that the band
structure in real space repeats itself for a displacement periodicity defined
by the unit vectors:

$\displaystyle\\{\tilde{\vec{D}}_{1},\tilde{\vec{D}}_{3}\\}$
$\displaystyle=\frac{1}{mn}\\{\vec{D}_{1},\vec{D}_{3}\\}$ (S13)

This result implies that the actual periodicity in real space is given by
$\ell_{m2}/(mn)$. A sketch of the different regions in the unit cell, and the
points where the band structure repeats itself is shown in Fig. S2.

Figure S3: Results for the magic angles in the chiral limit and
$\theta_{12}=2\theta_{0},\theta_{23}=\theta_{0}$. (a) \- (d): The local DOS
for a magic angle $\theta_{0}=1.163^{\circ}$ with three displacements and
their average, as in Fig. 3 of the main text. (e) - (h): As in the top panel,
but for a magic angle $\theta_{0}=1.14^{\circ}$.  Figure S4: As in Fig. S3 but
for $\theta_{12}=3\theta_{0},\theta_{23}=2\theta_{0}$. (a) - (d): The local
DOS for a magic angle $\theta_{0}=0.6^{\circ}$ with three displacements and
their average as in Fig. 3 of the main text. (e) - (h): Results for a magic
angle $\theta_{0}=0.59^{\circ}$.

Densities of states (DOS) for two choices of the magic angle in the chiral
limit, and $\\{m,n\\}=\\{2,1\\}$ are shown in Fig. S3. Similar results for
$\\{m,n\\}=\\{3,2\\}$ are shown in Fig. S4.

##  III Extension to multilayers with $N>3$.

The analysis discussed previously can be extended to multilayers with any
number of layers greater than three. Results for the density of states in the
chiral limit of an helical stack of four layers
($\theta_{12}=\theta_{23}=\theta_{34}=\theta_{0}$) at different magic angles
are shown in Fig. S5.

It is also instructive to generalize the alternative method discussed in [26].
The analysis addresses the magic angles, in the chiral limit, of a stack of
$N$ layers, where each layer is twisted by the same angle, $\theta$, and in
the same direction, with respect to the previous layer.

The magic angles for $N$ layers are determined by an $N-$component spinor
which is the solution of $N$ first order differential equations, similar to
the case of a 2-component spinor, $\\{\psi_{\vec{k}},\bar{\psi}_{\vec{k}}\\}$,
which determines the magic angles in twisted bilayer graphene. The solution
found in [26] describes the three component spinor for the trilayer in terms
of three combinations of products of combinations of these wavefunctions:

$\displaystyle\Psi_{3}$
$\displaystyle\equiv\left\\{\begin{array}[]{c}\psi_{\vec{k}_{1}}\psi_{\vec{k}_{2}}\\\
\frac{1}{\sqrt{2}}\left(\psi_{\vec{k}_{1}}\bar{\psi}_{\vec{k}_{2}}+\bar{\psi}_{\vec{k}_{1}}\psi_{\vec{k}_{2}}\right)\\\
\bar{\psi}_{\vec{k}_{1}}\bar{\psi}_{\vec{k}_{2}}\end{array}\right.$ (S17)

A similar equation can be written for an $N-$component spinor:

$\displaystyle\Psi_{3}$
$\displaystyle\equiv\left\\{\begin{array}[]{c}\prod_{i=1}^{i=N}\psi_{\vec{k}_{i}}\\\
\frac{1}{\sqrt{N}}\sum_{j=1}^{j\leq N}\bar{\psi}_{\vec{k}_{j}}\prod_{i=1,i\neq
j}^{i=N}\psi_{\vec{k}_{i}}\\\ \frac{1}{\sqrt{N(N-1)}}\sum_{j=1,k=1,k\neq
j}^{j=N,k=N}\bar{\psi}_{\vec{k}_{j}}\bar{\psi}_{\vec{k}_{k}}\prod_{i=1,i\neq
j,i\neq k}^{i=N}\psi_{\vec{k}_{i}}\\\ \cdots\\\
\prod_{i=1}^{i=N}\bar{\psi}_{\vec{k}_{i}}\end{array}\right.$ (S23)

We defined $\alpha_{2}=[v_{F}(4\pi)/(2d\sin\theta_{M}/2)]/t_{AB}$, where
$\theta_{M}$ is the magical angle for a twisted bilayer. Then, Eq. (S23) has a
solution for:

$\displaystyle\alpha_{N}$ $\displaystyle=\frac{4\pi
v_{F}}{dt_{AB}}\left[2\sin\left(\frac{\theta_{M}^{N}}{2}\right)\right]=\alpha_{2}\left(\frac{1}{(N-1)!}\right)^{\frac{1}{N-1}}$
(S24)

We have checked that this equation gives the correct solution for
$\alpha_{3}=\alpha_{2}/\sqrt{2}$ and $\alpha_{4}=\alpha_{2}\times
1/\sqrt[3]{6}$ agrees with our numerical results for a stack of four layers.
For $N\rightarrow\infty$ we find:

$\displaystyle\lim_{N\rightarrow\infty}\theta_{M}^{N}$
$\displaystyle\approx\frac{dt_{AB}}{4\pi v_{F}}\frac{e}{N-1}$ (S25)

where $e$ is Euler’s number. This result implies that a small shear force
applied to the top and bottom layers of a graphite stack can lead to a peak at
the density of states near the Fermi level.

Figure S5: DOS for magic angles in the chiral limit of four layer stacks with
$\theta_{12}=\theta_{23}=\theta_{34}=\theta_{0}$ (helical arrangement). (a)
$AAAA$ stacking, $\theta_{0}=0.55^{\circ}$. (b) $AABB$ stacking,
$\theta_{0}=0.395^{\circ}$. (c) $ABAB$ stacking, $\theta_{0}=0.299^{\circ}$.

##  IV Validity of the local Hamiltonian expansion

The moiré-of-moiré length scale emerges due to the misalignment of the moiré
patterns between layers 1 and 2, and layers 2 and 3,
$(\theta_{12}+\theta_{23})/2$, as shown in the Fig. 1(c) of the main text.
Starting from a particular momentum point in the second layer, we may reach
other points in the second layer by hopping through either the first or third
layer. Hoppings between only 2 of the layers allow us to define a momentum
space moiré lattice, a subset of which are shown as red points along the black
lines. The dotted line marks the closest approach between middle layer states
connected by hoppings that only pass through either the first or third layer,
which defines the moiré-of-moiré modulation lengthscale. The sketch
specifically shows the situation for ${m,n}={3,4}$, $\theta_{0}=1^{\circ}$.

The distance between a momentum state reached by $a$ hops through the first
layer and another state reached by $b$ hops through the third layer is
$\sqrt{3}k_{\mathrm{D}}P_{ab}$, where

$\displaystyle P_{ab}^{2}$
$\displaystyle=(a\theta_{12})^{2}+(b\theta_{23})^{2}-2ab\left\lvert\theta_{12}\theta_{23}\right\rvert\cos\left(\frac{\theta_{12}+\theta_{23}}{2}\right)$
$\displaystyle\approx(a\left\lvert\theta_{12}\right\rvert-b\left\lvert\theta_{23}\right\rvert)^{2}+ab\left\lvert\theta_{12}\theta_{23}\right\rvert\left(\frac{\theta_{12}+\theta_{23}}{2}\right)^{2}$
$\displaystyle\rightarrow(am-b\left\lvert
n\right\rvert)^{2}\theta_{0}^{2}+abm\left\lvert
n\right\rvert\left(\frac{m+n}{2}\right)^{2}\theta_{0}^{4},$ (S26)

and we take the small angle approximation in the second line and commensurate
angles ($m\geq 0$ WLOG) in the third line. The moiré-of-moiré lengthscale is
then $\l_{m2}=2\pi/\sqrt{3}k_{\mathrm{D}}\mathrm{min}(P_{ab})$, where we
minimise $P_{ab}$ over all integers $a\geq 0$ and $b\geq 0$, except for the
trivial solution $P_{00}=0$. For small commensurate angles with $m$ and $n$
small, this generally occurs at $P_{\left\lvert n\right\rvert
m}=\frac{\left\lvert mn(m+n)\right\rvert}{2}\theta_{0}^{2}$, however, if $m$
and $n$ are large and of the same sign, $P_{ab}$ may be minimised at some
other $a\neq\left\lvert n\right\rvert$, $b\neq m$. For example, at
$\theta_{12}=1.0^{\circ}$, $\theta_{23}=0.9^{\circ}$, $P_{1,1}<P_{9,10}$. For
completeness, we note that this description only considers lattice points
lying on the two lines as shown in the sketch, when in fact there are others
in the 2D plane. Provided the angles and thus the misalignment is small, the
points on these lines remain the closest to each other.

Numerically calculated results for $l_{m2}$ are shown in the Fig. 1(d) of the
main text. As expected, $l_{m2}$ diverges along the line
$\theta_{12}=-\theta_{23}$, and is generally larger for
$\theta_{23}/\theta_{12}$ negative than for $\theta_{23}/\theta_{12}$ positive
due to the reduced misalignment of the moiré patterns. $l_{m2}$ is also
particularly pronounced at commensurate angles with simple ratios, where it is
easy to minimise the first term of Eq. (S26) without the second term getting
too large. For experimental samples smaller than $l_{m2}$, a particular local
Hamiltonian expansion is expected to be a good description of the bands while
for samples much larger than $l_{m2}$, it is necessary to average over
expansion centers, or equivalently, over $\mathbf{D}_{1}-\mathbf{D}_{3}$. We
note that the form for $l_{m2}$ derived and calculated here agrees with that
previously derived [37] for the special case of
$\theta_{23}/\theta_{12}\approx 1$.

Neglecting the slow periodicity in the local expansion, as is done in the
numerical calculations of us and [37], neglects the misalignment and so forces
the black lines of the upper panel sketch to coincide. In the case of
commensurate angles, this leads to $\mathrm{min}(P_{ab})=0$ and thus $l_{m2}$
diverging, leading us to ignore the modulation of
$\mathbf{D}_{1}-\mathbf{D}_{3}$.

##  V Finding the magic surface in trilayer configuration space

Our numerical calculations have shown that a magic angle pair may be found
seemingly for any numerically accessible (commensurate) twist angle ratio, and
so the natural question would be whether this may be extended to the continuum
in angle-pair-space. We start by writing the trilayer Hamiltonian in layer
space,

$\displaystyle H=\begin{pmatrix}D_{1}&{\cal H}_{T}^{12}&0\\\ {\cal
H}_{T}^{12\dagger}&D_{2}&{\cal H}_{T}^{23}\\\ 0&{\cal
H}_{T}^{23\dagger}&D_{3}\end{pmatrix}_{L},$ (S27)

with $D_{l}$ the Dirac cones of layer $l$ and the ${\cal H}_{T}^{ij}$ as in
Eq. (S9). The subscript $L$ denotes layer space, and each “matrix element”
here is itself a matrix in sublattice space. In the chirally symmetric limit,
$t_{AA}=0$, a simple permutation allows us to rewrite the Hamiltonian in the
dimensionless form

$\displaystyle
H_{\mathrm{chiral}}/t_{AB}=\begin{pmatrix}0&\mathcal{D}^{\dagger}\\\
\mathcal{D}&0\end{pmatrix}_{\mathrm{AB}},$ (S28)

where the subscript AB indicates sublattice space and the zeroes on the
diagonals are the natural consequence of using the chirally symmetric model.
The layer space submatrix is given by

$\displaystyle\mathcal{D}=\begin{pmatrix}-i\frac{\hbar
v_{\mathrm{F}}}{t_{AB}}\bar{\partial}&U_{12}(\mathbf{r})&0\\\
U_{12}(-\mathbf{r})&-i\frac{\hbar
v_{\mathrm{F}}}{t_{AB}}\bar{\partial}&U_{23}(\mathbf{r})\\\
0&U_{23}(-\mathbf{r})&-i\frac{\hbar
v_{\mathrm{F}}}{t_{AB}}\bar{\partial}\end{pmatrix}_{L},$ (S29)

where $\bar{\partial}_{i}=\tfrac{1}{2}(\partial_{x}+i\partial_{y})$ and their
relative twists have been compensated for by corresponding countertwists of
the spinors.
$U_{12}(\mathbf{r})=\sum_{j=1,2,3}e^{i\frac{2\pi(j-1)}{3}}e^{im\vec{g}_{j}\cdot\left(\vec{r}-\vec{\mathbf{D}}_{1}\right)}$
is the moiré potential between layers 1 and 2, and similarly for
$U_{23}(\mathbf{r})$ replacing $m$ with $n$ and $\mathbf{D}_{1}$ with
$\mathbf{D}_{3}$.

The off-diagonal elements of the layer space submatrix have the form

$\displaystyle\mathcal{D}_{off}=\begin{pmatrix}0&a&0\\\ b&0&c\\\
0&d&0\end{pmatrix}_{L},\quad\\{a,b,c,d\\}\in\mathbb{C},$ (S30)

which may be collected into a single off-diagonal through the similarity
transformation

$\displaystyle S^{-1}\mathcal{D}_{off}S=\begin{pmatrix}0&0&0\\\
0&0&a\sqrt{1+\frac{cd}{ab}}\\\
0&b\sqrt{1+\frac{cd}{ab}}&0\end{pmatrix}_{L^{\prime}},\quad
S=\frac{1}{\sqrt{2}}\begin{pmatrix}-\frac{c}{b}&0&\frac{a}{d}\\\
0&\sqrt{2}&0\\\ 1&0&1\end{pmatrix},$ (S31)

where the subscript $L^{\prime}$ reminds us of the nontrivial mixing of the
layers. We are thus able to re-express the interfering double-moiré problem
into a single effective moiré potential. We note that the transformation $S$
is not unitary in general and so does not conserve the inner products of
spinors, but does conserve the bandstructure. We further note that there is
some freedom to choose the form of the single off-diagonal element; our choice
most easily connects to previous work [19] where in the case of alternating
twist angles, $U_{12}(\mathbf{r})=U_{23}(-\mathbf{r})$ and so $S$ is unitary
and the form factor $\sqrt{1+\frac{cd}{ab}}=\sqrt{2}$, resulting in a constant
scaling of the bilayer moiré potential and thus of the magic angle.

Applying the transformation of Eq. (S31) to Eq. (S29) then allows us to obtain
a single effective moiré potential, from which an estimate of the magic angle
pair may be made by comparison with the bilayer moiré potential, a nontrivial
extension of previous work [19] which focused on the special case of
alternating twist angle multilayers.

Our estimate for the magic angle and the error of our estimate is calculated
from the average and standard deviation of the form factor
$\sqrt{1+\frac{cd}{ab}}$ over space. While the form of the effective moiré
potential we have chosen does indeed have singularities e.g. at the origin,
these diverge as $\tfrac{1}{x+iy}$ and so the integral over the plane
converges, which is unsurprising as it is possible to choose a form of
effective moiré potential with no singularities at all, though in that case
the analogy to twisted bilayer graphene is weaker.

Figure S6: Top: Predicted magic angle pairs
$(\theta_{12}^{\mathrm{M}},\theta_{23}^{\mathrm{M}})$ for $\mathbf{D}=0$.
Shaded regions denote the error in the magic angle pair prediction. The
configurations studied in [26] and [19] are indicated with a black and purple
diamond respectively, and those corresponding to our numerical calculations
are marked with black diamonds. Dotted lines of constant
$\theta_{12}/\theta_{23}$ are marked as guides to the eye. Bottom: Predicted
range of magic angles over all $\mathbf{D}$. The shaded region indicates pairs
of angles for which there exists some $\mathbf{D}$ that would give flat bands.
The configurations studied in [49, 50, 51] are marked in green, red and blue
respectively.

Results are shown in the upper panel of Fig. S6 for zero shift,
$\mathbf{D}=0$, where we indicate with a red line the predicted magic angle
pairs according to the average-value analysis described above, with the
uncertainty in the prediction shown by the shaded pink region. The alternating
twist case previously studied [19] is marked with a blue diamond, and we note
that this is the only point where the analysis is exact. We mark also our
numerically calculated magic angle pairs with black diamonds, and note the
excellent agreement with the more general analysis, giving us confidence that
the results are reliable even for the numerically inaccessible case of
incommensurate angles. Our analysis predicts two branches of magic angle pairs
from which an infinite set descends, as in the bilayer case [18]. At most two
sets of numerically calculated magic angles per angle ratio are shown to
connect with the two branches, though other values have been found.

We repeat the calculation over all moiré shifts to obtain a range of magic
angle pairs, shown in the lower panel of Fig. S6. The magic angles vary
continuously with shift and so the manifold of magic configurations is a
three-dimensional surface within the four-dimensional trilayer configuration
space of two angles and a two-dimensional relative shift. This surface divides
the configuration space into unconnected sectors, as seen in the lower panel
of Fig. S6 where it is not possible to move from $\theta_{12}=\theta_{23}=0$
to $\theta_{12}=\theta_{23}=2^{\circ}$ without passing through the magic
surface. We remind the reader here that the plot is symmetric about the
origin. This may have consequences on the topology of the bands as a function
of trilayer configuration, as magic configurations may arise as the
“midpoints” of band crossings, with the flat bands flattened by repulsion of
the crossing bands. Previous work on alternating-twist magic angle trilayer
graphene [49], with trilayer configuration
$(\theta_{12}=1.57^{\circ},\theta_{23}=1.57^{\circ},\mathbf{D}=0)$, a twisted
monolayer–$AA$-bilayer system [50], configuration
$(\theta_{12}=1.22^{\circ},\theta_{23}=0,\mathbf{D}=0)$, and a moiré
quasiperiodic crystal [51], with configuration
$(\theta_{12}=1.41^{\circ},\theta_{23}=1.88^{\circ},\mathbf{D}\>\mathrm{unknown})$,
may be placed as single points within this extended magic phase.

##  VI Tight-binding results

We use a round disk method to construct the twisted graphene trilayers with
arbitrary twist angles. The two independent twist angles $\theta_{12}$ and
$\theta_{23}$ are chosen to be the rotation of the second layer 2 relative to
the first layer 1 and the rotation of the third layer 3 relative to the second
layer 2, respectively. The rotation origin is chosen at an atom site. We use a
twist angel pair ($\theta_{12}$, $\theta_{23}$) as the notation for different
twist angle configurations. Positive (negative) values of the twist angle
denotes counterclockwise (clockwise) rotations. The sample with (-$\theta$,
$\theta$) has a mirror symmetry with the middle layer as the mirror plane. To
calculate the property of these large scale systems with arbitrary twist
angles, we construct the system in a large round disk. The radius of the disk
should be set sufficiently large to rid the effects of edge states [44]. In
the actual calculation, the disk with radius of $172.2$ nm ($700d$ being $d$
the lattice constant of graphene) and contains 10 million carbon atoms are
utilized for the twist angles investigated in this work. Figure S7 shows the
($21.8^{\circ}$, $21.8^{\circ}$) configuration of twisted trilayer graphene
with $AAA$ stacking.

Figure S7: Schematic of the $AAA$ trilayer graphene with
$\theta_{12}=\theta_{23}=\theta_{0}=21.8^{\circ}$. The $\theta_{12}$ and
$\theta_{23}$ are the twist angles between layer 1 and layer 2, layer 2 and
layer 3, respectively. The orange and purple rhombi are the moiré unit cells
with $\theta_{12}$ and $\theta_{23}$, respectively, which has length
$\ell_{m}=0.65$ nm. The black rhombus is the unit cell of the moiré-of-moiré
with $\ell_{m2}=1.7$ nm. The green square shows an $ABA$ stacking region.

A parameterized full tight-binding (TB) model is used. In the tight-binding
model, only $p_{z}$ orbitals are taken into account, we construct the
Hamiltonian of the twisted graphene trilayers as

$H=\sum_{i}\epsilon_{i}|i\rangle\langle i|+\sum_{\langle
i,j\rangle}t_{ij}|i\rangle\langle j|,$ (S32)

where $|i\rangle$ is the $p_{z}$ orbital located at $\mathbf{r}_{i}$, and
$\langle i,j\rangle$ is the sum over index $i$ and $j$ with $i\neq j$.
According to the Slater-Koster (SK) formalism, the hopping integral $t_{ij}$,
interaction between two $p_{z}$ orbitals located at $\mathbf{r}_{i}$ and
$\mathbf{r}_{j}$ has the form [45]

$t_{ij}=n^{2}V_{pp\sigma}(r_{ij})+(1-n^{2})V_{pp\pi}(r_{ij}),$ (S33)

where $r_{ij}=|\mathbf{r}_{j}-\mathbf{r}_{i}|$ is the distance between $i$ and
$j$ sites, with $n$ as the direction cosine along the direction $\bm{e_{z}}$
perpendicular to the graphene layer . The Slater and Koster parameters
$V_{pp\pi}$ and $V_{pp\sigma}$ follow

$\displaystyle V_{pp\pi}(r_{ij})=-t_{0}e^{q_{\pi}(1-r_{ij}/a)}F_{c}(r_{ij}),$
(S34) $\displaystyle
V_{pp\sigma}(r_{ij})=t_{1}e^{q_{\sigma}(1-r_{ij}/h)}F_{c}(r_{ij}),$

where $a=1.42$ Å and $h=3.349$ Å are the nearest in-plane and out-of-plane
carbon-carbon distance, respectively, $t_{0}=3.1$ eV and $t_{1}=0.43$ eV are
the TB hopping parameters taken from the moiré quasiperiodic crystal paper
[51]. The parameters $q_{\sigma}$ and $q_{\pi}$ satisfy
$\frac{q_{\sigma}}{h}=\frac{q_{\pi}}{a}=2.218$ Å-1, and the smooth function is
$F_{c}(r)=(1+e^{(r-r_{c})/l_{c}})^{-1}$, in which $l_{c}$ and $r_{c}$ are
chosen as $0.265$ and 5.0 Å, respectively. In the twisted graphene trilayer
calculations, we only consider the interlayer hoppings between adjacent
layers. An open boundary condition is used in the disk model. The electronic
properties of the twisted graphene trilayer with arbitry angles
($\theta_{12}$,$\theta_{23}$) are calculated by the tight-binding propagation
method implemented in the TBPLaS simulator [43]. For example, the detailed
formula of the density of states is

$D(\varepsilon)=\frac{1}{2\pi
S}\displaystyle\sum_{p=1}^{S}\int_{-\infty}^{\infty}e^{i\varepsilon
t}\langle\varphi_{p}(0)|e^{-iHt}|\varphi_{p}(0)\rangle dt,$ (S35)

where $|\varphi_{p}(0)\rangle$ is one initial state which is the random
superposition of all basis states, $S$ is the number of random initial states.
The distribution of states in real space can be obtained by calculating the
quasieigenstates [43] (a superposition of degenerate eigenstates with certain
energy). The quasieigenstates has the expression:

$|\Psi(\varepsilon)\rangle=\frac{1}{\sqrt{\sum_{n}|A_{n}|^{2}\delta(\varepsilon-
E_{n})}}\sum_{n}A_{n}\delta(\varepsilon-E_{n})|n\rangle,$ (S36)

where $A_{n}$ are random complex numbers with $\sum_{n}|A_{n}|^{2}=1$, $E_{n}$
is the eigenvalue and $|n\rangle$ is the corresponding eigenstate. The local
density of states (LDOS) mapping calculated from the quasieigenstates is
highly consistent with the experimentally scanning tunneling microscopy
$\mathrm{d}I/\mathrm{d}V$ mapping. We also consider the lattice relaxation
effect on the electronic properties of twisted trilayer graphene. In the round
disk sample, the edge carbon atoms possessing dangling bond are passivated by
placing in-plane hydrogen atoms to saturate the dangling $\sigma$ edge bonds.
The carbon-hydrogen bond length is assumed to be 0.1 nm. We employ the
classical molecular dynamics simulation package LAMMPS to do the full (both
in-plane and out-of-plane) lattice relaxation [46]. The intralayer C-C and C-H
interactions and interlayer C-C interactions are simulated with REBO [47] and
kolmogorov/crespi/z version of Kolmogorov-Crespi (KC) [48] potentials,
respectively.

Figure S8: The TB DOS of $AAB$ trilayer graphene with twist angles
$\theta_{12}=\theta_{23}=\theta_{0}$. The hopping parameters are $t_{0}=2.8$
eV and $t_{1}=0.44$ eV. The displacement $\mathbf{D}$ is zero.

We focus on the helical twisted trilayer graphene configuration. Firstly, we
consider the particular case that the two twist angles are equal. In the TB
calculation, we set $t_{0}=2.8$ eV and $t_{1}=0.44$ eV, which give a magic
angle of $1.05^{\circ}$ in the twisted bilayer graphene case. The evolution of
the density of states with twist angles of $AAB$ twisted trilayer graphene is
shown in Fig. S8. Note the discontinuity of the van Hove singularity (VHS)
evolves with the twist angle is due to the low resolution of the twist angle,
which will not change our conclusion. The system can be considered as a $AAA$
trilayer with a displacement of the carbon-carbon distance $a$ between the top
and bottom layers. The tendency is similar to the $AAA$ trilayer graphene case
[25, 36]. That is, with $\theta_{0}$ decreases, the van Hove singularity gap
(the energy difference between the first van Hove singularity on the valence
and conduction bands) narrows and reaches minimum value at around $2^{\circ}$.
This result is reasonable since the moiré unit cell with periodicity
$\ell_{m2}$ contains all the configurations, for instance, the $AAA$ and $ABA$
stackings. In this case, the first ”magic angle” is $2^{\circ}$. The VHSs
merge again at the angle $\theta_{0}=1^{\circ}$. In the magic angle of
trilayer in chiral limit in Fig. (3) of the main text, there is a magic angle
at $\theta_{0}=0.72^{\circ}$. In the case of $\theta_{12}=\theta_{0}$,
$\theta_{23}=2\theta_{0}$, shown in Fig. S9, the first magic angle is
$\theta_{0}=1.25^{\circ}$, which is consistent with the result in the
continuum limit in Ref. [52].

Figure S9: The TB DOS of $AAB$ trilayer graphene with $\theta_{0}$ where twist
angles $\theta_{12}=\theta_{0}$, $\theta_{23}=2\theta_{0}$. The hopping
parameters are $t_{0}=2.8$ eV and $t_{1}=0.44$ eV. The displacement
$\mathbf{D}$ is zero. Curves are vertically shifted for clarity. Figure S10:
The lattice relaxation effect. (a) The TB DOS of $AAA$ trilayer graphene with
magic twist pair $\theta_{12}=\theta_{23}=\theta_{0}=1.57^{\circ}$. (b) The
local DOS at zero energy in real space with a zoom-in of the states at
different displacement $\mathbf{D}$. The colors of the squares correspond to
the ones in Fig. 2(a) of the main text. The hopping parameters are $t_{0}=2.8$
eV and $t_{1}=0.44$ eV.

As discussed in the Ref. [25], in the case of $AAA$ twisted trilayer graphene
with $\theta_{12}=\theta_{23}$, the first magic angle decreases from
$2^{\circ}$ to $1.57^{\circ}$ when the lattice relaxation is considered. A
sharp peak is located at the charge neutrility point, shown in Fig. S10(a).
The states of the peak in the real space are different for different
displacement $\mathbf{D}$ are different, which is in agreement with the result
in Fig. 2(a) of the main text. Such state feature can be detected by the
scanning tunneling microscopy in experiment. However, For the relaxation of
the system with the REBO and KC potentials, we could not find a strong
reconstruction of the structure with the $ABA$ region relax to large area and
clear domain wall region, which may achieve via a stronger intralyer and
interlayer potential [39].

##  VII Dielectric screening in trilayers

It has been widely shown that the dielectric response for multi-twisted
materials increases drastically as the bands become flat [53]. The plasmon
screening on electron-electron interaction is believed to have a significant
effect on the strong correlation phenomena. One potential advantage of magic
multilayers over tBG is that the calculated noninteracting bandstructures are
expected to be more resilient against Hartree effects due to the higher
dielectric function from plasmon screening.

We may verify this by calculating the static component of the dielectric
function. For simplicity, we ignore Umklapp processes, which would give a
matrix structure to the polarization, with entries labelled by moiré
reciprocal lattice vectors, $\vec{G},\vec{G}^{\prime}$. Within the framework
of random phase approximation (RPA), the dynamic dielectric function can be
calculated via

$\epsilon(\mathbf{k},\omega)=1-V(\mathbf{k})\Pi(\mathbf{k},\omega)$ (S37)

where $V(\mathbf{k})=2\pi e^{2}/\kappa k$ is the Fourier component of the two-
dimensional Coulomb interaction, with $\kappa$ being the background dielectric
constant. We choose $\kappa=4$ in our calculation, in effect assuming the
substrate is hexagonal boron nitride (hBN). $\Pi(\mathbf{k},\omega)$ is the
dynamic polarization function given by

$\Pi(\mathbf{k},\omega)=2\sum_{\mathbf{q},\xi}\sum_{m,n}\frac{\left(f_{\mathbf{q}+\mathbf{k},\xi}^{n}-f_{\mathbf{q},\xi}^{m}\right)F_{\mathbf{q},\mathbf{q}+\mathbf{k},\xi}^{nm}}{E_{\mathbf{q}+\mathbf{k},\xi}^{n}-E_{\mathbf{q},\xi}^{m}-\omega-i0}$
(S38)

here $m,n$ are band index, $\xi$ is valley index, $f_{\mathbf{q}}^{m}$ is the
Fermi-Dirac distribution, and
$F_{\mathbf{q},\mathbf{q}+\mathbf{k},\xi}^{nm}=\left|\psi_{n,\mathbf{q}+\mathbf{k},\xi}^{\dagger}\psi_{m,\mathbf{q},\xi}\right|^{2}$
is the form factor of two different Bloch states.

Figure S11: Static component of the dielectric function,
$\epsilon(\mathbf{k},\omega=0)$ evaluated over the Brillouin zone for
alternating-twist magic angle trilayer graphene at $AAA$ (red) and $AAB$
(blue) stacking, compared to that of magic angle tBG (black). The value of
$\epsilon$ for trilayers is greatly increased compared to that of tBG
everywhere except the $\Gamma$ point, indicating an overall greater degree of
plasmon screening in trilayers and thus greater robustness of the calculated
noninteracting bandstructures to Hartree effects.

Results are shown in Fig. S11, which shows $\epsilon(\mathbf{k},\omega=0)$
over the moiré Brillouin zone for alternating-twist magic angle trilayer
graphene in the $AAA$ and $AAB$ configurations, and for tBG. We choose to
compare with alternating twist trilayers to remove moiré-of-moiré effects, and
to be able to make direct comparisons between the Brillouin zones of the two
materials. The value of $\epsilon$ is seen to be $\mathbf{D}$ dependent, as
expected from the $\mathbf{D}$ dependance of the bandstructure, and is
noticeably higher than that of tBG everywhere except at $\Gamma$, indicating
that there is a greater degree of plasmon screening in the trilayer, as
expected, and that the calculated noninteracting bands are therefore expected
to be more robust against Coulomb interaction effects than those of tBG.

##  VIII Electronic bands and Berry curvature.

The flat electronic bands at magic angles studied in the main text show finite
Berry curvatures, and, generally, non zero Chern numbers, see also [41, 39].
Examples are given in given in Fig.[S12] and Fig.[S13].

Figure S12: Top: Low energy electronic bands for an $ABA$ trilayer with
$\\{m,n\\}=\\{3,-4\\}$ and
$\\{\theta_{12},\theta_{23}\\}=\\{0.85^{\circ},-1.14^{\circ}\\}$. Using the
parameters described in the main text, two flat bands are obtained at
$\epsilon=0$. Top: Low energy electronic bands. Bottom: Berry curvature
normalized to $2\pi{\cal C}$ when integrated over the Brillouin Zone, where
${\cal C}=\\{2,-1\\}$ is the Chern number of each band.

Figure S13: As in Fig.[S12] but for $\\{m,n\\}=\\{2,-3\\}$ and
$\\{\theta_{12},\theta_{23}\\}=\\{0.54^{\circ},-0.81^{\circ}\\}$. The Chern
numbers are ${\cal C}=\\{1,-2\\}$.

## References

  * Lopes dos Santos  _et al._ [2007] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Andrei and MacDonald [2020] E. Y. Andrei and A. H. MacDonald, Nature Materials 19, 1265 (2020).
  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Nature Physics 16, 725 (2020).
  * Chen and Lado [2021] G. Chen and J. L. Lado, Phys. Rev. Res. 3, 033276 (2021).
  * Zhang  _et al._ [2021] L. Zhang, F. Wu, S. Hou, Z. Zhang, Y.-H. Chou, K. Watanabe, T. Taniguchi, S. R. Forrest, and H. Deng, Nature 591, 61 (2021).
  * Khosravian and Lado [2022] M. Khosravian and J. L. Lado, Phys. Rev. Mater. 6, 094010 (2022).
  * Scammell and Scheurer [2023] H. D. Scammell and M. S. Scheurer, Phys. Rev. Lett. 130, 066001 (2023).
  * Li  _et al._ [2021] T. Li, S. Jiang, L. Li, Y. Zhang, K. Kang, J. Zhu, K. Watanabe, T. Taniguchi, D. Chowdhury, L. Fu, J. Shan, and K. F. Mak, Nature 597, 350 (2021).
  * F. Ferreira, S. Magorrian, V. Enaldiev, D. Ruiz-Tijerina, V. Fal’ko [2021] F. Ferreira, S. Magorrian, V. Enaldiev, D. Ruiz-Tijerina, V. Fal’ko, Appl. Phys. Lett. 118, 241602 (2021).
  * Vogl  _et al._ [2022] M. Vogl, S. Chaudhary, and G. A. Fiete, Journal of Physics: Condensed Matter 35, 095801 (2022).
  * Enaldiev  _et al._ [2022] V. V. Enaldiev, F. Ferreira, J. G. McHugh, and V. I. Fal’ko, npj 2D Materials and Applications 6, 74 (2022).
  * Chen  _et al._ [2022] G. Chen, M. Khosravian, J. L. Lado, and A. Ramires, 2D Materials 9, 024002 (2022).
  * Törmä  _et al._ [2022] P. Törmä, S. Peotta, and B. A. Bernevig, Nature Reviews Physics 4, 528 (2022).
  * Liu  _et al._ [2022] L. Liu, S. Zhang, Y. Chu, C. Shen, Y. Huang, Y. Yuan, J. Tian, J. Tang, Y. Ji, R. Yang, K. Watanabe, T. Taniguchi, D. Shi, J. Liu, W. Yang, and G. Zhang, Nature Communications 13, 3292 (2022).
  * Burg  _et al._ [2022] G. W. Burg, E. Khalaf, Y. Wang, K. Watanabe, T. Taniguchi, and E. Tutuc, Nature Materials 21, 884 (2022).
  * San-Jose  _et al._ [2012] P. San-Jose, J. González, and F. Guinea, Phys. Rev. Lett. 108, 216802 (2012).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Phys. Rev. Lett. 122, 106405 (2019).
  * Khalaf  _et al._ [2019] E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Phys. Rev. B 100, 085109 (2019).
  * Cao  _et al._ [2021] Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 595, 523 (2021).
  * Kim  _et al._ [2022] H. Kim, Y. Choi, C. Lewandowski, A. Thomson, Y. Zhang, R. Polski, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Nature 606, 494 (2022).
  * Park  _et al._ [2022] J. M. Park, Y. Cao, L. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature Mat. 21, 877 (2022).
  * Sheffer and Stern [2021] Y. Sheffer and A. Stern, Phys. Rev. B 104, L121405 (2021).
  * Eaton  _et al._ [2022] A. Eaton, Y. Li, H. A. Fertig, and B. Seradjeh, Phys. Rev. B 106, 045117 (2022).
  * Meng  _et al._ [2023] H. Meng, Z. Zhan, and S. Yuan, Phys. Rev. B 107, 035109 (2023).
  * Popov and Tarnopolsky [2023] F. K. Popov and G. Tarnopolsky,  (2023), arXiv:arXiv:2303.15505 [cond-mat.str-el] .
  * Morales-Durán  _et al._ [2023] N. Morales-Durán, J. Wang, G. R. Schleder, M. Angeli, Z. Zhu, E. Kaxiras, C. Repellin, and J. Cano,  (2023), arXiv:2304.06669 [cond-mat.str-el] .
  * Li  _et al._ [2022a] Y. Li, A. Eaton, H. A. Fertig, and B. Seradjeh, Phys. Rev. Lett. 128, 026404 (2022a).
  * Uri  _et al._ [2023] A. Uri, S. C. de la Barrera, M. T. Randeria, D. Rodan-Legrain, T. Devakul, P. J. D. Crowley, N. Paul, K. Watanabe, T. Taniguchi, R. Lifshitz, L. Fu, R. C. Ashoori, and P. Jarillo-Herrero 10.48550/arXiv.2302.00686 (2023).
  * Koshino [2015] M. Koshino, New Journal of Physics 17, 015014 (2015).
  * Amorim and Castro [2018] B. Amorim and E. V. Castro 10.48550/arXiv.1807.11909 (2018).
  * Cea  _et al._ [2019] T. Cea, N. R. Walet, and F. Guinea, Nano Letters 19, 8683 (2019).
  * Li  _et al._ [2019] X. Li, F. Wu, and A. H. MacDonald 10.48550/arXiv.1907.12338 (2019).
  * Mora  _et al._ [2019] C. Mora, N. Regnault, and B. A. Bernevig, Phys. Rev. Lett. 123, 026402 (2019).
  * Wu  _et al._ [2020] F. Wu, R.-X. Zhang, and S. Das Sarma, Phys. Rev. Res. 2, 022010 (2020).
  * Zhu  _et al._ [2020a] Z. Zhu, S. Carr, D. Massatt, M. Luskin, and E. Kaxiras, Phys. Rev. Lett. 125, 116404 (2020a).
  * Mao  _et al._ [2023] Y. Mao, D. Guerci, and C. Mora, Phys. Rev. B 107, 125423 (2023).
  * [38] See Supplementary Information.
  * Devakul  _et al._ [2023] T. Devakul, P. J. Ledwith, L.-Q. Xia, A. Uri, S. de la Barrera, P. Jarillo-Herrero, and L. Fu https://doi.org/10.48550/arXiv.2305.03031 (2023), arXiv:arXiv:2305.03031 [cond-mat.str-el] .
  * Nakatsuji  _et al._ [2023] N. Nakatsuji, T. Kawakami, and M. Koshino, arXiv.2305.13155 https://doi.org/10.48550/arXiv.2305.13155 (2023), arXiv:2305.13155 [cond-mat.mes-hall] .
  * Guerci  _et al._ [2023] D. Guerci, Y. Mao, and C. Mora, arXiv:2305.03702 https://doi.org/10.48550/arXiv.2305.03702 (2023), arXiv:arXiv:2305.03702 [cond-mat.mes-hall] .
  * 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).
  * Li  _et al._ [2023] Y. Li, Z. Zhan, X. Kuang, Y. Li, and S. Yuan, Computer Physics Communications 285, 108632 (2023).
  * Yu  _et al._ [2019] G. Yu, Z. Wu, Z. Zhan, M. I. Katsnelson, and S. Yuan, npj Comput. Mater. 5, 1 (2019).
  * Slater and Koster [1954] J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, _et al._ , Computer Physics Communications 271, 108171 (2022).
  * Brenner  _et al._ [2002] D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, Journal of Physics: Condensed Matter 14, 783 (2002).
  * Kolmogorov and Crespi [2005] A. N. Kolmogorov and V. H. Crespi, Phys. Rev. B 71, 235415 (2005).
  * J. M. Park, Y. Cao, L. Xia, S. Sun, K. Watanabe, T. Taniguchi and P. Jarillo-Herrero [2022] J. M. Park, Y. Cao, L. Xia, S. Sun, K. Watanabe, T. Taniguchi and P. Jarillo-Herrero, Nat. Mater. 21, 10.1038/s41563-022-01287-1 (2022).
  * Li  _et al._ [2022b] Y. Li, A. Eaton, H. A. Fertig, and B. Seradjeh, Phys. Rev. Lett. 128, 026404 (2022b).
  * A. Uri, S. C. de la Barrera, M. T. Randeria, D. Rodan-Legrain, T. Devakul, P. J. D. Crowley, N. Paul, K. Watanabe, T. Taniguchi, R. Lifshitz, L. Fu, R. C. Ashoori, and P. Jarillo-Herrero [2023] A. Uri, S. C. de la Barrera, M. T. Randeria, D. Rodan-Legrain, T. Devakul, P. J. D. Crowley, N. Paul, K. Watanabe, T. Taniguchi, R. Lifshitz, L. Fu, R. C. Ashoori, and P. Jarillo-Herrero, arXiv:2302.00686 10.48550/arXiv.2302.00686 (2023).
  * Zhu  _et al._ [2020b] Z. Zhu, S. Carr, D. Massatt, M. Luskin, and E. Kaxiras, Physical review letters 125, 116404 (2020b).
  * Goodwin  _et al._ [2019] Z. A. H. Goodwin, F. Corsetti, A. A. Mostofi, and J. Lischner, Phys. Rev. B 100, 235424 (2019).
  * Fedor K. Popov and Grigory Tarnopolsky [2023] Fedor K. Popov and Grigory Tarnopolsky, arXiv.2305.16385 10.48550/arXiv.2305.16385 (2023).


# Triple point fermions in ferroelectric GeTe

Juraj Krempaský1 Laurent Nicolaï2 Martin Gmitra3 Houke Chen4 Mauro Fanciulli5
Eduardo B. Guedes1 Marco Caputo1 Milan Radović1 V. V. Volobuev6,7 Ondřej Caha8
Gunther Springholz9 Jan Minár2 J. Hugo Dil1,10 1Photon Science Division, Paul
Scherrer Institut, CH-5232 Villigen, Switzerland
3New Technologies-Research Center University of West Bohemia, Plzeň, Czech
Republic
3Institute of Physics, P. J. Šafárik University in Košice, Park Angelinum 9,
040 01 Košice, Slovakia
4Department of Physics, Tsinghua University, Beijing 100084, China
5LPMS, Université de Cergy-Pontoise, 95031 Cergy-Pontoise, France
6International Research Centre MagTop, Institute of Physics, Polish Academy of
Sciences, Aleja Lotnikow 32/46, PL-02668 Warsaw, Poland
7National Technical University ”KhPI”, Kyrpychova Str. 2, 61002 Kharkiv,
Ukraine
8Department of Condensed Matter Physics, Masaryk University, Kotlářská 267/2,
61137 Brno, Czech Republic
9Institut für Halbleiter-und Festkörperphysik, Johannes Kepler Universität,
A-4040 Linz, Austria
10Institut de Physique, École Polytechnique Fédérale de Lausanne, CH-1015
Lausanne, Switzerland

###### Abstract

Ferroelectric $\alpha$-GeTe is unveiled to exhibit an intriguing multiple non-
trivial topology of the electronic band structure due to the existence of
triple-point and type-II Weyl fermions, which goes well beyond the giant
Rashba spin splitting controlled by external fields as previously reported.
Using spin- and angle-resolved photoemission spectroscopy combined with ab
initio density functional theory, the unique spin texture around the triple
point caused by the crossing of one spin degenerate and two spin-split bands
along the ferroelectric crystal axis is derived. This consistently reveals
spin winding numbers that are coupled with time reversal symmetry and Lorentz
invariance, which are found to be equal for both triple-point pairs in the
Brillouin zone. The rich manifold of effects opens up promising perspectives
for studying non-trivial phenomena and multi-component fermions in condensed
matter systems.

The exploration of topological materials commenced with the discovery of
topological insulators Kane and Mele (2005); Hasan and Kane (2010a); Hsieh
_et al._ (2009a); Nishide  _et al._ (2010); Hsieh  _et al._ (2009b); Pan  _et
al._ (2011); Miyamoto  _et al._ (2012); Landolt  _et al._ (2014); Hasan and
Kane (2010b); Bradlyn  _et al._ (2016), followed by topological semimetals Lv
_et al._ (2015); Soluyanov  _et al._ (2015); Armitage  _et al._ (2018); Xu
_et al._ (2015), topological superconductors Hao and Hu (2018) and nodal line
semimetals Fang  _et al._ (2016). Recently, topological semimetals with
triply-degenerate points were predicted in materials such as molybdenum
phosphide Lv  _et al._ (2017), tungsten carbide Ma  _et al._ (2018),
InAs1-xSbx alloys Winkler  _et al._ (2016), and strained HgTe Zaheer  _et al._
(2013). In contrast to an even number of band crossings, band structures with
triple band crossing are rather rare in condensed matter. Such triple-point
(TP) fermions can be found along the high symmetry axis in the Brillouin zone
(BZ) allowing for two and one dimensional double group representations Winkler
_et al._ (2016, 2019); Zhu  _et al._ (2016). This is also the case for
$\alpha$-GeTe(111) which was so far extensively studied in terms of the
topologically trivial Rashba-type band structure Di Sante  _et al._ (2013);
Krempaský  _et al._ (2016a, 2017).

The origin of these intriguing features in $\alpha$-GeTe lies in the
spontaneous structural phase transition of the inversion symmetric rock-salt
$\beta$-GeTe structure (FCC, space group #225) into a ferroelectric
rhombohedral structure with Ge-atoms displaced by $\approx$0.3 $\mathrm{\AA}$
and distorted by $\approx$2∘ (space group #160) Krempaský  _et al._ (2016a);
Kriegner  _et al._ (2019). This inversion symmetry breaking introduces a spin-
splitting except for the bands that line up with the high-symmetry
ferroelectric $\langle 111\rangle$ axis, which is the main symmetry axis of
the $C_{3v}$ point group. Along this axis TP fermions can form which are
classified as type-A or B, depending on the number of nodal lines that connect
a pair of TP.

Figure 1: (a) Bulk BZ and surface BZ of rhombohedral $\alpha$-GeTe with
$\langle 111\rangle$ ferroelectric axis along the $k_{\mathrm{z}}$ direction.
(b) BSF initial state calculations along $k_{\mathrm{x}}$ for sixteen selected
$k_{\mathrm{z}}$ values denoted in panel (a). (c) Band structure along
$k_{\mathrm{z}}$ with position of $\mathrm{TP}^{\sigma}_{-}$ and
$\mathrm{TP}^{\sigma}_{+}$ indicated by orange arrows.

In this Letter we provide theoretical and experimental evidence that
$\alpha$-GeTe possess type-B triple point fermions inside its valence band
states along the ferroelectric $\langle 111\rangle$ axis. To this date such
non-trivial topology in $\alpha$-GeTe failed to be noticed. The spin texture
around the TP is unravelled and the same spin winding number of $\pm 1$ is
found for both TP pairs in the 3D BZ. Furthermore, we identified additional
Dirac and type-II Weyl fermions that coexist with these three-component
fermions.

Epitaxial films of 200 nm of $\alpha$-GeTe for our (spin-) and angle-resolved
photoemission spectroscopy ((S)ARPES) studies were grown by molecular beam
epitaxy on BaF2(111) substrates and in situ transferred for high-resolution
ARPES at the UE112-PGM2a beamline of BESSY II. SARPES measurements were
performed at the Swiss Light Source of the Paul Scherrer Institut at the
COPHEE end station equipped with two orthogonal Mott polarimeters to measure
the three spatial components of spinors in arbitrary reciprocal space points
Hoesch  _et al._ (2002). We compared our experimental data with calculated
initial state band-maps presented in terms of Bloch spectral functions (BSF)
as implemented in the fully relativistic spin-polarised Korringa-Kohn-Rostoker
theory Ebert  _et al._ (2011) and with the full potential all-electron Wien2K
density functional code. Considering the relatively small energy scales
involved, this allowed us to check the results for consistency.

Figure 2: (a) Calculated band structure along the
$\mathrm{U^{\prime}}$-$\mathrm{Z}$-$\mathrm{U}$ mirror plane, with triple
point (TP), type-II Weyl fermion (WP II) and Dirac point (DP). (b) Formation
of TP pair ($\mathrm{TP}^{-}_{-}$, $\mathrm{TP}^{+}_{-}$) by spin-degenerate
band $\Lambda_{6}$ and two spin split bands $\Lambda_{4,5}$. (c,d,e) Close-up
of band structure along $k_{\mathrm{y}}$ at three different values of
$k_{\mathrm{z}}$ across the TP pair. Three types of touching points are
indicated with black, green, and magenta markers, black arrows indicate the
$k_{\mathrm{z}}$ dispersion of the latter. (f) Band structure around the TP
under 10 T B-field along the z-direction. The TP splits into a pair of Weyl
points (WP). (g) Band structure for rock-salt $\beta$-GeTe above ferroelectric
phase transition. The triple fermions merges into a single Dirac point (DP).
(inset) band structure along $k_{\mathrm{x}}$.

Fig. 1a shows the BZ with two pairs of triple points
$\mathrm{TP}^{\sigma}_{z}$, where the $\sigma=\pm$ index denotes the TP within
a single pair, and $z=\pm$ denotes the two TP pairs with opposite $k_{z}$
within the BZ. Panel (b) shows valence band (VB) dispersions for sixteen
selected $k_{z}$ values in the $\mathrm{Z}$-$\mathrm{A}$-$\Gamma$ plane. The
top-left image presents the Rashba-type spin splitting of the uppermost VB at
the Z-point. When moving $k_{z}$ towards $\Gamma$ these bands disperse and
finally meet with the lower spin-split bands at approximately $\frac{1}{3}$BZ
distance from the $\Gamma$-point (orange arrow). Another visualization of the
TP formation is in panel (c) plotting the dispersion along $\Gamma\mathrm{Z}$.
At first glance this crossing point is reminiscent of the topologically
protected touching point between electron and hole pockets in a double-Weyl
scenarioTakane  _et al._ (2019); Sanchez  _et al._ (2019); Rao  _et al._
(2019), but Fig. 2b-e unambiguously resolve a three-component fermion matching
the space group symmetry. Besides this, the bulk electronic structure in the
$\mathrm{U^{\prime}}$-$\mathrm{Z}$-$\mathrm{U}$ mirror plane for a $k_{z}$
right at the TP in Fig. 2a shows two additional band crossings with non-
trivial topology. First, there is a Dirac point (DP) near 2 eV binding energy.
Second, a type-II Weyl point (WP II) appears around -0.4 Å-1 as indicated in
the figure. The presence of this Weyl point influences the spin texture and
changes the spin winding number from $\pm$1 to $\pm$2, as seen in Fig. 5e,f
and discussed in more detail below.

Figure 3:  (top) High-resolution ARPES data for various photon energies around
the TP. (bottom) Corresponding curvature maps showing the TP formation at 30.5
eV photon energy inside the orange frame. The arrows indicate surface states
SS, $\mathrm{SS_{2}}$, surface resonances (SR) and bulk bands (B).

Following the TP-classification by Winkler $et~{}al.$ Winkler  _et al._
(2019), $\alpha$-GeTe should host a type-B TP with four nodal lines. Fig. 2
summarizes this scenario: panel (b) resolves the bands $\Lambda_{4,5,6}$
dispersing in $k_{\mathrm{z}}$, whereby the two spin-split bands
$\Lambda_{4}$, $\Lambda_{5}$ cross the doubly degenerate $\Lambda_{6}$,
forming a pair of triple fermion quasiparticles. In contrast to type-A TPs
with only one nodal line, the triple points within each pair
$\mathrm{TP}^{\sigma}_{-}$ and $\mathrm{TP}^{\sigma}_{+}$ are connected by
four nodal lines. One nodal line is lined up with the $C_{3v}$ axis and
therefore shows up at zero in-plane momentum only (purple markers in Fig.
2c-e). On the other hand the three remaining nodal lines are extending in
three vertical mirror planes, which in our case line up along the three
$\mathrm{Z}$-$\mathrm{U}$ directions (orange lines in Fig. 1a). Panels (c-e)
are close-ups of the bands structure across the single TP${}^{\pm}_{-}$ pair
for three $k_{z}$ values. For the intermediate $k_{z}$ between
$\mathrm{TP}^{-}_{-}$ and $\mathrm{TP}^{+}_{-}$ there are two crossing points
forming two of the four nodal lines, whereby the one at finite $k_{\parallel}$
(black and green markers) disperse with $k_{\mathrm{z}}$ (see animated Figure
S1 in Supplemental Material).

As a whole, the triple fermion forms an intermediate state between Dirac and
Weyl fermions. Actually, as indicated in Fig. 2f, $\alpha$-GeTe inside a 10 T
magnetic field (which is many orders of magnitude lower than the exchange
field in (GeMn)Te Krempaský  _et al._ (2016b)) the TP splits into a pair of
Weyl points. Furthermore, above the ferroelectric transition temperature
($\mathrm{T_{C}}$ $\approx$700 K) the TP transforms into a Dirac point as seen
in panel (g). This renders GeTe rather unique as all three phases can be
achieved under realistic conditions in a single material system.

In our earlier $\alpha$-GeTe ARPES studies we extensively described the
measured electronic structure in terms of surface states (SS), surface
resonances (SR), and bulk states (B) Krempaský  _et al._ (2016a, b, a, 2017);
Kremer  _et al._ (2020); Krempaský  _et al._ (2020). Such a decomposition also
applies to the ARPES data in Fig. 3 and allows to find the appropriate region
for measuring the TP and its spin texture. Below we demonstrate that inside
the $\mathrm{-0.2<k_{y}<0.2\,\AA^{-1}}$ momentum span the additional spin tag
in SARPES data enables us to unambiguously resolve the individual bulk bands
around a TP-pair.

Based on the $k_{\mathrm{z}}$-series of band maps in Fig. 1b, the ARPES data
in Fig. 3c, and the schematic $k_{\mathrm{z}}$ vs. photon energy diagram
(Supplemental Figure S2), the TP is predicted, and experimentally verified,
for photon energies around $\mathrm{h\nu}$ = 30, 58 and 89 eV. The latter two
values correspond to consecutive TP pairs at opposite parts of the BZ, thus
enabling us to experimentally examine their symmetry connection.

Figure 4:  (a) ARPES band maps along $\overline{\mathrm{M}\Gamma\mathrm{M}}$
measured at the COPHEE setup with $\mathrm{h\nu}$=58 and 89 eV. Orange bullets
indicate $\mathrm{TP}^{\sigma}_{-}$ and $\mathrm{TP}^{\sigma}_{+}$. MDC A and
MDC B indicate where data in Fig. 5 is obtained. (b) CEM measured at
$\mathrm{h\nu}$=89 eV and -0.4 eV binding energy at COPHEE. Marker A indicates
the $k_{\mathrm{x,y}}$ locus used for SARPES EDC, next to the bulk bands
denoted with dashed red-blue lines. (c,d) SARPES EDC cut A total intensity and
related $P_{x,y,z}$ spin polarizations with corresponding
$\mathrm{B_{1,2,3,4}}$ bulk band fitting indicated in red-blue and surface-
derived bands in gray. (e) 3D TP cartoon illustrating considered EDC cut A
with corresponding in-plane spin-textures, showing spin vectors projected
mainly along $P_{x}$ with a characteristic (+ - - +) pattern around the TP.
Figure 5: (a) Tabulated spin textures for $\mathrm{TP}^{\sigma}_{-}$ and
$\mathrm{TP}^{\sigma}_{+}$ measured along the MDC A (MDC B) cuts 250 meV above
(below) the TP, as denoted in Fig. 4a. (b,c) Corresponding SARPES intensities
and 3D spin polarization fitting for $\mathrm{TP}^{\sigma}_{-}$ and
$\mathrm{TP}^{\sigma}_{+}$ respectively. Red and blue curves indicate the
$\mathrm{B_{1,2,3,4}}$ peak fitting, surface-derived bands are in gray. (d)
Analogous spin textures measured by rotating the sample by 180∘, the frames in
green and magenta indicate two inequivalent directions denoted in Fig. 1a. (e)
BSF constant energy maps, black arrow indicates the Weyl-II fermion, dashed
line its energy dispersion and hybridization with TP bands. (f) Close ups of
spin-resolved CEMs at -1.55 and -1.6 eV, respectively, showing a change in
spin winding number from $\pm$1 to $\pm$2 around the WP II point inclusion
with conspicuous ”spin swirl” denoted by dashed circle.

A notable feature of the TP electronic structure is that the bands gradually
shift between $\mathrm{TP}^{-}_{-}$ and $\mathrm{TP}^{+}_{-}$ as marked by
black arrows in Fig. 2c-e. This spectral weight transfer across the TP is
manifested in a series of high-resolution ARPES band maps in Fig. 3c, and
their curvature plots Zhang  _et al._ (2011) in the bottom panels. In accord
with the theoretical predictions (Fig. 1b), the $k_{\mathrm{z}}$-dispersion of
the spectral weight is progressing toward the TP around 30.5 eV, where it
becomes equally distributed above and below the TP. Interestingly, the band
maps inside the orange frame for $\mathrm{h\nu}$=30.5 eV appear to mimic the
TP band mixing seen in Fig. 2c-e, however with much larger momentum and energy
scale. Because the bulk spectral intensity around the TP totally suppresses
the weakly contributing surface states $\mathrm{SS_{2}}$ (horizontal dashed
line in Fig. 3), the contribution of the latter can be neglected in our SARPES
analysis. We note that surface states are typically associated with nodal
points in topological metals and semimetals, but $\mathrm{SS_{2}}$ is not
interconnecting the TP pairs because it is clearly off-set in energy at zero
momentum, as indicated inside the orange frame. In fact, because both TP pairs
line up in the $\Gamma\mathrm{Z}$-direction, the associated surface states
cannot show up on the (111) surface.

Similar to the Rashba model, the spin polarization of bulk states around a TP
is expected to be tangentially momentum-locked in accord with time-reversal
symmetry (TRS) properties, as also shown in the calculated spin texture in
Fig.5e. Our SARPES data is visualized as spin-resolved energy distribution
(EDC) and momentum distribution curves (MDC) summarized in Figures 4 and 5,
respectively. For all measurements our analysis concentrates on the in-plane
spin polarisation $P_{x,y}$ of the four spin-split bulk bands denoted by
$B_{1,2,3,4}$. Their properties are determined using a 3D vectorial fit of the
measured spin expectation values with a well established fitting routine that
extracts their positions, widths, intensities, and 3D spin vector Meier  _et
al._ (2009). For the EDCs a Shirley background is subtracted from the total
intensity, whereas for the MDCs a constant background is used. To simplify
comparison between theory and experiment, we aligned the sample in such a way
that the out-of-plane spin polarization $P_{z}$ is oriented along the
$k_{\mathrm{z}}$-direction and $P_{x}$ ($P_{y}$) along the
$\overline{\Gamma\mathrm{K}}$ ($\overline{\Gamma\mathrm{M}}$) direction.
Within this experimental geometry a tangential spin winding is determined by
$P_{x}$.

We first test the spin texture in terms of TRS and Lorentz invariance (LI) for
the selected $k_{\mathrm{x,y}}$ point labeled A in Fig. 4b, with SARPES data
summarized in panels (c-d). The resulting spin texture is summarized in panel
(e) next to a model visualizing the EDC cut through the four $B_{1-4}$ bulk
bands around the TP. We note that the spin helicity dictated by TRS and LI are
preserved: the spins are nearly tangentially locked in momentum, whereby bands
$B_{1,4}$ and $B_{2,3}$ do not change the spin helicity. The spin texture
imposed by LI, i.e. symmetric spin in energy with regard to the TP, becomes
apparent when considering the (+ - - +) pattern in Fig. 4d; in contrast to the
Rashba-type (+ - + -) pattern seen in Supplemental Figure S3, measured at the
BZ boundary Z-point (see Fig. 1a).

This experimental result is further confirmed in Fig. 5. We measured each TP-
pair with two MDC cuts 250 meV above (A) and below (B) the TP, as depicted in
Fig. 4a. In addition, because there are two inequivalent directions (Fig. 1a
green and magenta lines), each TP-pair was measured along U’-A-U and after
rotation by 180∘along U-A-U’ (see SARPES data in Supplemental Figure S4). The
resulting in-plane spin textures shown in Fig. 5d are consistent with a spin
winding number of $\pm 1$. To best of our knowledge this is the first
experimental evidence of the spin texture of a triple fermion, featuring a
winding number of $\pm$1 in accord with both TRS and LI and theoretical
predictions. Here it should be noted that the inclusion of WP II induces an
additional spin swirl and changes the winding number to $\pm 2$ as illustrated
by the calculations in Fig. 5e. Furthermore, the two triple point pairs
TP${}^{\sigma}_{-}$ and TP${}^{\sigma}_{+}$ show the same spin helicity in
accordance with the absence of a $(x,y)$ mirror plane.

In conclusion, we demonstrated that $\alpha$-GeTe is a type-B triple fermion
topological system enabled by the ferroelectric non-centrosymmetric atomic
arrangement of Ge and Te atoms along the $\langle 111\rangle$ axis.
Furthermore, we measured for the first time the spin winding number for a
triple point system. Moreover, $\alpha$-GeTe is currently unique in the
possibility to manipulate topological properties by temperature variation
Krempaský  _et al._ (2016a), ferroelectric switching Krempaský  _et al._
(2018), and/or magnetic doping Krempaský  _et al._ (2016b).

M.F. acknowledges support from the SNF Project No. P2ELP2_181877 and G.S. from
the Project No. P30960-N27 of the Austrian Science Fund (FWF). J.M. and L.N.
would like to thank the CEDAMNF Project financed by the Ministry of Education,
Youth and Sports of Czech Republic, Project No. CZ.02.1.010.00.015_003/0000358
and the Czech Science Foundation (GACR), Project No. 2018725S. M.G.
acknowledges VEGA 1010520 and Internal Research Grant System VVGS-2019-1227,
V.V.V. the Foundation for Polish Science through the IRA Programme co-financed
by EU within SG OP, and finally O.C. acknowledges the financial support by the
European Regional Development Fund Project CEITEC Nano+ (No.
CZ.021.010.00.016_0130001728).

## References

  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * Hasan and Kane (2010a) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010a).
  * Hsieh  _et al._ (2009a) D. Hsieh, Y. Xia, L. Wray, D. Qian, A. Pal, J. H. Dil, J. Osterwalder, F. Meier, G. Bihlmayer, C. L. Kane, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Science 323, 919 (2009a).
  * Nishide  _et al._ (2010) A. Nishide, A. A. Taskin, Y. Takeichi, T. Okuda, A. Kakizaki, T. Hirahara, K. Nakatsuji, F. Komori, Y. Ando, and I. Matsuda, Phys. Rev. B 81, 041309(R) (2010).
  * Hsieh  _et al._ (2009b) D. Hsieh, Y. Xia, D. Qian, L. Wray, J. H. Dil, F. Meier, J. Osterwalder, L. Patthey, J. G. Checkelsky, N. P. Ong, A. V. Fedorov, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Nature 460, 1101 EP (2009b).
  * Pan  _et al._ (2011) Z.-H. Pan, E. Vescovo, A. V. Fedorov, D. Gardner, Y. S. Lee, S. Chu, G. D. Gu, and T. Valla, Phys. Rev. Lett. 106, 257004 (2011).
  * Miyamoto  _et al._ (2012) K. Miyamoto, A. Kimura, T. Okuda, H. Miyahara, K. Kuroda, H. Namatame, M. Taniguchi, S. V. Eremeev, T. V. Menshchikova, E. V. Chulkov, K. A. Kokh, and O. E. Tereshchenko, Phys. Rev. Lett. 109, 166802 (2012).
  * Landolt  _et al._ (2014) G. Landolt, S. Schreyeck, S. V. Eremeev, B. Slomski, S. Muff, J. Osterwalder, E. V. Chulkov, C. Gould, G. Karczewski, K. Brunner, H. Buhmann, L. W. Molenkamp, and J. H. Dil, Phys. Rev. Lett. 112, 057601 (2014).
  * Hasan and Kane (2010b) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010b).
  * Bradlyn  _et al._ (2016) B. Bradlyn, J. Cano, Z. Wang, M. G. Vergniory, C. Felser, R. J. Cava, and B. A. Bernevig, Science 353 (2016).
  * Lv  _et al._ (2015) B. Q. Lv, H. M. Weng, B. B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. C. Huang, L. X. Zhao, G. F. Chen, Z. Fang, X. Dai, T. Qian, and H. Ding, Phys. Rev. X 5, 031013 (2015).
  * Soluyanov  _et al._ (2015) A. A. Soluyanov, D. Gresch, Z. Wang, Q. S. Wu, M. Troyer, X. Dai, and B. A. Bernevig, Nature 527, 495 (2015).
  * Armitage  _et al._ (2018) N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Xu  _et al._ (2015) S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, S.-M. Huang, H. Zheng, J. Ma, D. S. Sanchez, B. Wang, A. Bansil, F. Chou, P. P. Shibayev, H. Lin, S. Jia, and M. Z. Hasan, Science 349, 613 (2015).
  * Hao and Hu (2018) N. Hao and J. Hu, National Science Review 6, 213 (2018).
  * Fang  _et al._ (2016) C. Fang, H. Weng, X. Dai, and Z. Fang, Chinese Physics B 25, 117106 (2016).
  * Lv  _et al._ (2017) B. Q. Lv, Z.-L. Feng, Q.-N. Xu, X. Gao, J.-Z. Ma, L.-Y. Kong, P. Richard, Y.-B. Huang, V. N. Strocov, C. Fang, H.-M. Weng, Y.-G. Shi, T. Qian, and H. Ding, Nature 546, 627 (2017).
  * Ma  _et al._ (2018) J.-Z. Ma, J.-B. He, Y.-F. Xu, B. Q. Lv, D. Chen, W.-L. Zhu, S. Zhang, L.-Y. Kong, X. Gao, L.-Y. Rong, Y.-B. Huang, P. Richard, C.-Y. Xi, E. S. Choi, Y. Shao, Y.-L. Wang, H.-J. Gao, X. Dai, C. Fang, H.-M. Weng, G.-F. Chen, T. Qian, and H. Ding, Nature Physics 14, 349 (2018).
  * Winkler  _et al._ (2016) :w G. W. Winkler, Q. S. Wu, M. Troyer, P. Krogstrup, and A. A. Soluyanov, Phys. Rev. Lett. 117, 076403 (2016).
  * Zaheer  _et al._ (2013) S. Zaheer, S. M. Young, D. Cellucci, J. C. Y. Teo, C. L. Kane, E. J. Mele, and A. M. Rappe, Phys. Rev. B 87, 045202 (2013).
  * Winkler  _et al._ (2019) G. W. Winkler, S. Singh, and A. A. Soluyanov, Chinese Physics B 28, 077303 (2019).
  * Zhu  _et al._ (2016) Z. Zhu, G. W. Winkler, Q. S. Wu, J. Li, and A. A. Soluyanov, Phys. Rev. X 6, 031003 (2016).
  * Di Sante  _et al._ (2013) D. Di Sante, P. Barone, R. Bertacco, and S. Picozzi, Advanced Materials 25, 509 (2013).
  * Krempaský  _et al._ (2016a) J. Krempaský, H. Volfová, S. Muff, N. Pilet, G. Landolt, M. Radović, M. Shi, D. Kriegner, V. Holý, J. Braun, H. Ebert, F. Bisti, V. A. Rogalev, V. N. Strocov, G. Springholz, J. Minár, and J. H. Dil, Phys. Rev. B 94, 205111 (2016a).
  * Krempaský  _et al._ (2017) J. Krempaský, M. Fanciulli, N. Pilet, J. Minár, W. Khan, M. Muntwiler, F. Bertran, S. Muff, A. Weber, V. Strocov, V. Volobuiev, G. Springholz, and J. Dil, Journal of Physics and Chemistry of Solids (2017).
  * Kriegner  _et al._ (2019) D. Kriegner, G. Springholz, C. Richter, N. Pilet, E. M uller, M. Capron, H. Berger, V. Holý, J. H. Dil, and J. Krempaský, Crystals 9, 335 (2019).
  * Hoesch  _et al._ (2002) M. Hoesch, T. Greber, V. N. Petrov, M. Muntwiler, M. Hengsberger, W. Auwärter, and J. Osterwalder, J. Electron Spectrosc. Relat. Phenom. 124, 263 (2002).
  * Ebert  _et al._ (2011) H. Ebert, D. Ködderitzsch, and J. Minár, Rep. Prog. Phys. 74 (2011).
  * Takane  _et al._ (2019) D. Takane, Z. Wang, S. Souma, K. Nakayama, T. Nakamura, H. Oinuma, Y. Nakata, H. Iwasawa, C. Cacho, T. Kim, K. Horiba, H. Kumigashira, T. Takahashi, Y. Ando, and T. Sato, Phys. Rev. Lett. 122, 076402 (2019).
  * Sanchez  _et al._ (2019) D. S. Sanchez, I. Belopolski, T. A. Cochran, X. Xu, J.-X. Yin, G. Chang, W. Xie, K. Manna, V. S uss, C.-Y. Huang, N. Alidoust, D. Multer, S. S. Zhang, N. Shumiya, X. Wang, G.-Q. Wang, T.-R. Chang, C. Felser, S.-Y. Xu, S. Jia, H. Lin, and M. Z. Hasan, Nature 567, 500 (2019).
  * Rao  _et al._ (2019) Z. Rao, H. Li, T. Zhang, S. Tian, C. Li, B. Fu, C. Tang, L. Wang, Z. Li, W. Fan, J. Li, Y. Huang, Z. Liu, Y. Long, C. Fang, H. Weng, Y. Shi, H. Lei, Y. Sun, T. Qian, and H. Ding, Nature 567, 496 (2019).
  * Krempaský  _et al._ (2016b) J. Krempaský, S. Muff, F. Bisti, M. Fanciulli, H. Volfová, A. P. Weber, N. Pilet, P. Warnicke, F. Bertran, H. Ebert, J. Braun, J. Minár, G. Springholz, J. Dil, and V. Strocov, Nat. Commun. 7, 13071 (2016b).
  * Kremer  _et al._ (2020) G. Kremer, T. Jaouen, B. Salzmann, L. Nicolaï, M. Rumo, C. W. Nicholson, B. Hildebrand, J. H. Dil, J. Minár, G. Springholz, J. Krempaský, and C. Monney, Phys. Rev. Research 2, 033115 (2020).
  * Krempaský  _et al._ (2020) J. Krempaský, M. Fanciulli, L. Nicolaï, J. Minár, H. Volfová, O. Caha, V. V. Volobuev, J. Sánchez-Barriga, M. Gmitra, K. Yaji, K. Kuroda, S. Shin, F. Komori, G. Springholz, and J. H. Dil, Phys. Rev. Research 2, 013107 (2020).
  * Zhang  _et al._ (2011) P. Zhang, P. Richard, T. Qian, Y.-M. Xu, X. Dai, and H. Ding, Review of Scientific Instruments 82, 043712 (2011).
  * Meier  _et al._ (2009) F. Meier, J. H. Dil, and J. Osterwalder, New J. Phys. 11, 125008 (2009).
  * Krempaský  _et al._ (2018) J. Krempaský, S. Muff, J. Minár, N. Pilet, M. Fanciulli, A. P. Weber, E. B. Guedes, M. Caputo, E. Müller, V. V. Volobuev, M. Gmitra, C. A. F. Vaz, V. Scagnoli, G. Springholz, and J. H. Dil, Physical Review X 8, 021067 (2018).