††thanks: These authors contributed equally.††thanks: These authors
contributed equally.††thanks: These authors contributed equally.

# Magic momenta and three dimensional Landau levels from a three dimensional
graphite moiré superlattice

Xin Lu  School of Physical Science and Technology, ShanghaiTech University,
Shanghai 201210, China  Bo Xie  School of Physical Science and Technology,
ShanghaiTech University, Shanghai 201210, China  Yue Yang  School of Physical
Science and Technology, ShanghaiTech University, Shanghai 201210, China  Xiao
Kong  Institute of Technology for Carbon Neutrality, Shenzhen Institute of
Advanced Technology, Chinese Academy of Sciences, Shenzhen, China  Shanghai
Institute of Microsystem and Information Technology, Chinese Academy of
Sciences, Shanghai, China  Jun Li  School of Physical Science and Technology,
ShanghaiTech University, Shanghai 201210, China  ShanghaiTech Laboratory for
Topological Physics, ShanghaiTech University, Shanghai 201210, China  Feng
Ding  Institute of Technology for Carbon Neutrality, Shenzhen Institute of
Advanced Technology, Chinese Academy of Sciences, Shenzhen, China  Shanghai
Institute of Microsystem and Information Technology, Chinese Academy of
Sciences, Shanghai, China  Zhu-Jun Wang<EMAIL_ADDRESS>School
of Physical Science and Technology, ShanghaiTech University, Shanghai 201210,
China  Jianpeng Liu<EMAIL_ADDRESS>School of Physical Science and
Technology, ShanghaiTech University, Shanghai 201210, China  ShanghaiTech
Laboratory for Topological Physics, ShanghaiTech University, Shanghai 201210,
China  Liaoning Academy of Materials, Shenyang 110167, China

###### Abstract

Twisted bilayer graphene (TBG) and other quasi-two-dimensional moiré
superlattices have attracted significant attention due to the emergence of
various correlated and topological states associated with the flat bands in
these systems. In this work, we theoretically explore the physical properties
of a new type of three dimensional graphite moiré superlattice, the bulk
alternating twisted graphite (ATG) system with homogeneous twist angle, which
is grown by in situ chemical vapor decomposition method. Compared to TBG, the
bulk ATG system is bestowed with an additional wavevector degrees of freedom
due to the extra dimensionality. As a result, we find that when the twist
angle of bulk ATG is smaller than twice of the magic angle of TBG, there
always exist “magic momenta” at which the in-plane Fermi velocities of the
moiré bands vanish. Moreover, topologically distinct flat bands of TBG at
different magic angles can even co-exist at different out-of-plane wavevectors
in a single bulk ATG system. Most saliently, when the twist angle is
relatively large, exactly dispersionless three dimensional zeroth Landau level
would emerge in the bulk ATG, which may give rise to robust three dimensional
quantum Hall effects over a large range of twist angles.

Twisted bilayer graphene (TBG) [1, 2, 3] has aroused significant interest in
recent years. In TBG, the mismatch between two honeycomb lattices forms the
moiré pattern with periodic modulations of $AA$ and $AB$ regions in real space
[4, 5], which enables both intra- and inter-sublattice interlayer couplings.
In particular, when the twist angle is around a discrete set of values known
as “magic angles” [1], the lowest two bands per spin per valley become ultra-
flat, and even turn out to be exactly flat in the chiral limit [6]. These flat
bands are further found to be topologically nontrivial with Landau-level like
wavefunctions [6, 7, 8, 9, 10, 11, 12, 13, 14]. The unprecedented flatness
combined with nontrivial band topology yields plentiful phenomena in magic-
angle TBG such as superconductivity, correlated insulators, quantum anomalous
Hall effect, orbital magnetic states, and so on [15, 16, 17, 18, 19, 20, 21,
22, 23, 18, 24, 25, 26, 27, 28, 19, 20, 21, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38], which have stimulated tremendous theoretical research [39, 40,
41, 11, 42, 12, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57,
58].

It is well appreciated that tuning a TBG sample to the vicinity of the first
magic angle ($\theta^{(1)}_{m}\\!\approx\\!1.05^{\circ}$) is extremely
challenging, especially if the device were fabricated based on the traditional
transferring and stacking technique [59, 60]. Nevertheless, the stringent
magic-angle condition for the emergence of flat bands in TBG can be somehow
released by stacking more twisted graphene layers in the out-of-plane
direction, forming twisted multilayer graphene (TMG) [61, 62, 63, 64, 65, 66,
67, 68, 69] and alternating twisted multilayer graphene (ATMG) [70, 71, 72,
73, 74] systems, where topological flat bands are robustly present over a
finite range of twist angles or at larger values of magic angles. Various
exotic states such as unconventional superconductivity [75], quantum anomalous
Hall effect [76], and generalized Wigner crystals [77] etc., have also been
observed in these TMG and ATMG systems. From the perspective of fundamental
science, it is then natural to ask what kind of new states and new phenomena
would emerge when an infinite number layers of twisted graphene are stacked on
top of each other, forming a three dimensional graphite moiré superlattice.

Recently, an origami–kirigami approach by chemical vapor deposition (CVD)
growth has been developed for the fabrication of graphene spiral [78]. Using
such a technique, double-helix structure of hundreds of layers of the
intertwined graphene with a uniform twisted interlayer stacking angle has been
successfully achieved (see Supplementary Video 1). The structural model of
such a graphene spiral, as depicted in Fig. 1(a), showcases this dual-helical
configuration. The top view of the dual-helical graphene model ((Fig. 1(b)))
indicates a periodic superlattice arrangement. Away from the central
dislocation line, the system exhibits an alternating twisted pattern between
each of the two adjacent graphene layers, wherein intra-plane moiré
superlattice vectors ($\mathbf{t}_{1},\mathbf{t}_{2}$ in Fig. 1(b)) can be
identified. In principle, it is feasible to achieve an infinite growth of
numerous graphene layers along the vertical direction defined by the helical
dislocation, which for the first time realizes a genuine three-dimensional
moiré graphite superlattice with alternating interlayer twist angle. In Fig.
1(b1), the Transmission Electron Microscopy (TEM) image of the twisted
graphene spiral sample is presented. Employing Mask-filtering and Inverse Fast
Fourier transform processing in Digital Micrograph, the TEM image distinctly
demonstrate a moiré periodic pattern characterized by an average superlattice
constant of $L_{s}=4.9\,$nm (as illustrated in the inset of Fig. 1(b1)). The
corresponding twisted angle $\theta\\!\approx\\!2.9\,^{\circ}$ is obtained
through the relationship $L_{s}=a/(2\sin{\theta/2})$ ($a=2.46$ Å is graphene’s
lattice constant). As depicted in Fig. 1(b2), the electron diffraction
patterns from two sets of twisted graphene layers are clearly distinguishable,
yielding a measured twist angle of approximately $3\,^{\circ}$. This result is
in concordance with the deduction based on the real-space moiré superlattice
constant. The cross-sectional scanning transmission electron microscopy (STEM)
image offers a contrasting feature of the double-helix graphene spiral. As
demonstrated in Fig. 1(c1), the helical dislocation nucleus is centrally
positioned, giving rise to the growth of two sets of alternatingly twisted
graphene layers following the helical dislocation path. Based on the
aforementioned details, the double-helix structural model for numerous
graphene layers was constructed as schematically shown in Fig. 1(c). To
further facilitate a comprehensive understanding of the cross-sectional
intricacies of the helical dislocation, simulations were conducted based on
the model depicted in Fig. 1(b) and Fig. 1(c), as presented in Fig. 1(d1).
Figs. 1(c2) and (d2) respectively illustrate the line profiles within the
blue, green, and red rectangle regions in Figs. 1(c1) and (d1). It is evident
that the profiles from the two sets exhibit quantitative consistency. In
comparison to the regions on both sides (blue and red rectangle), the helical
core region (green rectangle) exhibits an approximate vertical displacement of
0.18 nm (note that the horizontal axes in Figs. 1(c2) and (d2) denote the $z$
coordinate), equivalent to a half of graphene’s interlayer distance. A
noteworthy observation is that this screw dislocation itself hosts bound
states, elaborated further in the Supplementary Information [79]. The
successful realization of this macroscopically extensive bulk alternating
twisted graphite (ATG) system effectively extends the realm of intriguing
moiré graphene superlattices physics from two dimensions to three dimensions.

Figure 1: Structural characterizations of twisted graphene spiral (bulk
alternating twisted graphite). (a) The constructed structural model of a
graphene spiral with alternating interlayer twist angle. (b) The top view of
the helical graphene which reveals a periodic moiré pattern. Away from the
central dislocation line, the system exhibits bulk alternating twisted pattern
between every two adjacent graphene layers, wherein in-plane moiré
superlattice vectors $\mathbf{t}_{1}$, $\mathbf{t}_{2}$ are discernible. b1
presents a high-resolution TEM image of the twisted spiral graphene sample,
from which the moiré period is measured, corresponding to an interlayer
twisted angle of $2.9\,^{\circ}$. The angular measurement obtained from
electron diffraction in b2 corroborates with the angle depicted in b1. (c)
Schematic illustration of the cross-section of dual-helical twisted graphene
spiral. The annotations $\mathbf{t}_{B1}$ and $\mathbf{t}_{B2}$ indicate the
Burgers vectors of the two inequivalent sets of spiral graphene layers twisted
with respect each other by a finite angle. (c1) The STEM image of the twisted
graphene spiral sample. (d1) The simulated STEM image of the structural model
in (c). Comparative contrast line profile analysis (c2 and d2) reveals that
the simulated results are quantitatively consistent with the experimentally
measured ones.

In this work, we theoretically study the electronic structures, topological
properties, and Landau levels of the bulk ATG system. The additional
dimensionality endows bulk ATG an additional wavevector degrees of freedom to
explore flat-band physics of magic-angle TBG. In particular, we find that when
the twist angle $\theta\leq 2\theta^{(1)}_{m}\approx 2.1\,^{\circ}$, one can
always find a “magic quasi momentum” $\hbar k_{z}^{*}$ at which the 2D band
structures and wavefunctions within the $k_{x}$-$k_{y}$ plane are exactly the
same as that of TBG at the first magic angle. Moreover, when the twist angle
$\theta$ is further smaller than twice of the $n$th magic angle
$\theta_{m}^{(n)}$, i.e., $\theta\leq 2\theta^{(n)}_{m}$ (with $n=2,3,...$),
then for each $\theta$, there are at least $n$ magic momenta
$\\{k_{z}^{(s),*},s=1,...,n\\}$ at which the band structures and wavefunctions
within the $k_{x}$-$k_{y}$ planes are exactly the same as those of TBG at the
$s$th magic angle $\theta_{m}^{(s)}$. In other words, when
$\theta\leq\theta_{m}^{(2)}\approx 1\,^{\circ}$, multiple flat bands with
distinct topological properties would co-exist in the same bulk ATG system,
which may give rise to unprecedented correlated and topological phases of
matter. On the other hand, for larger twist angles
$\theta\\!>\\!2\theta^{(1)}_{m}\approx 2.1^{\circ}$, Dirac-cone type band
structures with finite Fermi velocities persist for every $k_{z}$ in the
Brillouin zone of bulk ATG. The Dirac points spectacularly align on a chain
along the $k_{z}$ direction, forming an exactly straight Weyl nodal chain for
each spin and valley species. Such quasi-2D Fermi surface would give rise to
weakly dispersive three dimensional Landau levels that may lead to 3D quantum
Hall effects, especially for twist angle $\theta\gtrapprox 3\,^{\circ}$. Most
saliently, given that the zeroth Landau level of graphene is exactly pinned to
the Dirac point, and that the Dirac points in bulk ATG are energetically
aligned at different $k_{z}$, there exists a branch of exactly dispersionless
zeroth 3D Landau band emerging from such 3D bulk ATG system under finite
vertical magnetic fields. Unusual correlated states may emerge at partial
fillings of such massively degenerate 3D zeroth Landau levels.

#### Magic momenta –

Neglecting the spiral dislocation at the core, we can consider bulk ATG as an
infinite aligned stacking of TBG in the out-of-plane direction with lattice
constant $2d_{0}\\!=\\!6.7\,$Å. The coupling between two nearest neighbor TBG
motifs is just the interlayer moiré potential. Such picture leads
straightforward to a low-energy effective Hamiltonian for bulk ATG [79] based
on the continuum model [1, 2, 80] of TBG with twist angle $\theta$:

$\displaystyle\textrm{H}^{\mu}_{\theta}=\hbar
v_{F}\left[\begin{array}[]{cc}-(\mathbf{k}-\mathbf{K}_{1}^{\mu})\cdot\bm{\sigma}^{\mu}&T(1+e^{i\phi})\\\
T^{\dagger}(1+e^{-i\phi})&-(\mathbf{k}-\mathbf{K}_{2}^{\mu})\cdot\bm{\sigma}^{\mu}\end{array}\right]$
(3)

where the Fermi velocity is $\hbar v_{F}/a\\!=\\!2.1354\,$eV [81] with
graphene’s lattice constant $a=2.46\,$Å. Neglecting the intervalley coupling,
the diagonal blocks represent the $\mathbf{k}\cdot\textbf{p}$ Hamiltonian of
the two layers of graphene in valley $\mu=\pm 1$ near the Dirac points
$\mathbf{K}_{1/2}^{\mu}$ with the Pauli matrices
$\bm{\sigma}^{\mu}=(\mu\sigma_{x},\sigma_{y})$ defined in the sublattice
space, and $\mathbf{k}=(k_{x},k_{y})$. The off-diagonal term stands for the
moiré potential between two adjacent graphene layers. The phase factor $e^{\pm
i\phi}=e^{\pm i2k_{z}d_{0}}$ is the only new contribution arising from the
nearest neighbor interlayer couplings along the opposite out-of-plane
directions in such 3D bulk system. Here, it is naturally expected that the
distance between two adjacent graphene layers is fixed, denoted as $d_{0}$,
since there is no space for out-of-plane atomic corrugations in such 3D bulk
moiré superlattice. This postulation is verified by direct molecular dynamics
simulations of bulk ATG, which gives an equilibrium interlayer distance
$d_{0}=3.35\,$Å [79].

The moiré potential $T$ are characterized by intra- and inter-sublattice
interlayer tunnelling amplitudes $u_{0}\\!=\\!u^{\prime}_{0}\\!=\\!0.103\,$eV.
They are equal to each other due to the absence of out-of-plane corrugations
in the bulk ATG system [79]. Up to an overall bandwidth $\hbar
v_{F}k_{\theta}$ (with $k_{\theta}=4\pi/(3L_{s})$), when
$u_{0}=u_{0}^{\prime}$ the continuum model of TBG is fully characterized by a
dimensionless ratio $\alpha_{0}(\theta)=u^{\prime}_{0}/(\hbar
v_{F}k_{\theta})$. Similarly, since $k_{z}$ is a good quantum number and the
bulk ATG can be seen as a series of decoupled TBG indexed by $k_{z}$, the non-
interacting physics of each of them is totally governed by

$\alpha(k_{z},\theta)=\alpha_{0}(\theta)|1+e^{2ik_{z}d_{0}}|\;,$ (4)

where $k_{z}d_{0}\in[0,\pi]$. Only the modulus matters because the phase can
be removed by a gauge transformation. Therefore, given a twist angle $\theta$,
one can explore a range of $k_{z}$-dependent effective moiré potential, which
varies from zero ($k_{z}=\pi/2d_{0}$) to doubled moiré potential strength
($k_{z}=0$) of the corresponding isolated TBG with twist angle $\theta$. In
the small angle approximation, as long as $\theta/2$ is smaller than $s$th
magic angle $\theta_{m}^{(s)}$ of TBG, one can always find the corresponding
$s$th magic momentum $k_{z}^{(s),*}$ such that
$\alpha(k_{z}^{(s),*},\theta)=\alpha_{0}(\theta_{m}^{(s)})$, as shown in Fig.
2(a). For example, to achieve the flat bands at the first magic angle in TBG,
it only suffices to twist the bulk ATG of angle
$\theta\\!<\\!2\theta_{m}^{(1)}\\!\approx\\!2.1\,^{\circ}$. Furthermore, the
bulk ATG can simultaneously host flat bands of TBG at several magic angles if
$\theta\\!<\\!2\theta_{m}^{(n)}$ with $n\\!\geq\\!2$. Since the $s$th magic
angle of TBG approximately equals one $s$th of the first magic angle, the
total number of flat bands (denoted by $n$) found at different out-of-plane
wavevectors in bulk ATG at twist angle $\theta$ verifies a thumb rule:
$n=[2\theta_{m}^{(1)}/\theta]$ ($[\textrm{ }]$ denotes the round down
operation). Accordingly, we can define $n$ magic momenta
$\\{k_{z}^{(s),*},s=1,...,n\\}$ for bulk ATG at fixed twist angle $\theta$.
For example, if $\theta=0.45\,^{\circ}<2\times\theta_{m}^{(4)}$, the flat
bands of TBG from first to fourth magic angle emerge concomitantly in the bulk
ATG at four magic momenta, as shown in Fig. 2(b). In Fig. 2(c), we present the
moiré band structures of bulk ATG with $\theta=0.8\,^{\circ}$ at several
different $k_{z}$. In such a case,
$2\theta_{m}^{(3)}\\!<\\!\theta\\!<\\!2\theta^{(2)}_{m}$, topological flat
bands at the first and second magic angles of TBG co-exist at different
magnetic momenta $k_{z}^{(1),*}$ and $k_{z}^{(2),*}$, as plotted by cyan and
magenta lines in Fig. 2(c).

#### Band Topology –

Since bulk ATG integrates a continuum of $k_{z}$-indexed TBG into a single
bulk system, it is intriguing to explore the topological properties by varying
both twist angle $\theta$ and $k_{z}$. By introducing a small sublattice mass
term ($\sim 0.1\,$meV) into the Hamiltonian Eq. (3), all the degeneracy points
are gapped out such that every band can be associated with a Chern number. At
each $k_{z}$ and $\theta$, we calculate the Chern number of the valence flat
band from valley $K$, and the calculated topological phase diagram is shown in
Fig. 2(d). First, the Chern number can change only if the twist angle of bulk
ATG $\theta$ is smaller than $\theta_{c,0}=2.3^{\circ}$, which corresponds to
the apex point of the rightmost phase boundary. From right to left across
$\theta_{c,0}$, the Chern number at $k_{z}=0$ switches from 1 to 0.
Nevertheless, the Chern number shifts back to unity if one continues to cross
the apex point of the second rightmost phase boundary, indicating another
critical angle $\theta_{c,1}=2.15\,^{\circ}$. If one continues to decrease
$\theta$, a series of such phase boundaries are passed, resulting in a rich
topological phase diagram in the ($k_{z}$,$\theta$) space as shown in Fig.
2(d).

If the twist angle is fixed between $\theta_{c,0}$ and $\theta_{c,1}$, e.g.,
$\theta=2.2^{\circ}$, the Brillouin zone can be cut into two topologically
different parts: the Chern number is zero if $k_{z}\in[-k_{z,0},k_{z,0}]$ and
is one otherwise. A close look on the phase boundary shows that such
topological transition is triggered by an accidental crossing at $(0,0,\pm
k_{z,0})$ between the valence flat band and the remote band. Therefore, the
system behaves as a Weyl semimetal with a pair of Weyl nodes with opposite
topological charges at $(0,0,\pm k_{z,0})$ in each valley and spin if spatial
inversion symmetry is broken by a slight sublattice mass term. Note that in
the absence of small mass term, the band touching point is three-fold
degenerate forming a three-fold multifermion because the remote band is always
twice degenerate dictated by symmetry [82, 79]. Furthermore, with smaller
twist angle, bulk ATG (with small sublattice mass term) can host eight
($\theta=1.6^{\circ}$), twelve ($\theta=1.2^{\circ}$), and more generally,
$4n$ pairs of Weyl nodes with $n=1,2,3,4,\dots$. Most spectacularly, the Weyl
nodes in bulk ATG perfectly align on the $k_{z}$-axis. The topological charge
of Weyl node is not limited to unity. For example, one pair of Weyl nodes in
bulk ATG of $\theta=0.85^{\circ}$ have topological charges equal to $\pm 2$
since the Chern number of the $k_{z}$-indexed 2D flat band shifts by 2 across
the Weyl nodes.

Figure 2: (a) Evolution of magic momenta $k_{z}^{(m),*}$ as a function of
physical twist angle $\theta$. Higher $m$th magic momenta also exist but not
shown here for clarity. The green dots are the results from tight-binding
calculations [79] (b) Evolution of Fermi velocity near Moiré $K$ point as a
function of $k_{z}$ for $\theta=0.45^{\circ}$. The red dots indicate the four
magic momenta $k_{z}^{*}$, at which the Fermi velocity vanishes exactly and
the flat band emerges. (c) Evolution of bulk ATG’s 2D band structures in the
plane $(k_{x},k_{y})$ from $k_{z}d_{0}=0$ to $\pi/2$. The twist angle of bulk
ATG is chosen to be $0.8^{\circ}$, at which bulk ATG can host two magic
momenta $k_{z}^{(1),*}$ and $k_{z}^{(2),*}$, as shown in (a). (d) Topological
phase diagram of the flat bands of $k_{z}$-indexed TBG in bulk ATG at
different twist angle. The color indicates the Chern number of the valence
flat band in valley $K$. (e) Quasi-1D Fermi surface of ATG of
$\theta=7^{\circ}$ for $E_{F}=0.1$ eV.

#### 3D quantum Hall effect

–

The unusual band dispersions and nontrivial topological properties of bulk ATG
imply its unconventional responses to electromagnetic fields. It follows from
Eq. (3) that $k_{z}$ plays a special role compared to the other two in-plane
momenta, leading to a strongly anisotropic, quasi-1D Fermi surface as shown in
Fig. 2(e), which is desirable to realize 3D quantum Hall effect [83, 84]. To
this end, a key requirement needs to be fulfilled: the twist angle of the bulk
ATG should be sufficiently larger than $2\theta_{m}^{(1)}$, otherwise at some
$k_{z}$ the lowest bands would flatten, which suppresses Landau level spacing
thus disfavors quantum Hall effect. In the case of large twist angles, the
moiré potential can be treated by perturbation theory, and the bulk ATG around
each Dirac point can be described by a Dirac Hamiltonian with renormalized
Fermi velocity

$H_{\rm{eff}}=\hbar
v_{F}[1-9\alpha^{2}(k_{z},\theta)](k_{x}\sigma_{x}+k_{y}\sigma_{y})\;,$ (5)

which is valid if $\theta\\!\gtrapprox\\!6^{\circ}$ based on the criterion
$9\alpha^{2}(k_{z},\theta)<0.1$. In practice, 3D quantum Hall effect is
present at even smaller angles (numerically checked at least for
$\theta\\!\gtrapprox\\!4\,^{\circ}$ at $E_{F}\\!\sim\\!20\,$meV where the
trigonal warping effect is still weak). Based on the Hamiltonian Eq. (5), at
fixed $k_{z}$, the Landau level’s energy should be exactly the same as
graphene but with renormalized Landau level spacing, following $\sqrt{B}$
behavior. Remarkably, for every $k_{z}$ there exists the zeroth Landau level
pinned at the Dirac point, which remains robust even in the presence of slight
inhomogeneity and disorder by virtue of Atiyah–Singer index theorem [85]. This
gives rise to exactly dispersionless 3D zeroth Landau level in such bulk ATG
system.

Figure 3: (a) Overall gap between the zeroth and the first 3D Landau levels
under magnetic field $B=3$-$15\,$T for twist angle
$\theta=3.0$-$9.5\,^{\circ}$. The color coding indicates the amplitude of gap.
(b) Overall gap between the first and the second 3D Landau levels in the same
parameter space as (a). The gap is closed for $\theta\lessapprox 7\,^{\circ}$.
(c) Colormap of the filling factor of 3D Landau levels in bulk ATG as a
function of 3D electron density and magnetic field. The white dashed lines
indicates the fillings with fully occupied 3D Landau levels. (d) Hofstadter
butterfly of $k_{z}$-indexed TBG in bulk ATG for $k_{z}=0$ and
$\theta=7^{\circ}$ for one spin sector. Green dashed lines are the calculated
Laudau level dispersions for $n=1,2,3$ using the approximated Hamiltonian Eq.
(5). The Chern number of gaps are given. (e) Landau level dispersion along the
$k_{z}$-dispersion at $B=9.8\,$T. It turns out that the in-plane dispersion of
Hofstadter bands is too small to be visible in the given energy scale.
Therefore, here we only show the case for $k_{x}=k_{y}=0$ without losing
generality. (f) Wannier plot associated with Hofstadter butterfly in (d). In
our convention, the filling is normalized to unity when the first conduction
band is filled.

Technically, the large moiré periodicity makes Landau levels split in a
fractal manner forming the well-known Hofstadter butterfly [86, 87, 88, 89,
90]. Our calculations [79] show that for $\theta\\!=\\!7\,^{\circ}$, the
fractal splitting of the first Landau level starts at $B\sim 100\,$T, which is
way beyond the usual experimentally achievable magnitude [see Fig. 2(d)]. So,
in practice, the renormalized Dirac Hamiltonian Eq. (5) suffices to get the
qualitative Landau level structure of the bulk ATG if
$\theta\\!\gtrapprox\\!6\,^{\circ}$. Under $B=9.8\,$T, the first and zeroth
Landau levels of the bulk ATG at $\theta\\!=\\!7\,^{\circ}$, extracted from
Hofstadter butterfly calculations, are separated in energy by a large gap
$\sim 66\,$meV, and the first and second Landau levels are also separated by
an global gap $\sim 2\,$meV (see Fig. 2(e)), such that several quantized Hall
conductivities can be observed. This is guaranteed by the Dirac dispersion
near charge neutrality point which applies to all $k_{z}$. Specifically, when
the $n$th Landau level is occupied and is separated from the higher one by a
global gap, the Hall conductivity should follow the sequence

$\sigma_{xy}=2(4n+2)\frac{e^{2}}{h}\cdot\frac{1}{2d_{0}},$ (6)

where the prefactor of 2 comes from spin degeneracy. This can be shown by
Wannier plot Fig. 2(f), where we only consider one spin sector. According to
the Diophantine equation, the slope of each line is precisely the Chern
number, or TKNN invariant of the filled Hofstadter bands [91, 92, 93, 94]. The
eight-fold degeneracy is due to two Dirac cones in each of two valleys for two
spin sectors.

Unlike the quantum Hall problem in 2D electron gas or in graphene, here the
degeneracy of each zeroth Landau level is: $D=8\,B\,S\,N_{z}/\phi_{0}$, where
$B$, $S$, $N_{z}$ denote the magnetic field, total cross-sectional area of the
bulk system, and number of primitive cells in the out-of-plane direction,
respectively. The factor of $8$ comes from valley, spin, and layer degeneracy.
The additional $k_{z}$ degeneracy may lead to a variety of unconventional
correlated and topological states for the partially occupied 3D zeroth Landau
level. For example, at partial integer fillings a gap may be opened up due to
strong $e$-$e$ interaction effects, such that the system may undergo a field-
induced metal-to-insulator transition through the spontaneous breaking of
flavor (valley, spin, and layer) symmetry. The additional $k_{z}$ degeneracy
may lead to further out-of-plane charge modulation driven by $e$-$e$ and/or
electron-phonon couplings, which may give rise to potential 3D topological
charge density wave state. Other new states of matter are also anticipated at
fractional fillings of the dispersionless 3D zeroth Landau level.

#### Conclusions and outlook –

In this work, we theoretically study the electronic structures, topological
properties, and Landau levels of a CVD-grown 3D moiré graphite superlattice,
the bulk ATG system. We have shown that bulk ATG can be seen as a series of
decoupled TBG with $k_{z}$ dependent moiré potential. Consequently, the flat
bands of TBG at the first magic angle are always present in bulk ATG at some
magic $k_{z}$ as long as $\theta<2\theta_{m}^{(1)}$. Moreover, multiple
topological flat bands of TBG at different magic angles emerge concomitantly
in bulk ATG at a set of magic vertical momenta. By introducing a small mass
term,a Weyl semimetal phase may anticipated in bulk ATG with tunable Weyl
nodal charges. Most saliently, when the twist angle $\theta\gtrapprox
3\rm{-}4\,^{\circ}$, bulk ATG is bestowed with a Dirac-like dispersion and a
quasi-2D Fermi surface at low energy, which gives rise to exactly
dispersionless 3D zeroth Landau level and weakly dispersive higher Landau
levels. When these 3D Landau levels are fully occupied, the system is expected
to realize 3D integer quantum Hall effects.

The physics of bulk ATG has even more to explore than TBG. First, when
$\theta<2\theta_{m}^{(1)}$, flat bands of magic-angle TBG occur at some magic
$k_{z}$, which coexist with linearly dispersive Dirac fermions as shown in
Fig. 2(c). The $e$-$e$ interaction effects and electron-phonon coupling
effects of such peculiar band structures are open questions. Second, when
$\theta<2\theta_{m}^{(2)}$, topological flat bands of TBG at multiple magic
angles are integrated into a single bulk ATG system, and unusual interacting
ground states would be anticipated in such multi-flat-band system. Last, when
$\theta$ is sufficiently larger than $2\theta_{m}^{(1)}$, e.g.,
$\theta\gtrapprox 3$-$4\,^{\circ}$, exactly dispersionless 3D zeroth Landau
levels would emerge. As a result, the intriguing quantum Hall physics
associated with 2D Landau levels would be bestowed with an additional $k_{z}$
degrees of freedom. Unusual quantum states of matter may emerge in the
partially occupied 3D Landau levels, which deserves further explorations.

###### Acknowledgements.

This work is supported by the National Key R & D program of China (grant No.
2020YFA0309601), the National Natural Science Foundation of China (grant No.
12174257), and the start-up grant of ShanghaiTech University.

## References

  * Bistritzer and MacDonald [2011a] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011a).
  * Lopes dos Santos  _et al._ [2007a] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007a).
  * Shallcross  _et al._ [2010] S. Shallcross, S. Sharma, E. Kandelaki, and O. A. Pankratov, Phys. Rev. B 81, 165105 (2010).
  * Uchida  _et al._ [2014] K. Uchida, S. Furuya, J.-I. Iwata, and A. Oshiyama, Phys. Rev. B 90, 155451 (2014).
  * 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).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Phys. Rev. Lett. 122, 106405 (2019).
  * Liu  _et al._ [2019a] J. Liu, J. Liu, and X. Dai, Phys. Rev. B 99, 155415 (2019a).
  * Song  _et al._ [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, Phys. Rev. Lett. 123, 036401 (2019).
  * Ahn  _et al._ [2019] J. Ahn, S. Park, and B.-J. Yang, Phys. Rev. X 9, 021013 (2019).
  * Po  _et al._ [2019] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Phys. Rev. B 99, 195455 (2019).
  * Bultinck  _et al._ [2020a] N. Bultinck, S. Chatterjee, and M. P. Zaletel, Phys. Rev. Lett. 124, 166601 (2020a).
  * Bultinck  _et al._ [2020b] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Phys. Rev. X 10, 031034 (2020b).
  * Wang  _et al._ [2021] J. Wang, Y. Zheng, A. J. Millis, and J. Cano, Physical Review Research 3, 023155 (2021).
  * Shi and Dai [2022] H. Shi and X. Dai, Phys. Rev. B 106, 245129 (2022).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018a).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Codecido  _et al._ [2019] E. Codecido, Q. Wang, R. Koester, S. Che, H. Tian, R. Lv, S. Tran, K. Watanabe, T. Taniguchi, F. Zhang, M. Bockrath, and C. N. Lau, Sci. Adv. 5, 10.1126/sciadv.aaw9770 (2019).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Nature 574, 653 (2019).
  * Stepanov  _et al._ [2020] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Nature 583, 375 (2020).
  * Saito  _et al._ [2020] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F. Young, Nature Physics 16, 926 (2020).
  * Liu  _et al._ [2021] X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek, and J. Li, Science 371, 1261 (2021).
  * Cao  _et al._ [2021a] Y. Cao, D. Rodan-Legrain, J. M. Park, N. F. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, and P. Jarillo-Herrero, science 372, 264 (2021a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, _et al._ , Nature 556, 80 (2018b).
  * Kerelsky  _et al._ [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, _et al._ , Nature 572, 95 (2019).
  * Jiang  _et al._ [2019] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule, J. Mao, and E. Y. Andrei, Nature 573, 91 (2019).
  * Xie  _et al._ [2019] Y. Xie, B. Lian, B. Jäck, X. Liu, C.-L. Chiu, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 572, 101 (2019).
  * Choi  _et al._ [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. von Oppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge, Nat. Phys. 15, 1174 (2019).
  * Serlin  _et al._ [2019] M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Science (2019).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Nat. Phys. 16, 725 (2020).
  * Andrei  _et al._ [2021] E. Y. Andrei, D. K. Efetov, P. Jarillo-Herrero, A. H. MacDonald, K. F. Mak, T. Senthil, E. Tutuc, A. Yazdani, and A. F. Young, Nat. Rev. Mater. 6, 201 (2021).
  * Liu and Dai [2021a] J. Liu and X. Dai, Nature Reviews Physics 3, 367 (2021a).
  * Stepanov  _et al._ [2021] P. Stepanov, M. Xie, T. Taniguchi, K. Watanabe, X. Lu, A. H. MacDonald, B. A. Bernevig, and D. K. Efetov, Phys. Rev. Lett. 127, 197701 (2021).
  * Nuckolls  _et al._ [2020] K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 588, 610 (2020).
  * Wu  _et al._ [2021] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Nature Materials (2021).
  * Das  _et al._ [2021] I. Das, X. Lu, J. Herzog-Arbeitman, Z.-D. Song, K. Watanabe, T. Taniguchi, B. A. Bernevig, and D. K. Efetov, Nature Physics 17, 710 (2021).
  * Pierce  _et al._ [2021] A. T. Pierce, Y. Xie, J. M. Park, E. Khalaf, S. H. Lee, Y. Cao, D. E. Parker, P. R. Forrester, S. Chen, K. Watanabe, _et al._ , Nature Physics 17, 1210 (2021).
  * Shen  _et al._ [2021] C. Shen, J. Ying, L. Liu, J. Liu, N. Li, S. Wang, J. Tang, Y. Zhao, Y. Chu, K. Watanabe, T. Taniguchi, R. Yang, D. Shi, F. Qu, L. Lu, W. Yang, and G. Zhang, Chinese Physics Letters 38, 047301 (2021).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Phys. Rev. Lett. 122, 246401 (2019).
  * Seo  _et al._ [2019] K. Seo, V. N. Kotov, and B. Uchoa, Phys. Rev. Lett. 122, 246402 (2019).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Phys. Rev. Lett. 124, 097601 (2020).
  * Wu and Das Sarma [2020] F. Wu and S. Das Sarma, Phys. Rev. Lett. 124, 046403 (2020).
  * Liu and Dai [2021b] J. Liu and X. Dai, Phys. Rev. B 103, 035427 (2021b).
  * Zhang  _et al._ [2020] Y. Zhang, K. Jiang, Z. Wang, and F. Zhang, Phys. Rev. B 102, 035136 (2020).
  * Hejazi  _et al._ [2021] K. Hejazi, X. Chen, and L. Balents, Phys. Rev. Research 3, 013242 (2021).
  * Kang and Vafek [2020] J. Kang and O. Vafek, Phys. Rev. B 102, 035161 (2020).
  * Chen  _et al._ [2021] B.-B. Chen, Y. D. Liao, Z. Chen, O. Vafek, J. Kang, W. Li, and Z. Y. Meng, Nat. Commun. 12, 5480 (2021).
  * Lu  _et al._ [2022] C. Lu, Y. Zhang, Y. Zhang, M. Zhang, C.-C. Liu, Y. Wang, Z.-C. Gu, W.-Q. Chen, and F. Yang, Phys. Rev. B 106, 024518 (2022).
  * Da Liao  _et al._ [2021] Y. Da Liao, J. Kang, C. N. Breiø, X. Y. Xu, H.-Q. Wu, B. M. Andersen, R. M. Fernandes, and Z. Y. Meng, Phys. Rev. X 11, 011014 (2021).
  * Bernevig  _et al._ [2021] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Phys. Rev. B 103, 205413 (2021).
  * Lian  _et al._ [2021] B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Phys. Rev. B 103, 205414 (2021).
  * Xie  _et al._ [2021] F. Xie, A. Cowsik, Z.-D. Song, B. Lian, B. A. Bernevig, and N. Regnault, Phys. Rev. B 103, 205416 (2021).
  * Soejima  _et al._ [2020] T. Soejima, D. E. Parker, N. Bultinck, J. Hauschild, and M. P. Zaletel, Phys. Rev. B 102, 205111 (2020).
  * Potasz  _et al._ [2021] P. Potasz, M. Xie, and A. H. MacDonald, Phys. Rev. Lett. 127, 147203 (2021).
  * Zhang  _et al._ [2021] X. Zhang, G. Pan, Y. Zhang, J. Kang, and Z. Y. Meng, Chinese Physics Letters 38, 077305 (2021).
  * Hofmann  _et al._ [2022] J. S. Hofmann, E. Khalaf, A. Vishwanath, E. Berg, and J. Y. Lee, Phys. Rev. X 12, 011061 (2022).
  * Parker  _et al._ [2021] D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel, and N. Bultinck, Phys. Rev. Lett. 127, 027601 (2021).
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig, Phys. Rev. Lett. 129, 047601 (2022).
  * Dean  _et al._ [2010] C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard, _et al._ , Nature nanotechnology 5, 722 (2010).
  * Liu  _et al._ [2022] C. Liu, Z. Li, R. Qiao, Q. Wang, Z. Zhang, F. Liu, Z. Zhou, N. Shang, H. Fang, M. Wang, _et al._ , Nature Materials 21, 1263 (2022).
  * Liu  _et al._ [2020] X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, Y. Ronen, H. Yoo, D. Haei Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Nature 583, 221 (2020).
  * Cao  _et al._ [2020] Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 10.1038/s41586-020-2260-6 (2020).
  * Shen  _et al._ [2020] C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, K. Watanabe, T. Taniguchi, R. Yang, Z. Y. Meng, D. Shi, O. V. Yazyev, and G. Zhang, Nature Physics 16, 520 (2020).
  * Liu  _et al._ [2019b] J. Liu, Z. Ma, J. Gao, and X. Dai, Phys. Rev. X 9, 031021 (2019b).
  * Koshino [2019] M. Koshino, Phys. Rev. B 99, 235406 (2019).
  * Lee  _et al._ [2019] J. Y. Lee, E. Khalaf, S. Liu, X. Liu, Z. Hao, P. Kim, and A. Vishwanath, Nature Communications 10, 5333 (2019).
  * Ledwith  _et al._ [2022] P. J. Ledwith, A. Vishwanath, and E. Khalaf, Phys. Rev. Lett. 128, 176404 (2022).
  * Wang and Liu [2022] J. Wang and Z. Liu, Phys. Rev. Lett. 128, 176403 (2022).
  * Zhang  _et al._ [2023] S. Zhang, B. Xie, Q. Wu, J. Liu, and O. V. Yazyev, Nano Letters 23, 2921 (2023).
  * Khalaf  _et al._ [2019] E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Phys. Rev. B 100, 085109 (2019).
  * Xie  _et al._ [2022] B. Xie, R. Peng, S. Zhang, and J. Liu, npj Computational Materials 8, 110 (2022).
  * Ledwith  _et al._ [2021] P. J. Ledwith, E. Khalaf, Z. Zhu, S. Carr, E. Kaxiras, and A. Vishwanath, Tb or not tb? contrasting properties of twisted bilayer graphene and the alternating twist $n$-layer structures ($n=3,4,5,\dots$) (2021).
  * Leconte  _et al._ [2022] N. Leconte, Y. Park, J. An, A. Samudrala, and J. Jung, 2D Materials 9, 044002 (2022).
  * Cea  _et al._ [2019] T. Cea, N. R. Walet, and F. Guinea, Nano letters 19, 8683 (2019).
  * Cao  _et al._ [2021b] Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 595, 526 (2021b).
  * Polshyn  _et al._ [2020] H. Polshyn, J. Zhu, M. Kumar, Y. Zhang, F. Yang, C. Tschirhart, M. Serlin, K. Watanabe, T. Taniguchi, A. MacDonald, _et al._ , Nature , 1 (2020).
  * Li  _et al._ [2022] Q. Li, B. Cheng, M. Chen, B. Xie, Y. Xie, P. Wang, F. Chen, Z. Liu, K. Watanabe, T. Taniguchi, S.-J. Liang, D. Wang, C. Wang, Q.-H. Wang, J. Liu, and F. Miao, Nature 609, 479 (2022).
  * Wang  _et al._ [2023] Z.-J. Wang, X. Kong, Y. Huang, J. Li, L. Bao, K. Cao, Y. Hu, J. Cai, L. Wang, H. Chen, Y. Wu, Y. Zhang, F. Pang, Z. Cheng, P. Babor, M. Kolibal, Z. Liu, Y. Chen, Q. Zhang, Y. Cui, K. Liu, H. Yang, X. Bao, H.-J. Gao, Z. Liu, W. Ji, F. Ding, and M.-G. Willinger, Nature Materials 10.1038/s41563-023-01632-y (2023).
  * [79] See Supplemental Information for: (a) details of the continuum model for bulk alternating twisted graphite, (b) detailed band structures and topological transitions in bulk alternating twisted graphite, (c) lattice relaxations and atomistic tight-binding model for bulk alternating twisted graphite, (d) detailed formalism of calculating Hofstadter butterfly spectra based on the continuum model, (e) more results about the Hofstadter butterfly spectra, (f) bound states associated with the spiral dislocation line.
  * Lopes dos Santos  _et al._ [2012] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 86, 155449 (2012).
  * Moon and Koshino [2013] P. Moon and M. Koshino, Physical Review B 87, 205404 (2013).
  * Hejazi  _et al._ [2019a] K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Balents, Phys. Rev. B 99, 035111 (2019a).
  * Bernevig  _et al._ [2007] B. A. Bernevig, T. L. Hughes, S. Raghu, and D. P. Arovas, Phys. Rev. Lett. 99, 146804 (2007).
  * Tang  _et al._ [2019] F. Tang, Y. Ren, P. Wang, R. Zhong, J. Schneeloch, S. A. Yang, K. Yang, P. A. Lee, G. Gu, Z. Qiao, _et al._ , Nature 569, 537 (2019).
  * Novoselov  _et al._ [2005] K. Novoselov, A. K. Geim, S. Morozov, D. Jiang, M. Katsnelson, I. Grigorieva, S. Dubonos, and A. Firsov, nature 438, 197 (2005).
  * Bistritzer and MacDonald [2011b] R. Bistritzer and A. H. MacDonald, Phys. Rev. B 84, 035440 (2011b).
  * de Gail  _et al._ [2011] R. de Gail, M. O. Goerbig, F. Guinea, G. Montambaux, and A. H. Castro Neto, Phys. Rev. B 84, 045436 (2011).
  * Hejazi  _et al._ [2019b] K. Hejazi, C. Liu, and L. Balents, Phys. Rev. B 100, 035115 (2019b).
  * Zhang  _et al._ [2019] Y.-H. Zhang, H. C. Po, and T. Senthil, Phys. Rev. B 100, 125104 (2019).
  * Lian  _et al._ [2020] B. Lian, F. Xie, and B. A. Bernevig, Phys. Rev. B 102, 041402 (2020).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Reviews of modern physics 82, 1959 (2010).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
  * Kohmoto [1989] M. Kohmoto, Phys. Rev. B 39, 11943 (1989).
  * Dana  _et al._ [1985] I. Dana, Y. Avron, and J. Zak, Journal of Physics C: Solid State Physics 18, L679 (1985).
  * Lopes dos Santos  _et al._ [2007b] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007b).
  * 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, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Comp. Phys. Comm. 271, 108171 (2022).
  * Los and Fasolino [2003] J. H. Los and A. Fasolino, Phys. Rev. B 68, 024107 (2003).
  * Nam and Koshino [2017] N. N. T. Nam and M. Koshino, Phys. Rev. B 96, 075311 (2017).
  * Angeli  _et al._ [2018] M. Angeli, D. Mandelli, A. Valli, A. Amaricci, M. Capone, E. Tosatti, and M. Fabrizio, Phys. Rev. B 98, 235137 (2018).
  * Cahill and Glauber [1969] K. E. Cahill and R. J. Glauber, Phys. Rev. 177, 1857 (1969).
  * Streda [1982] P. Streda, Journal of Physics C: Solid State Physics 15, L717 (1982).

Supplementary Information for “Magic momenta and three dimensional Landau
levels from a three dimensional graphite moiré superlattice”

###### Contents

  1. I S1. Continuum model for bulk alternating twisted graphite
    1. I.1 Band structures for $\theta=1.05^{\circ}$ at different $k_{z}$
    2. I.2 Gap closing and topological transition at the vicinity of magic momentum
  2. II S2. Lattice relaxation of graphite spiral
  3. III S3. Landau quantization of continuum model
    1. III.1 Hamiltonian in Landau level basis
    2. III.2 Displacement operators
    3. III.3 Commensurate condition
    4. III.4 Numerical implementation
    5. III.5 Wannier plot
  4. IV S4. Hofstadter butterfly spectra
  5. V S5. Bound states of the spiral dislocation line

##  I S1. Continuum model for bulk alternating twisted graphite

Away from the dislocation line, graphite spiral (GS) can be thought as bulk
alternating twisted multilayer graphene (ATG), as shown in Fig. 1 of the main
text. Bulk ATG consists of an infinitely large number of graphene layers
stacked with an alternating twist angle $\pm\theta$ between each pair of
adjacent graphene layers. Equivalently, bulk ATG can be viewed as a periodic
stacking of twisted bilayer graphene (TBG) along the out-of-plane direction
with lattice constant $2d_{0}$, where $d_{0}=3.35\,$Å is the distance between
two adjacent graphene layers. Within each stacking motif, TBG forms a moiré
pattern with spatial periodicity $L_{s}=a/(2\sin(\theta/2))$, where
$a=2.46\,$Å is the lattice constant of monolayer graphene.

To obtain the low-energy electronic structure of bulk ATG, we start from the
continuum model of TBG for valley $\mu$ at twist angle $\theta$ [95, 1]:

$\displaystyle\textrm{H}^{\mu}_{\theta}=\left[\begin{array}[]{cc}-\hbar
v_{F}(\mathbf{k}-\mathbf{K}_{1}^{\mu})\cdot\bm{\sigma}^{\mu}&T\\\
T^{\dagger}&-\hbar
v_{F}(\mathbf{k}-\mathbf{K}_{2}^{\mu})\cdot\bm{\sigma}^{\mu}\end{array}\right]$
(9)

where the Fermi velocity is $\hbar v_{F}/a\\!=\\!2.1354\,$eV [81]. The
diagonal blocks represent the low-energy Hamiltonian of the two layers of
graphene in valley $\mu=\pm 1$ near its Dirac points $\mathbf{K}_{1/2,\mu}$
with the Pauli matrices $\bm{\sigma}^{\mu}=(\mp\sigma_{x},\sigma_{y})$ and
$\mathbf{k}=(k_{x},k_{y})$. The intervalley coupling is omitted for the low-
energy states on the moiré length scale so that $\mu$ is a good quantum
number. The off-diagonal term is the moiré potential.
$T=U(\mathbf{r})\exp(-i\Delta\mathbf{K}_{\mu}\cdot\mathbf{r})$ with the shift
between two Dirac points in the same valley from two adjacent layer
$\Delta\mathbf{K}_{\mu}=\mathbf{K}_{2,\mu}-\mathbf{K}_{1,\mu}=\mu(0,k_{\theta})$
and

$\displaystyle U(\textbf{r})=u_{0}^{\prime}\left[\begin{array}[]{cc}\kappa
g(\mathbf{r})&g(\mathbf{r}-\mathbf{r}_{\rm{AB}})\\\
g(\mathbf{r}-\mathbf{r}_{\rm{AB}})&\kappa g(\mathbf{r})\end{array}\right]$
(12)

where $\kappa=u_{0}/u^{\prime}_{0}$ with intra-sublattice and inter-sublattice
interlayer tunneling amplitudes denoted by $u_{0}$ and $u_{0}^{\prime}$,
respectively. Here $k_{\theta}=4\pi/3L_{s}$. If out-of-plane corrugations are
taken into account, one finds that $u_{0}^{\prime}=0.0975\,$eV and
$u_{0}=0.0797\,$eV [5]. The phase factor is defined as
$g(\mathbf{r})\\!=\\!\sum^{3}_{j=1}e^{i\mathbf{q}_{j}\cdot\mathbf{r}}$,
$\mathbf{q}_{1}=k_{\theta}(0,1)$,
$\mathbf{q}_{2}=k_{\theta}(-\sqrt{3}/2,-1/2)$ and
$\mathbf{q}_{3}=k_{\theta}(\sqrt{3}/2,-1/2)$ with
$\mathbf{r}_{\rm{AB}}=(\sqrt{3}L_{s}/3,0)$. Up to a multiplicative factor
$\hbar v_{F}k_{\theta}$, the continuum model of TBG is fully characterized by
the dimensionless constants $\alpha_{0}=u^{\prime}_{0}/(\hbar
v_{F}k_{\theta})$ and $\kappa=u_{0}^{\prime}/u_{0}$.

Bulk ATG is formed by stacking vertically TBG in the $z$-direction. The
tunnelling between two adjacent TBGs is given by a $k_{z}$-dependent moiré
potential $T^{\prime}$ with possibly different dimensionless ratio
$\alpha_{0}^{\prime}$ and $\kappa^{\prime}$. In the second quantization, the
Hamiltonian is written:

$\displaystyle\hat{H}^{\mu}=\sum_{z}\sum_{l=1,2}\left(H^{\mu}_{l}\right)_{\alpha\alpha^{\prime}}\hat{c}^{\dagger}_{l,\alpha}(z)\hat{c}_{l,\alpha^{\prime}}(z)+\sum_{z}\left(T\right)_{\alpha\alpha^{\prime}}\hat{c}^{\dagger}_{1,\alpha}(z)\hat{c}_{2,\alpha^{\prime}}(z)+\left(T^{\prime}\right)_{\alpha\alpha^{\prime}}\hat{c}^{\dagger}_{1,\alpha}(z)\hat{c}_{2,\alpha^{\prime}}(z+d_{0}+d_{0}^{\prime})+\rm{H.c.}$
(13)

where $l=1,2$ and $\alpha=A,B$ represent the layer and the sublattice indices
within TBG unit, and $H^{\mu}_{l}$ is the diagonal term of Eq. (9).
$T_{\alpha,\alpha^{\prime}}$ in Eq. (13) is just the off-diagonal moiré
potential term within a TBG unit as expressed in Eqs. (9)-(12). Without
$T^{\prime}$, the Hamiltonian Eq. (13) is just a sum of decoupled TBG at
different $z$. This form is reminiscent of the tight-binding model of quasi-1D
chain, where the hopping between nearest primitive cells is precisely given by
$T^{\prime}$ and the lattice constant in the out-of-plane direction is
$d_{0}+d_{0}^{\prime}$. Here $d_{0}$ denotes the interlayer distance within a
TBG unit, and $d_{0}^{\prime}$ denotes the interlayer distance between two
neighboring TBG units. By virtue of lattice translational symmetry in the $z$
direction, we use the Fourier transform

$\displaystyle\hat{c}_{l,\alpha}(z)=\frac{1}{\sqrt{N_{z}}}\sum_{k_{z}}e^{ik_{z}z}\hat{c}_{l,\alpha}(k_{z})$
(14)

to rewrite the Hamiltonian such that

$\displaystyle\hat{H}^{\mu}=\sum_{k_{z}}\sum_{l=1,2}\left(H^{\mu}_{l}\right)_{\alpha\alpha^{\prime}}\hat{c}^{\dagger}_{l,\alpha}(k_{z})\hat{c}^{\dagger}_{l,\alpha^{\prime}}(k_{z})+\sum_{k_{z}}\left(T\right)_{\alpha\alpha^{\prime}}\hat{c}^{\dagger}_{1,\alpha}(k_{z})\hat{c}_{2,\alpha^{\prime}}(k_{z})+\left(T^{\prime}\right)_{\alpha\alpha^{\prime}}\hat{c}^{\dagger}_{1,\alpha}(k_{z})\hat{c}_{2,\alpha^{\prime}}(k_{z})e^{ik_{z}(d_{0}+d_{0}^{\prime})}+\rm{H.c.}\
.$ (15)

In principal, dimerization could occur such that $d_{0}\neq d_{0}^{\prime}$,
such that the system at each in-plane moiré wavevector resembles a Su-
Schrieffer-Hegger chain. However, our transmission electron microscopy
measurements indicate that dimerization is absent in the bulk ATG system.
Therefore, hereafter we take $d_{0}=d_{0}^{\prime}=3.35\,$Å, and
$T=T^{\prime}$, thus we obtain the Hamiltonian Eq. (1) in the main text.

Figure 4: (a) The Fermi velocity of TBG near moiré $K$ point as a function of
twist angle. The red spots mark the magic angles. The $k_{z}$-dependent in-
plane Fermi velocities of bulk ATG with twist angle (b) $0.60\,^{\circ}$, (c)
$0.906\,^{\circ}$, and (d) $2.00\,^{\circ}$. The red spots mark the magic
momenta at which the in-plane Fermi velocities vanish.

When $u_{0}=u_{0}^{\prime}$, the magic angle of TBG is solely determined by a
dimensionless coupling constant $\alpha_{0}(\theta)=u_{0}^{\prime}/(\hbar
v_{F}k_{\theta})$. In bulk ATG, at each $k_{z}$ one can define a renormalized
coupling constant
$\alpha(\theta,k_{z})=u_{0}^{\prime}|1+e^{i2k_{z}d_{0}}|/(\hbar
v_{F}k_{\theta})$, such that at twist angle $\theta$, when $k_{z}=k_{z}^{*}$
satisfies

$\alpha(\theta,k_{z}^{*})=\alpha_{0}(\theta_{m}^{(1)})\;,$ (16)

one would obtain the flat bands at “magic momentum” $k_{z}^{*}$ which are
exactly the same as those of TBG at the first magic angle $\theta_{m}^{(1)}$.
Equivalently, the magic momentum condition is simply written as

$2\cos(k_{z}^{*}d_{0})=\frac{\sin{\theta/2}}{\sin{\theta_{m}^{(1)}/2}}$ (17)

For small twist angle $\theta$, $\sin{\theta/2}\approx\theta/2$. In order to
have real valued solution of $k_{z}^{*}$ of Eq. (17), the twist angle of bulk
ATG needs to satisfy $\theta<2\theta_{m}^{(1)}$. More generally, when $\theta$
is smaller than twice of the $n$th magic angle of TBG, i.e.,
$\theta<2\theta_{m}^{(n)}$, one can always find the $s$th magic momentum
$k_{z}^{(s),*}$ which satisfies

$2\cos(k_{z}^{(s),*}d_{0})=\frac{\sin{\theta/2}}{\sin{\theta_{m}^{(s)}/2}}$
(18)

where $\theta^{(s)}_{m}$ is the $s$th magic angle of TBG with $s=1,2,...,n$.
Technically, we define magic angles of TBG as the twist angles at which the
Fermi velocity of Dirac cone vanishes. In Fig. 4(a), we present the Fermi
velocity of TBG near moiré $\mathbf{K}$ point as a function of $\theta$. The
red spot represents the magic angles from the first to the fifth.

In Fig. 2(a) in the main text, we plot the values of the magic momenta
$\\{k_{z}^{(s),*}\\}$ (for $s=1,...,5$) as a function of twist angle $\theta$.
We also calculate the $k_{z}$-dependent in-plane Fermi velocity of bulk ATG
with three specific twist angles $\theta=0.60\,^{\circ}$, $0.906\,^{\circ}$
and $2.00\,^{\circ}$ in Fig. 4(b), (c) and (d), respectively. When
$\theta=0.60\,^{\circ}$, from Fig. 2(a) in the main text, we can learn that
there will be three $k_{z}$ values which satisfy the magic-momenta condition
(Eq. (18)), which are marked by red dots in Fig. 4(b). Same in Fig. 4(c),
there are two spots representing the two magic momenta at which the in-plane
Fermi velocities vanish. In Fig. 4(d), there is only one spot representing one
magic momentum.

###  I.1 Band structures for $\theta=1.05^{\circ}$ at different $k_{z}$

Figure 5: Band structures in bulk ATG with twist angle $\theta=1.05^{\circ}$
in the $(k_{x},k_{y})$-plane for (a) $k_{z}d_{0}\\!=\\!0$, (b) $\pi/6$, (c)
$\pi/4$, (d) $\pi/3$, (e) $5\pi/12$ and (f) $\pi/2$.

In Fig. 5(a)-(f), we explicitly show the evolution of 2D band structures in
the $(k_{x},k_{y})$-plane for $k_{z}d_{0}$ varying from 0 to $\pi/2$ in bulk
ATG with twist angle $\theta=\theta_{m}^{(1)}=1.05\,^{\circ}$. The previous
argument shows that the flat bands of TBG at the first magic angle emerge at
$k_{z}^{(1),*}d_{0}=\pi/3$, which is confirmed by numerical calculations as
shown in Fig. 5(d). At $k_{z}d_{0}=\pi/2$, the phase factor
$(1+e^{2ik_{z}d_{0}})$ vanishes (see Eq. (1) in the main text), such that the
band structure is just those of two decoupled single layer graphene folded
into the moiré Brillouin zone, as shown in Fig. 5(f).

Figure 6: (a) Band gaps $\Delta_{1,2}$ at $\Gamma$ point for bulk ATG with
twist angle $2.2\,^{\circ}$ at different $k_{z}$. The blue line represents the
band gap at $\Gamma$ point between the flat bands and the remote bands. The
orange line represent the band gap at $\Gamma$ point between two remote bands.
In (b), we show a zoomed-in view near the critical point. The dashed line and
the black dots mark the topological transition point of the flat band. The
Chern numbers are labeled as $C$. We also mark with red dots the kink spot in
the evolution of $\Delta_{1}$, which is due to band order switching. This is
clearly shown in (c) the evolution of the energy of three relevant bands at
$\Gamma$ point, where the band inversion points are annotated by red and black
dot, one-to-one corresponding with those in (b). The flat bands and remotes
bands near $\Gamma$ point for bulk ATG with (d) $k_{z}d_{0}/\pi=0.093$, (e)
$0.097$ and (f) $0.1$. The solid lines represent energy band from the $\mu=+$
valley and the dashed lines represent energy band from the other valley.

###  I.2 Gap closing and topological transition at the vicinity of magic
momentum

Figure 7: (a) Band gaps $\Delta_{1,2}$ at $\Gamma$ point for bulk ATG with
twist angle $2\,^{\circ}$ at different $k_{z}$. The blue line represents the
band gap at $\Gamma$ point between the flat bands and the remote bands. The
orange line represent the band gap at $\Gamma$ point between two remote bands.
In (b), we show a zoomed-in view near the two critical points. The dashed
lines and the black dots mark the topological transition points of the flat
band. The Chern numbers are labeled as $C$. We also mark with red dots one of
the kink spots in the evolution of $\Delta_{1}$, which is due to band order
switching. This is shown in (c) the evolution of the energy of three relevant
bands at $\Gamma$ point, where three band inversion points are annotated by
red and black dot, one-to-one corresponding with those in (b). The flat bands
and remotes bands near $\Gamma$ point for bulk ATG with (c)
$k_{z}d_{0}/\pi=0.10$, (d) $0.11$ and (e) $0.14$. The solid lines represent
energy band from the $\mu=+$ valley and the dashed lines represent energy band
from the other valley.

In this section, we provide the detailed information of the band closings at
the phase boundaries of the topological transitions. As we mentioned in the
main text, we establish an intriguing topological phase diagram by varying
both twist angle $\theta$ and $k_{z}$. When the twist angle is fixed between
$\theta_{c,0}\\!=\\!2.3^{\circ}$ and $\theta_{c,1}\\!=\\!2.15^{\circ}$, the
Brillouin zone can be divided into two topologically different parts. When the
twist angle is smaller than $\theta_{c,1}$, there will be more topological
transitions while varying $k_{z}$ inside the 3D Brillouin zone. Here, we
provide the detailed electronic band structures at the near the topological
transitions.

For bulk ATG with a twist angle of $\theta=2.2^{\circ}$, the Brillouin zone
can be partitioned into two distinct regions: the Chern number of the highest
valence band of $K$ valley equals zero if
$k_{z}\in[-k^{(1),*}_{z},k^{(1),*}_{z}]$ and equals one otherwise. We
introduce a small mass term with a value of $m_{z}\\!=\\!0.001\,$eV and obtain
the band structures of bulk ATG at the vicinity of $k_{z}^{(1),*}$. In Fig.
6(a), we present the direct band gaps at $\Gamma$ point for bulk ATG at a
function of $k_{z}$. The blue line represents the gap between the highest
valence band and the band immediately below it (dubbed as the first remote
band) at $\Gamma$ point, denoted as $\Delta_{1}$ in Fig. 6(a); while the
orange line represents the gap between the first remote band and the second
remote band at $\Gamma$ point, denoted as $\Delta_{2}$ in Fig. 6(a). Fig. 6(b)
is a zoomed-in view of the two direct gaps at $\Gamma$ near the gap closing
points. The dashed line marks the critical $k_{z}$ value marking the
topological transition where $\Delta_{1}$ vanishes. The corresponding Chern
number $C$ of the highest valence band of valley $K$ aside the transition
point is also shown. We see that $\Delta_{1}$ vanishes at
$k_{z}\\!\approx\\!0.097\pi/d_{0}$, across which the Chern number changes by
unity. Furthermore, near the phase boundary, another gap closure occurs
between two remote bands, but the Chern number of the highest valence band
remains unaffected as expected. In Fig. 6(d-f), we explicitly present the band
structures in the $(k_{x},k_{y})$-plane of the bulk ATG with
$k_{z}d_{0}/\pi\\!=\\!0.093$, $0.097$ and $0.1$, respectively, which
explicitly demonstrates the gap closure at $k_{z}d_{0}/\pi\\!=\\!0.097$.

Furthermore, in the case of bulk ATG with a twist angle of
$\theta=2.0^{\circ}$, there are four topologically different regions. In Fig.
7(a), we show the behavior of the direct band gap at $\Gamma$ point. Again,
the blue line represents the gap between the highest valence band and the
first remote band (denoted as $\Delta_{1}$), while the orange line represents
the band gap between the two remote bands (denoted by $\Delta_{2}$). Fig. 7(b)
is a magnified view near two critical points, with the values of
$k^{(2),*}_{z}d_{0}/\pi=0.115$ and of $k^{(1),*}_{z}d_{0}/\pi=0.166$ at which
$\Delta_{1}$ vanishes. While increasing $k_{z}$ from zero, the Chern number
undergoes a transition from one to zero, and subsequently form zero back to
one. In Fig. 7(d-f), we present the band structures in the
$(k_{x},k_{y})$-plane with $k_{z}d_{0}/\pi\\!=\\!0.10$, $0.11$ and $0.14$,
respectively.

##  II S2. Lattice relaxation of graphite spiral

In this section we study the lattice distortions of bulk ATG using direct
molecular dynamics simulations. Specifically, the structural relaxation is
calculated utilizing Large-scale Atomic-Molecular Massively Parallel
Simulation (LAMMPS) [96], employing Long-range Carbon Bond-Order Potential
(LCBOP) [97]. The initial structure is a moiré supercell with periodic
boundary condition applied in both the in-plane and out-of-plane directions.
The twist angles between two adjacent layers, according to the commensurate
condition, are set to $\pm 1.05\,^{\circ}$, $\pm 1.31\,^{\circ}$, $\pm
1.58\,^{\circ}$ and $\pm 1.83\,^{\circ}$. The lattice relaxation calculations
are carried out using the steepest descent algorithm. In Fig. 8(a), we present
the relative interlayer shift in the moiré supercell. The lateral lattice
distortions in bulk ATG are qualitatively the same as those in TBG [98, 99].
Due to the presence of the translational symmetry along the $z$-direction and
mirror symmetry, it turns out that there is no out-of-plane distortion in bulk
ATG. The system under different twist angle exhibits the lowest total energy
always at an interlayer distance of $d_{0}=3.35$ Å so that the vertical
lattice constant is $2d_{0}=6.7\,$Å. We have not observed any tendency of
dimerization between neighboring layers, consistent with the TEM results.

Based on the relaxed lattice structure, we further calculate the electronic
band structures of bulk ATG at the first magic angle $\theta=1.05\,^{\circ}$
using an atomistic Slater-Koster tight-binding model. To be specific, the
hopping amplitude between two $p_{z}$ orbitals at different sites is expressed
as [81]

$t(\textbf{d})=V_{\sigma}\left(\frac{\textbf{d}\cdot\hat{\textbf{z}}}{d}\right)+V_{\pi}\left[1-\left(\frac{\textbf{d}\cdot\hat{\textbf{z}}}{d}\right)^{2}\right],$
(19)

where $V_{\sigma}\\!=\\!V^{0}_{\sigma}e^{-(r-d_{0})/\delta_{0}}$ and
$V_{\pi}\\!=\\!V^{0}_{\pi}e^{-(r-a_{0})/\delta_{0}}$.
$\textbf{d}=(d_{x},d_{y},d_{z})$ is the displacement vector between two sites.
$d_{c}=3.35\,$Å is the interlayer distance. $a_{0}=a/\sqrt{3}=1.42\,$Å is the
distance between two nearest neighbor carbon atoms. $\delta_{0}=0.184a$, with
$a=2.46\,$Å. $V^{0}_{\sigma}=0.48\,$eV and $V^{0}_{\pi}=-2.7\,$eV. The band
structures of bulk ATG at $\theta=1.05\,^{\circ}$ at several fixed values of
$k_{z}$ are showed in Fig. 8(b)-(d). The band width of flat bands varies with
$k_{z}$, and reaches a minimum at $k_{z}d_{0}\approx 0.35\pi$. This is
slightly different from the value obtained from the magic-momentum condition
$k_{z}d_{0}\approx\pi/3$ (Eq. (1) in the main text) because the criteria for
the magic angle is slightly different, i. e., the vanishing Fermi velocity
around Dirac points does not exactly coincide with the condition of minimal
overall bandwidth.

Figure 8:  (a) The relative interlayer in-plane distortion field of bulk ATG
with twist angle $\theta=1.05^{\circ}$. The colorbar represents the amplitude
of the distortion fields. The arrows represent the direction of the distortion
fields. We present the band structure of bulk ATG at $\theta=1.05^{\circ}$
with fully relaxed lattice structure at (b) $k_{z}d_{0}/\pi=0$, (c) $0.35$ (c)
and (d) $0.5$.

##  III S3. Landau quantization of continuum model

We consider the Landau quantization of the continuum model for bulk ATG in the
presence of an out-of-plane magnetic field $B\hat{\textbf{z}}$. The associated
vector potential can always be chosen to be in-plane such that $k_{z}$ remains
a good quantum number. Therefore, the prefactor $1+\exp(i\phi)$ with
$\phi=2k_{z}d_{0}$ is totally inert in the following derivations. Therefore,
the derivations of the Landau quantization in bulk ATG are exactly the same as
those of TBG but with $k_{z}$-dependent moiré potential. Technically, we
follow the procedures given in [86], which are also adopted in [88, 89]. The
only difference is the gauge choice in the initial Hamiltonian.

###  III.1 Hamiltonian in Landau level basis

In our set-up, we start with the continuum Hamiltonian

$\displaystyle\textrm{H}^{\mu}_{\theta}=\left[\begin{array}[]{cc}-\hbar
v_{F}(\mathbf{k}-\mathbf{K}_{1,\mu})\cdot\bm{\sigma}^{\mu}&T(1+e^{i\phi})\\\
T^{\dagger}(1+e^{-i\phi})&-\hbar
v_{F}(\mathbf{k}-\mathbf{K}_{2,\mu})\cdot\bm{\sigma}^{\mu}\end{array}\right]$
(22)

where $\mathbf{k}=(k_{x},k_{y})$ and $k_{z}$ appearing in $\phi$ is irrelevant
in the following derivations. We restore all the units and expand the moiré
potential

$\begin{split}\bm{\sigma}^{\mu=\pm}&=(\mp\sigma_{x},\sigma_{y})\\\
\mathbf{K}_{l,\mu}&=k_{\theta}\left(\frac{\sqrt{3}}{2},\frac{(-1)^{l}\mu}{2}\right)\\\
T_{\mu}&=\sum^{3}_{j=1}T_{j,\mu}e^{i\mu\mathbf{Q}_{j}\cdot\mathbf{r}}\\\
T_{j,\mu}&=\left[\begin{array}[]{cc}u_{0}&u_{0}^{\prime}w^{\mu(j-1)}\\\
u_{0}^{\prime}w^{-\mu(j-1)}&u_{0}\end{array}\right]\\\
\mathbf{Q}_{1}&=(0,0),\quad\mathbf{Q}_{2}=\sqrt{3}k_{\theta}(-1/2,-\sqrt{3}/2),\quad\mathbf{Q}_{3}=\sqrt{3}k_{\theta}(1/2,-\sqrt{3}/2)\end{split}$
(23)

with $w=\exp(i2\pi/3)$. In the following, we only show the details of
calculations for $\mu=+$ while that of $\mu=-$ is exactly the same. We will
give the final results for both valleys at the end of this section.

Now we introduce an out-of-plane magnetic field via vector potential in the
Landau gauge $\mathbf{A}=(-By,0,0)$. The diagonal part Eq. (22) can be easily
written in the Landau level basis while the off-diagonal term cannot. Let us
first focus on the diagonal part. For $\mu=+$, we Peierls substitute $k_{x}$
by $k_{x}-eBy/\hbar$ then

$\displaystyle-\hbar
v_{F}(\mathbf{k}-\mathbf{K}_{l,+})\cdot\bm{\sigma}^{+}=\frac{\hbar\sqrt{2}v_{F}}{\ell_{\rm{B}}}\left[\begin{array}[]{cc}0&\hat{a}^{\dagger}(k_{x})\\\
\hat{a}(k_{x})&0\end{array}\right]+\hbar
v_{F}\left[\begin{array}[]{cc}0&-K_{l,+,x}-iK_{l,+,y}\\\
-K_{l,+,x}+iK_{l,+,y}&0\end{array}\right]$ (28)

where the magnetic length $\ell_{\rm{B}}=\sqrt{\hbar/eB}$ and the ladder
operators

$\begin{split}\hat{a}(k_{x})&=\frac{\ell_{\rm{B}}}{\sqrt{2}}\left(k_{x}-\frac{y}{\ell_{\rm{B}}^{2}}-ik_{y}\right)\\\
\hat{a}^{\dagger}(k_{x})&=\frac{\ell_{\rm{B}}}{\sqrt{2}}\left(k_{x}-\frac{y}{\ell_{\rm{B}}^{2}}+ik_{y}\right)\\\
\rm{and}&\quad[\hat{a}(k_{x}),\hat{a}^{\dagger}(k_{x})]=1\end{split}.$ (29)

We introduce the Landau level basis

$\begin{split}|l,\alpha;n,k_{x}\rangle&=\hat{f}^{\dagger}_{l,\alpha;n,k_{x}}|0\rangle\\\
\rm{with}&\quad|l,\alpha;n,k_{x}\rangle=\frac{(\hat{a}^{\dagger}(k_{x}))^{n}}{\sqrt{n!}}|l,\alpha;0,k_{x}\rangle\\\
&\quad\hat{a}(k_{x})|l,\alpha;0,k_{x}\rangle=0\end{split}$ (30)

where $\alpha=A,B$ is the sublattice and $l=1,2$ is the layer index. Then, the
diagonal part can be suitably written as

$\begin{split}\hat{H}^{\mu=+}_{D}=&\sum_{l,n,k_{x}}\hbar
v_{F}(-K_{l,+,x}-iK_{l,+,y})\
\hat{f}^{\dagger}_{l,A;n,k_{x}}\hat{f}_{l,B;n,k_{x}}+\rm{H.c.}\\\
&+\sum_{l,n,k_{x}}\frac{\sqrt{2(n+1)}\hbar v_{F}}{\ell_{\rm{B}}}\
\hat{f}^{\dagger}_{l,A;n+1,k_{x}}\hat{f}_{l,B;n,k_{x}}+\rm{H.c.}\end{split}.$
(31)

Now we write the off-diagonal part $\hat{H}^{\mu=+}_{M}$ in the Landau level
basis. To prepare the basis transformation, it is worthwhile to rewrite
$\hat{H}^{\mu=+}_{M}$ in the second-quantized form explicitly showing the
plane-wave basis, which is defined as

$\begin{split}|l,\alpha;\mathbf{k}\rangle&=\hat{c}^{\dagger}_{l,\alpha}(\mathbf{k})|0\rangle\\\
\langle\mathbf{r}|l,\alpha;\mathbf{k}\rangle&=\frac{1}{\sqrt{S}}e^{i\mathbf{k}\cdot\mathbf{r}}\chi_{l,\alpha}\end{split}$
(32)

with $\chi_{l,\alpha}$ represents the spinor part of the wavefunction and $S$
in the total area. In this basis, $\hat{H}^{\mu=+}_{M}$ is written as:

$\hat{H}^{\mu=+}_{M}=\sum_{\mathbf{k},\alpha\beta}\sum_{j=1}^{3}(1+e^{i\phi})(T_{j,+})_{\alpha\beta}\
\hat{c}^{\dagger}_{1,\alpha}(\mathbf{k}+\mathbf{Q}_{j})\hat{c}_{2,\alpha}(\mathbf{k})+\rm{H.c.}.$
(33)

By the completeness of the Landau level basis,

$|l,\alpha;\mathbf{k}\rangle=\sum_{n}|l,\alpha;n,k_{x}\rangle\langle
l,\alpha;n,k_{x}|l,\alpha;\mathbf{k}\rangle.$ (34)

The formula for the basis transformation is thus

$\hat{c}_{l,\alpha}(\mathbf{k})=\sum_{n}\hat{f}_{l,\alpha;n,k_{x}}\tilde{H}_{n}(k_{y})e^{-ik_{x}k_{y}\ell_{\rm{B}}^{2}}$
(35)

with $\tilde{H}_{n}(k_{y})$ is the Fourier transform of the $n$-th Hermite
function $H_{n}(y)$, i.e., the wavefunction of 1D quantum harmonic oscillator

$\tilde{H}_{n}(k_{y})=\int\frac{dy}{\sqrt{L_{y}}}e^{-ik_{y}y}H_{n}(y).$ (36)

So, the moiré potential in the Landau level basis reads

$\hat{H}^{\mu=+}_{M}=\sum_{j=1}^{3}\sum_{m,n,k_{x}}(1+e^{i\phi})F_{mn,k_{x}}(\mathbf{Q}_{j})\sum_{\alpha\beta}\left(T_{j,+}\right)_{\alpha\beta}\
\hat{f}^{\dagger}_{1,\alpha;m,k_{x}+Q_{j,x}}\hat{f}_{2,\alpha;n,k_{x}}+\rm{H.c.}$
(37)

where we use a shorthand notation

$F_{mn,k_{x}}(\mathbf{Q}_{j})=e^{i(k_{x}+Q_{j,x})Q_{j,y}\ell_{\rm{B}}^{2}}\
\sum_{k_{y}}e^{ik_{y}Q_{j,x}\ell_{\rm{B}}^{2}}\
\tilde{H}_{m}^{*}(k_{y}+Q_{j,y})\tilde{H}_{n}(k_{y}).$ (38)

###  III.2 Displacement operators

We now calculate the phase factor $F_{mn,k_{x}}(\mathbf{Q}_{j})$. Transforming
back $\tilde{H}_{n}(k_{y})$ to $H_{n}(y)$, the factor
$F_{mn,k_{x}}(\mathbf{Q}_{j})$ becomes

$F_{mn,k_{x}}(\mathbf{Q}_{j})=e^{i(k_{x}+Q_{j,x})Q_{j,y}\ell_{\rm{B}}^{2}}\
\langle m,0|e^{iQ_{j,y}}|n,-Q_{j,x}\rangle$ (39)

where we use the notation in Eq. (30) omitting irrelevant the $(l,\alpha)$
spinor part. Here, $k_{x}$ in $|n,k_{x}\rangle$ plays the role of the
cyclotron center, which sits at $y=k_{x}\ell_{\rm{B}}^{2}$. One can move the
cyclotron center back to the origin $y=0$ by applying a so-called displacement
operator on the state. It turns out the phase factor itself can be written in
terms of displacement operator such that one can calculate it analytically.

We provide here a brief self-contained introduction on how displacement
operator enters in our calculations. More details can be found in the classic
paper [100]. For the basis $|n,k_{x}\rangle$, we define

$D(\alpha)=e^{\alpha\hat{a}^{\dagger}(0)-\alpha^{*}\hat{a}(0)}$ (40)

with $\alpha\in\mathbb{C}$. One can check that

$\hat{a}(0)D(\alpha)|0,0\rangle=\alpha D(\alpha)|0,0\rangle$ (41)

which means that $D(\alpha)|0,0\rangle$ is an eigenstate of $\hat{a}(0)$ with
eigenvalue $\alpha$. The two above definitions for displacement operator are
equivalent. Other properties of $D$-operator are

$D(\alpha)D(\beta)=D(\alpha+\beta)e^{\frac{\alpha\beta^{*}-\alpha^{*}\beta}{2}}$
(42)

and the commutation relation

$D(\alpha)(\hat{a}^{\dagger}(0)+\alpha^{*})=\hat{a}^{\dagger}(0)D(\alpha),$
(43)

both of which can be proved by direct calculations [100].

In our problem, we notice that

$\hat{a}(-Q_{j,x})|0,-Q_{j,x}\rangle=0\quad\Leftrightarrow\quad\hat{a}(0)|0,-Q_{j,x}\rangle=\frac{Q_{j,x}\ell_{\rm{B}}}{\sqrt{2}}|0,-Q_{j,x}\rangle\quad\Leftrightarrow\quad|0,-Q_{j,x}\rangle=D\left(\frac{Q_{j,x}\ell_{\rm{B}}}{\sqrt{2}}\right)|0,0\rangle.$
(44)

So, as implied by the name, $|0,-Q_{j,x}\rangle$ is obtained by moving the
cyclotron center of $|0,0\rangle$ from the origin to
$y=-Q_{j,x}\ell_{\rm{B}}^{2}$. Also,

$\begin{split}D\left(\frac{Q_{j,x}\ell_{\rm{B}}}{\sqrt{2}}\right)&=e^{ik_{y}Q_{j,x}\ell_{\rm{B}}^{2}}\\\
D\left(-i\frac{Q_{j,y}\ell_{\rm{B}}}{\sqrt{2}}\right)&=e^{iQ_{j,y}y}.\end{split}$
(45)

Meanwhile,

$D\left(\frac{Q_{j,x}\ell_{\rm{B}}}{\sqrt{2}}\right)\hat{a}^{\dagger}(0)=\hat{a}^{\dagger}\left(\frac{Q_{j,x}\ell_{\rm{B}}}{\sqrt{2}}\right)D\left(\frac{Q_{j,x}\ell_{\rm{B}}}{\sqrt{2}}\right).$
(46)

Therefore, the matrix element in Eq. (39) can be written in terms of
displacement operators

$\displaystyle\langle m,0|e^{iQ_{j,y}}|n,-Q_{j,x}\rangle$
$\displaystyle=\langle
m,0|D\left(-i\frac{Q_{j,y}\ell_{\rm{B}}}{\sqrt{2}}\right)\frac{[\hat{a}^{\dagger}(-Q_{j,x})]^{n}}{\sqrt{n!}}D\left(\frac{Q_{j,x}\ell_{\rm{B}}}{\sqrt{2}}\right)|0,0\rangle$
$\displaystyle=\langle
m,0|D\left(-i\frac{Q_{j,y}\ell_{\rm{B}}}{\sqrt{2}}\right)D\left(\frac{Q_{j,x}\ell_{\rm{B}}}{\sqrt{2}}\right)|n,0\rangle$
$\displaystyle=\langle
m,0|D\left(\frac{(Q_{j,x}-iQ_{j,y})\ell_{\rm{B}}}{\sqrt{2}}\right)|n,0\rangle
e^{-i\frac{Q_{j,x}Q_{j,y}\ell_{\rm{B}}^{2}}{2}}.$

So, the phase factor becomes

$F_{mn,k_{x}}(\mathbf{Q}_{j})=e^{i\frac{(2k_{x}+Q_{j,x})Q_{j,y}\ell_{\rm{B}}^{2}}{2}}\
\langle m,0|D(z_{j})|n,0\rangle\equiv
e^{i\frac{(2k_{x}+Q_{j,x})Q_{j,y}\ell_{\rm{B}}^{2}}{2}}F^{0}_{mn}(\mathbf{Q}_{j})$
(47)

with $z_{j}=\ell_{\rm{B}}/\sqrt{2}(Q_{j,x}-Q_{j,y})$. The matrix element
$F^{0}_{mn}(\mathbf{Q}_{j})$ has an analytical form involving generalized
Laguerre polynomial $\\{L^{(\alpha)}_{n}\\}$ with $\alpha>-1$:

$F^{0}_{mn}(\mathbf{Q}_{j})=\left\\{\begin{array}[]{ll}e^{-\frac{|z_{j}|^{2}}{2}}z_{j}^{m-n}\sqrt{\frac{n!}{m!}}L^{(m-n)}_{n}(|z_{j}|^{2})&\mbox{if
$m\geq n$, }\\\
e^{-\frac{|z_{j}|^{2}}{2}}(-z_{j}^{*})^{n-m}\sqrt{\frac{m!}{n!}}L^{(n-m)}_{m}(|z_{j}|^{2})&\mbox{if
$m<n$.}\end{array}\right.$ (48)

The moiré potential Eq. (37) becomes

$\hat{H}^{\mu=+}_{M}=\sum_{j=1}^{3}\sum_{mn,\alpha\beta}(1+e^{i\phi})\left(T_{j,+}\right)_{\alpha\beta}F^{0}_{mn}(\mathbf{Q}_{j})e^{i\frac{Q_{j,x}Q_{j,y}\ell_{\rm{B}}^{2}}{2}}\sum_{k_{x}}e^{ik_{x}Q_{j,y}\ell_{\rm{B}}^{2}}\
\hat{f}^{\dagger}_{1,\alpha;m,k_{x}+Q_{j,x}}\hat{f}_{2,\alpha;n,k_{x}}+\rm{H.c.}$
(49)

###  III.3 Commensurate condition

Omitting layer, sublattice and Laudau level indices, the sum over $k_{x}$ in
Eq. (49) is reminiscent of 1D tight-binding model if we identify in the
following way:

  * •

The position of each repetition motif is $k_{x}\ell_{\rm{B}}^{2}$.

  * •

The lattice spacing is given by $Q^{m}_{x}\ell_{\rm{B}}^{2}>0$, which is the
greatest common divisor of $\\{\forall j,Q_{j,x}\ell_{\rm{B}}^{2}\\}$.

  * •

The nearest hopping is given by an oscillating phase
$e^{ik_{x}Q_{j,y}\ell_{\rm{B}}^{2}}$.

Such 1D chain has lattice translational symmetry if, given $Q_{j,y}$,

$\exists p\in\mathbb{N}\quad\mbox{s.t.}\quad
e^{i(k_{x}+pQ^{m}_{x})Q_{j,y}\ell_{\rm{B}}^{2}}=e^{ik_{x}Q_{j,y}\ell_{\rm{B}}^{2}}\quad\Leftrightarrow\quad\exists(p,q)\in\mathbb{N}^{2}\quad\mbox{s.t.}\quad
Q^{m}_{x}Q_{j,y}\ell_{\rm{B}}^{2}=2\pi\frac{q}{p}.$ (50)

To fulfill this commensurate condition for all $Q_{j,y}$, it suffice that
$Q^{m}_{y}>0$, i.e., the greatest common divisor of $\\{\forall j,Q_{j,y}\\}$,
satisfies Eq. (50). So, we obtain the version of the commensurate condition in
reciprocal space

$Q^{m}_{x}Q^{m}_{y}\ell_{\rm{B}}^{2}=2\pi\frac{q}{p}\quad\mbox{with $p,q$
coprime}.$ (51)

Accordingly, the commensurate condition in real space with
$Q^{m}_{x}=2\pi/\sqrt{3}L_{s}$ and $Q^{m}_{y}=2\pi/L_{s}$ is given by the
magnetic flux $\Phi$ in the Moiré pattern unit-cell of size
$\Omega_{0}=\sqrt{3}L_{s}^{2}/2$

$\frac{\Phi}{\Phi_{0}}=\frac{p}{2q}\quad\mbox{with $p,q$ coprime}$ (52)

where $\Phi_{0}$ is the magnetic flux quantum $h/e$.

Once the commensurate condition is satisfied, the unit-cell of the effective
1D chain in our problem has $p$ motifs. So, we can write

$k_{x}=\tilde{k}_{x}+(sp+r)Q^{m}_{x}$ (53)

with $\tilde{k}_{x}\in[0,Q^{m}_{x}[$ and $(s,r)$ are integers such that
$r=0,\dots,p-1$. Then, we define $Q_{j,x}\equiv\Delta r_{j}Q^{m}_{x}$ and

$\hat{f}_{i,\alpha;n,sp+r}(\tilde{k}_{x})\ \equiv\
\hat{f}_{i,\alpha;n,k_{x}}.$ (54)

So, the sum over $k_{x}$ in Eq. (49) becomes

$\displaystyle\sum_{k_{x}}e^{ik_{x}Q_{j,y}\ell_{\rm{B}}^{2}}\
\hat{f}^{\dagger}_{1,\alpha;m,k_{x}+Q_{j,x}}\hat{f}_{2,\alpha;n,k_{x}}$
$\displaystyle=\sum_{\tilde{k}_{x}}e^{i\tilde{k}_{x}Q_{j,y}\ell_{\rm{B}}^{2}}\sum_{s,r}e^{irQ^{m}_{x}Q_{j,y}\ell_{\rm{B}}^{2}}\
\hat{f}^{\dagger}_{1,\alpha;m,sp+r+\Delta
r_{j}}(\tilde{k}_{x})\hat{f}_{2,\alpha;n,sp+r}(\tilde{k}_{x})$
$\displaystyle=\sum_{\tilde{k}_{x}}e^{i\tilde{k}_{x}Q_{j,y}\ell_{\rm{B}}^{2}}\sum_{r=0}^{p-1}e^{irQ^{m}_{x}Q_{j,y}\ell_{\rm{B}}^{2}}\sum_{\tilde{k}_{y}}\
e^{-i\tilde{k}_{y}\Delta
r_{j}Q^{m}_{x}\ell_{\rm{B}}^{2}}\hat{f}^{\dagger}_{1,\alpha;m,r+\Delta
r_{j}}(\widetilde{\mathbf{k}})\hat{f}_{2,\alpha;n,r}(\widetilde{\mathbf{k}})$

where $\widetilde{\mathbf{k}}=(\tilde{k}_{x},\tilde{k}_{y})$. In the last
line, we use periodic boundary conditions and the following Fourier
transformation

$\displaystyle\hat{f}_{l,\alpha;n,sp+r}(\tilde{k}_{x})$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{\tilde{k}_{y}}\hat{f}_{l,\alpha;n,r}(\widetilde{\mathbf{k}})e^{i\tilde{k}_{y}(sp+r)Q^{m}_{x}\ell_{\rm{B}}^{2}}\quad\mbox{with}\quad\tilde{k}_{y}Q^{m}_{x}\ell_{\rm{B}}^{2}\in[0,\frac{2\pi}{p}[$
(55) $\displaystyle\hat{f}_{l,\alpha;n,r+p}(\widetilde{\mathbf{k}})$
$\displaystyle=\hat{f}_{l,\alpha;n,r}(\widetilde{\mathbf{k}})$ (56)

The moiré potential in this new basis becomes

$\begin{split}\hat{H}^{\mu=+}_{M}=\sum_{j=1}^{3}&\sum_{mn,\alpha\beta}(1+e^{i\phi})\left(T_{j,+}\right)_{\alpha\beta}F^{0}_{mn}(\mathbf{Q}_{j})e^{i\frac{Q_{j,x}Q_{j,y}\ell_{\rm{B}}^{2}}{2}}\\\
&\times\sum_{\widetilde{\mathbf{k}}}\sum_{r=0}^{p-1}e^{i\tilde{k}_{x}Q_{j,y}\ell_{\rm{B}}^{2}}e^{i(rQ_{j,y}-\Delta
r_{j}\tilde{k}_{y})Q^{m}_{x}\ell_{\rm{B}}^{2}}\
\hat{f}^{\dagger}_{1,\alpha;m,r+\Delta
r_{j}}(\widetilde{\mathbf{k}})\hat{f}_{2,\alpha;n,r}(\widetilde{\mathbf{k}})+\rm{H.c.}\end{split}.$
(57)

Therefore, the Dirac Hamiltonian part Eq. (31) becomes

$\begin{split}\hat{H}^{\mu=+}_{D}=&\sum_{l,n,\widetilde{\mathbf{k}}}\sum_{r=0}^{p-1}\hbar
v_{F}(-K_{l,+,x}-iK_{l,+,y})\
\hat{f}^{\dagger}_{l,A;n,r}(\widetilde{\mathbf{k}})\hat{f}_{l,B;n,r}(\widetilde{\mathbf{k}})+\rm{H.c.}\\\
&+\sum_{l,n,\widetilde{\mathbf{k}}}\sum_{r=0}^{p-1}\frac{\sqrt{2(n+1)}\hbar
v_{F}}{\ell_{\rm{B}}}\
\hat{f}^{\dagger}_{l,A;n+1,r}(\widetilde{\mathbf{k}})\hat{f}_{l,B;n,r}(\widetilde{\mathbf{k}})+\rm{H.c.}\end{split}.$
(58)

The calculations for $\mu=-$ is exactly parallel to what we have shown for
$\mu=+$. So, we give directly the results. The moiré potential is

$\begin{split}\hat{H}^{\mu=-}_{M}=\sum_{j=1}^{3}&\sum_{mn,\alpha\beta}(1+e^{i\phi})\left(T_{j,-}\right)_{\alpha\beta}F^{0}_{mn}(-\mathbf{Q}_{j})e^{i\frac{Q_{j,x}Q_{j,y}\ell_{\rm{B}}^{2}}{2}}\\\
&\times\sum_{\widetilde{\mathbf{k}}}\sum_{r=0}^{p-1}e^{-i\tilde{k}_{x}Q_{j,y}\ell_{\rm{B}}^{2}}e^{-i(rQ_{j,y}-\Delta
r_{j}\tilde{k}_{y})Q^{m}_{x}\ell_{\rm{B}}^{2}}\
\hat{f}^{\dagger}_{1,\alpha;m,r+\Delta
r_{j}}(\widetilde{\mathbf{k}})\hat{f}_{2,\alpha;n,r}(\widetilde{\mathbf{k}})+\rm{H.c.}\end{split}$
(59)

and the Dirac Hamiltonian part reads

$\begin{split}\hat{H}^{\mu=-}_{D}=&\sum_{l,n,\widetilde{\mathbf{k}}}\sum_{r=0}^{p-1}\hbar
v_{F}(K_{l,-,x}-iK_{l,-,y})\
\hat{f}^{\dagger}_{l,A;n,r}(\widetilde{\mathbf{k}})\hat{f}_{l,B;n,r}(\widetilde{\mathbf{k}})+\rm{H.c.}\\\
&+\sum_{l,n,\widetilde{\mathbf{k}}}\sum_{r=0}^{p-1}-\frac{\sqrt{2(n+1)}\hbar
v_{F}}{\ell_{\rm{B}}}\
\hat{f}^{\dagger}_{l,B;n+1,r}(\widetilde{\mathbf{k}})\hat{f}_{l,A;n,r}(\widetilde{\mathbf{k}})+\rm{H.c.}\end{split}.$
(60)

###  III.4 Numerical implementation

The commensurate condition gives rise to Hofstadter butterfly in the Landau
level spectrum as a function of magnetic field. In numerical implementation,
we have to introduce a cut-off on the number of Landau levels. As the Dirac
part is (almost) diagonalized in the Landau level basis, the coupling between
Landau levels comes from two places: one is $v_{F}k_{\theta}$ in the diagonal
term and the other is the moiré potential, which has an amplitude
$u_{0}^{\prime}$. So, denoting $E_{c}=\max(v_{F}k_{\theta},u_{0}^{\prime})$, a
viable cut-off needs to encounter the coupling between $n=0$ to $n_{0}$ such
that

$E_{c}\sim\frac{\sqrt{2n_{0}}\hbar v_{F}}{\ell_{\rm{B}}}\Leftrightarrow
n_{0}\sim\frac{\sqrt{3}E_{c}^{2}a^{2}}{8\pi\hbar^{2}v_{F}^{2}\theta^{2}}\cdot\frac{\Phi_{0}}{\Phi}.$
(61)

Note that $E_{c}=u_{0}^{\prime}$ only if $\theta<0.64^{\circ}$. In our plots,
the twist angle is larger than $1^{\circ}$, so $E_{c}=\hbar v_{F}k_{\theta}$
and

$n_{0}\sim\frac{2\sqrt{3}\pi}{9}\cdot\frac{\Phi_{0}}{\Phi}=\frac{1.21}{\frac{\Phi}{\Phi_{0}}},$
(62)

which is independent of $\theta$. Therefore, to achieve the numerical
convergence of the Hofstadter butterfly, one needs a cut-off $N_{c}$, which is
several times larger than $n_{0}$. Since $n_{0}\sim 100$ for
$\Phi/\Phi_{0}=0.01$, we choose in our numerical calculations $N_{c}=500$. For
even smaller $\Phi/\Phi_{0}$, we use $N_{c}=1600$.

Another important point in the numerical calculations is that one should
choose different cut-off, differed by one, for two sublattices. For example,
in valley $\mu=+$, the zeroth Landau level is purely $B$-sublattice polarized
and topologically protected. The higher $n$-th Landau level involves the
$n$-th state of sublattice $B$ and the $(n-1)$-th state of sublattice $B$. If
we choose the same cut-off $N_{c}$ for two sublattices, then the $N_{c}$-th
state of sublattice $A$ has no partner to couple and induces a floating zero
energy state. This would lead to non-physical numerical error hindering the
convergence of Hofstadter butterfly.

###  III.5 Wannier plot

Once we obtain the Laudau levels of a given $k_{z}$-indexed TBG, we can
extract the Chern number of gaps between Landau levels following the Streda
formula. This information is precisely provided by Wannier plots. In this
subsection, we explain how to draw the Wannier plots shown in Fig. 3 of the
main text, inspired from [88].

In Wannier plot, we actually show the density of states $\rho$ as a function
of filling factor $n$ and magnetic field $B$. For illustrative purposes, we
have to broaden the Dirac $\delta$-like peaks in $\rho$ at a given energy $E$
using Lorentz distributions, namely:

$\rho(E)=\sum_{i}\delta(E-E_{i})\to\sum_{i}\frac{1}{\pi}\frac{\gamma}{(E-E_{i})^{2}+\gamma^{2}}$
(63)

where the parameter $\gamma$ mimics Landau level broadening due to thermal or
disorder effects. The filling factor $n$ at energy $E$ is obtained by counting
the number of Landau levels below this energy, namely:

$n(E)=\frac{1}{2q}\sum_{i}\theta(E-E_{i})\to\sum_{i}\frac{1}{2q\pi}\arctan\left(\frac{E-E_{i}}{\gamma}\right)$
(64)

where we add a normalization factor $1/2q$ since the number of magnetic bands
is proportional to $q$, where $p$ and $q$ are co-prime numbers satisfying
$\Phi/\Phi_{0}=p/2q$. In this way, we can achieve the conventional definition
of filling factor, for example, $n\in[-2,2]$ in Wannier plot if one consider
the two flat bands in two valleys only in one spin sector.

Technically, the empirical parameter $\gamma$ has been thus adjusted for each
specific case to achieve best resolution of Wannier plots (usually
$\gamma<10^{-4}$ eV). In our plots, we also rescale the density of states by
its maximum $\rho_{\rm{max}}$. Then, we only plot the points having
$\rho/\rho_{\rm{max}}\leq 1\%$ to focus on the topological nature of large
gaps.

##  IV S4. Hofstadter butterfly spectra

In this section, we provide a detailed analysis on the evolution of Hofstadter
butterfly spectra of $k_{z}$-indexed TBG in bulk ATG of three twist angles
$\theta=2.0^{\circ},\ 2.2^{\circ}$ and $2.4^{\circ}$. These three twist angles
are the representative of the three topological phases, respectively, as shown
in Fig. 9. For each chosen twist angle, we select three representative $k_{z}$
momenta (annotated by green dots in Fig. 9) and account for the Hofstadter
butterfly spectrum of their corresponding $k_{z}$-index TBG. In the following
discussions, we always consider only one spin sector and neglecting Zeeman
effect on Landau levels. The magnetic flux ratio $\Phi_{0}/\Phi$ is defined
with respect to the physical twist angle and the corresponding 2D moiré
superlattice.

Figure 9: Topological phase diagram. The green dots indicate the momenta
$k_{z}$ that we discuss. Their associated twist angle are marked by dashed
lines.

Let us first study the case of $\theta=2.40^{\circ}$. According to the
topological phase diagram Fig. 9, all the $k_{z}$-indexed TBG in bulk ATG is
adiabatically connected to TBG of twist angle larger than the first magic
angle, i.e., there is no magic momentum for $\theta=2.40^{\circ}$. The three
representative $k_{z}$ are $k_{z}d_{0}/\pi=0,1/3,1/2$ (see Fig. 10(a-c),
respectively), which effectively correspond to TBG of twist angle
$\theta=1.2^{\circ},\ 2.4^{\circ}$ and TBG with vanishing moiré potential,
respectively. Here, “effectively” means that $k_{z}$-indexed TBG in bulk ATG
of $\theta=2.40^{\circ}$ share the same dimensionless ration $\alpha$ (defined
in the main text) with 2D TBG of certain twist angle.

For $k_{z}d_{0}/\pi=1/2$, the interlayer moiré coupling is completely
suppressed so that the $k_{z}$-indexed TBG is basically a trivial addition of
two layers of graphene forbidding any interlayer coupling. As shown in Fig.
10(c), Landau levels follow exactly the same $\sqrt{B}$-dispersion as in
graphene (green dashed line). The fractal pattern cannot be seen even for
$\Phi/\Phi_{0}=0.2$ in the given energy scale because the corresponding
magnetic flux in monolayer graphene’s small unit-cell is only $10^{-4}|$.
Because of topological pinning of zeroth Landau level around Dirac cones at
zero energy, the Chern number of the gaps follows the $(2n+1)$-rule in
monolayer graphene if considering only one spin sector. Therefore, the Chern
number of two decoupled graphene monolayers, which both follows the
$(2n+1)$-rule, exhibits the $(4n+2)$-rule, as indicated in Fig. 10(c).

For $k_{z}d_{0}/\pi=1/3$, we effectively study a TBG of twist angle
$2.40^{\circ}$. Although this angle is twice larger than the first magic angle
$\sim 1.1^{\circ}$, the moiré potential already substantially renormalize the
Fermi velocity around Dirac cone. Combining with the large moiré superlattice,
the zeroth Landau level starts to split at $\phi/\phi_{0}=0.15$, i.e.,
$B=20.8$T, as shown in Fig. 10(b). When $B<20.8$T, Landau levels still follow
the $\sqrt{B}$-dispersion except that the Fermi velocity is strongly
suppressed with $v_{F}^{*}/v_{F}\sim 0.4$. So, the Chern number of gaps still
follow the $(4n+2)$-rule, as for $k_{z}d_{0}/\pi=1/2$.

When $k_{z}\to 0$, the effective TBG under study is twisted by $1.2^{\circ}$,
which is close to the first magic angle. The bandwidth and the Fermi velocity
are then radically reduced so that the Landau levels are hard to resolve in
the given energy scale as shown in the Hofstadter spectrum Fig. 10(a). In the
inset, we zoom in on the low-energy Landau levels to energy scale of $\sim 1$
meV. A rich and complex fractal pattern of Landau levels appears at sub Tesla
magnetic field. The cluster of these tightly distributed Landau levels stems
from the two (nearly-)flat bands of TBG. The big gap between this cluster and
other Landau levels turns out to have a vanishing Chern number. The cluster
preserves its energy range while varying magnetic field, as stipulated by
Streda’s formula [101].

The Hofstadter spectrum at low field can be understood by looking into the 2D
band structure at $k_{z}=0$ of bulk ATG, as shown in Fig. 11. The two flat
bands are separated from the four remote bands by a gap at $\Gamma_{s}$ point,
where the remote bands in the same valley are twice degenerate thanks to the
$C_{3}$ and $M_{y}$ symmetries [82]. Most saliently, while extrapolating the
gap to zero magnetic field, the separation between the first high-energy
Landau level and the low-energy cluster is precisely the gap between the flat
bands and the two high-energy bands at $\Gamma_{s}$ as shown in Fig. 11.
Therefore, the Landau levels aside the big gap with zero Chern number have to
stem from the $\Gamma_{s}$ point of the flat bands and the remote bands. At
the vicinity of $\Gamma_{s}$, the six bands in study can be described by a
six-band $k\cdot p$ Hamiltonian [82], which is reminiscent of that found in
semiconductors with spin-orbit coupling. The off-diagonal part of the
Hamiltonian attributes non-zero Berry curvature to the bands. Therefore, the
Landau levels of these bands cannot naively classified into topologically
trivial Schrödinger fermion with quadratic dispersion or conventional massive
Dirac fermion with $\pi$-Berry phase (e.g., that found in hBN). What we have
here is something in between such that:

  * •

Landau levels from the remote bands follow the $\sqrt{B}$ dispersion while the
zeroth Landau level does not keep constant but vary with magnetic field.

  * •

The Landau levels from $\Gamma_{s}$ of the flat bands follows the linear $B$
dispersion.

Furthermore, the Chern number of gaps follows the $2n$-rule. This can be
understood by looking into Fig. 11. The factor 2 comes from the valley
degeneracy. Although the first high energy band is twice degenerate at
$\Gamma_{s}$, their distinct dispersion away from $\Gamma_{s}$ induces
different Landau level splitting so that the increment of the Chern number is
2 instead of 4.

Figure 10: Hofstadter butterfly spectra for $\theta=2.4^{\circ}$ at
$k_{z}d_{0}/\pi=0$, 1/3 and 1/2. Landau level dispersion in monolayer graphene
is marked with green dashed lines. The Chern number of gap are indicated by
red. Figure 11: Band structure zoomed at $\Gamma_{s}$ for
$\theta=2.40^{\circ}$ at $k_{z}=0$.

In the case of $\theta=2.20^{\circ}$ and $2.00^{\circ}$, the Hofstadter
spectra at $k_{z}d_{0}/\pi=1/2$ [see Fig. 12(c) and Fig. 14(c), respectively]
are exactly the same as for $\theta=2.40^{\circ}$ because they correspond to
the same decoupled bilayer graphene with vanishing interlayer coupling.

When $\theta=2.20^{\circ}$, we choose two other $k_{z}$ as $k_{z}d_{0}/\pi=0$
and 0.10. According to the topological phase diagram, the former is in the
$C=0$ region and the latter is on the phase boundary. The Landau levels of the
flat bands are zoomed in on the energy scale, shown in the inset of Fig.
12(a,b). When $k_{z}d_{0}/\pi=0$, the Chern number and the energy dispersion
of Landau levels originated from remote bands have the same behavior as for
$\theta=2.40^{\circ}$ at $k_{z}=0$ for the same reasons. However, in contrast
to $\theta=2.40^{\circ}$, the energy of Landau level of the flat band edges
has the $\sqrt{B}$ dispersion, which may be seen as a signature of topological
transition by exchange of Berry curvature at band touching points [see Fig.
13(left)]. When $k_{z}d_{0}/\pi=0.10$, gap reopens at $\Gamma_{s}$ but the
Landau levels are hard to resolve because of completely flattened bands [see
Fig. 13(right)] close to the topological transition boundary.

Figure 12: Hofstadter butterfly spectra for $\theta=2.2^{\circ}$ at
$k_{z}d_{0}/\pi=0$, 0.10 and 1/2. Landau level dispersion in monolayer
graphene is marked with green dashed lines. The Chern number of gap are
indicated by red. Figure 13: Band structure zoomed at $\Gamma_{s}$ for
$\theta=2.20^{\circ}$ at $k_{z}d_{0}/\pi=0$ and 0.10.

When $\theta=2.00^{\circ}$, we choose two other $kz$ as $k_{z}d_{0}/\pi=0$ and
0.13. According to the topological phase diagram, the former is in the $C=1$
region but smaller than the first magic momentum and the latter is in the
small region of $C=0$. The Landau levels of the flat bands are zoomed in on
the energy scale, shown in the inset of Fig. 14(a,b). It turns out that the
Hofstadter spectra Fig. 14(b) and the band structure Fig. 15 for
$\theta=2.00^{\circ}$ at $k_{z}d_{0}/\pi=0.13$ are similar to the case for
$\theta=2.20^{\circ}$ at $k_{z}d_{0}/\pi=0$. Therefore, the understanding of
them should be also the same since they belongs to the same topological phase.
However, the case for $\theta=2.00^{\circ}$ at $k_{z}d_{0}/\pi=0$ is
different. Although the band structure for $k_{z}d_{0}/\pi=0$ (left panel of
Fig. 15) is similar to $k_{z}d_{0}/\pi=0.13$ (right panel of Fig. 15), they
belongs to distinct topological phase, according to the phase diagram Fig. 9.
This can be seen from the zeroth Landau level of the remote band in Fig.
14(a), which follows a dispersion linear in $B$ instead of $\sqrt{B}$, i.e., a
hint of change in Berry curvature.

Figure 14: Hofstadter butterfly spectra for $\theta=2.0^{\circ}$ at
$k_{z}d_{0}/\pi=0$, 0.13 and 1/2. Landau level dispersion in monolayer
graphene is marked with green dashed lines. The Chern number of gap are
indicated by red. Figure 15: Band structure zoomed at $\Gamma_{s}$ for
$\theta=2.40^{\circ}$ at $k_{z}d_{0}/\pi=0$ and 0.13.

##  V S5. Bound states of the spiral dislocation line

In the main text, we have studied the physical properties of twisted graphene
spiral in regions far from the dislocation line, which forms a bulk
alternating twisted graphite. We now investigate the properties of the
dislocation line. We start with the lattice structure of twisted graphene
spiral. The dislocation line is located at the $AA$ point in the moiré
supercell. We consider the lattice distortion in the $z$ direction at position
r, while neglecting the in-plane displacement, i.e. $u_{x}=u_{y}=0$. The
lattice distortion in the $z$ direction at a lateral position r (by setting
the dislocation line as the origin) can be expressed as:

$u_{z}(\textbf{r})=\frac{2d_{0}}{2\pi}\arg(r_{x}+\,i\,r_{y}),$ (65)

where the function arg denotes the argument of the complex number. The
presence of the dislocation line breaks the moiré superlattice translational
symmetry. Thus, we set up an open boundary condition in the $x$-$y$ direction
and a periodic boundary in the $z$ direction. We perform the lattice
relaxation calculation utilizing Large-scale Atomic-Molecular Massively
Parallel Simulation (LAMMPS) [96], employing long-rang bond-order potential
for carbon (LCBOP) [97]. The initial structure is an $n\times\,n$ ($n=2,4,6$)
supercell with the dislocation line at the center of the sample. Based on the
relaxed lattice structure, we calculate the band structure of twisted graphene
spiral including the screw dislocation line using the atomistic tight binding
model described by Eq. (19). Then, we project the energy band to the nearest
neighbor carbon atoms around the dislocation line, in order to evaluate the
single-particle spectral function of the bounded states associated with the
spiral dislocation line. In Fig. 16, we present the spectral function of the
bounded states near the dislocation line for twisted graphene spiral with
different twist angles. Specifically, We construct a $2\times 2$ supercell for
twisted graphene spiral with twist angle $\theta=1.05^{\circ}$, a $2\times 2$
supercell for GS with $\theta=1.47^{\circ}$, a $4$ supercell with
$\theta=2.28^{\circ}$ and a $6\times 6$ supercell with $\theta=3.48^{\circ}$.
The dashed lines mark the energy of the Dirac point. We observe the presence
of several bounded states that cross the energy level of the Dirac point.
These bound states do not show any qualitative difference at different twist
angles.

Figure 16: The spectral function of the bounded states near the spiral
dislocation line for twisted graphene spiral with (a) $2\times 2$ supercell
and twist angle $\theta=1.05^{\circ}$; (b) $2\times 2$ supercell and twist
angle $\theta=1.47^{\circ}$; (c) $4\times 4$ supercell and twist angle
$\theta=2.28^{\circ}$; and (d) $6\times 6$ supercell and twist angle
$\theta=3.48^{\circ}$.