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

# Moiré Fractional Chern Insulators I: First-principles calculations and
Continuum Models of Twisted Bilayer MoTe2

Yujin Jia  Beijing National Laboratory for Condensed Matter Physics and
Institute of physics, Chinese academy of sciences, Beijing 100190, China
University of Chinese academy of sciences, Beijing 100049, China  Jiabin Yu
Department of Physics, Princeton University, Princeton, New Jersey 08544, USA
Jiaxuan Liu  Beijing National Laboratory for Condensed Matter Physics and
Institute of physics, Chinese academy of sciences, Beijing 100190, China
University of Chinese academy of sciences, Beijing 100049, China  Jonah
Herzog-Arbeitman  Department of Physics, Princeton University, Princeton, New
Jersey 08544, USA  Ziyue Qi  Beijing National Laboratory for Condensed Matter
Physics and Institute of physics, Chinese academy of sciences, Beijing 100190,
China  University of Chinese academy of sciences, Beijing 100049, China
Nicolas Regnault  Department of Physics, Princeton University, Princeton, New
Jersey 08544, USA  Laboratoire de Physique de l’Ecole normale supérieure, ENS,
Université PSL, CNRS, Sorbonne Université, Université Paris-Diderot, Sorbonne
Paris Cité, 75005 Paris, France  Hongming Weng  Beijing National Laboratory
for Condensed Matter Physics and Institute of physics, Chinese academy of
sciences, Beijing 100190, China  University of Chinese academy of sciences,
Beijing 100049, China  Songshan Lake Materials Laboratory, Dongguan, Guangdong
523808, China  B. Andrei Bernevig<EMAIL_ADDRESS>Department of
Physics, Princeton University, Princeton, New Jersey 08544, USA  Donostia
International Physics Center, P. Manuel de Lardizabal 4, 20018 Donostia-San
Sebastian, Spain  IKERBASQUE, Basque Foundation for Science, Bilbao, Spain
Quansheng Wu<EMAIL_ADDRESS>Beijing National Laboratory for
Condensed Matter Physics and Institute of physics, Chinese academy of
sciences, Beijing 100190, China  University of Chinese academy of sciences,
Beijing 100049, China

###### Abstract

Recent experiments observed fractional Chern insulators (FCI) in twisted
bilayer MoTe2 at zero magnetic field, yet even the single-particle model of
this material is controversial, leading to unreliable predictions of the
experimental phase diagram as discussed in [Yu et al., 2023]. In this light,
we revisit the single-particle model of twisted bilayer MoTe2. Utilizing
large-scale density functional theory, we calculate the band structure of
twisted AA-stacked bilayer MoTe2 at various twist angles relevant to
experiment. We find that a band inversion occurs near $4.41^{\circ}$ between
the second and third bands. Our _ab initio_ band structure is in qualitative
agreement with [Wang et al., 2023], but shows important differences in the
remote bands and in the $\Gamma$ valley. We incorporate two higher harmonic
terms into the continuum model to capture the highest 3 valence bands per
valley. We confirm that the two highest valence bands per valley have opposite
Chern numbers with $|C|=1$ for small angles, and also use our model to predict
a variety of Chern states in the remote bands accessible by displacement
field. Finally, we perform DFT calculations and build models for the AB
stacking configuration. Our work serves as a foundation for accurate
determination of the correlated phases in twisted bilayer MoTe2.

##  I Introduction

Fractional Chern insulators (FCI) [1, 2, 3] were show to naturally appear in
zero magnetic field when nearly flat Chern bands [4, 5] are fractionally
filled. Over the last several years, there have been extensive theoretical [6,
7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26,
27, 28, 29, 30, 31, 32] and experimental [33, 34, 35, 36, 37] studies on zero-
field FCIs in moiré materials [38, 39], as well as studies of fractional
quantum Hall (FQH)-like states under nonzero external magnetic field [40]
(which were called FCIs because the Chern bands where they appeared did not
have flat Berry curvature and quantum geometry), and small B-field-induced
FCIs [41]. Remarkably, FCIs without any external magnetic fields have recently
been observed in twisted bilayer MoTe2 ($t$MoTe2) [33, 34, 35, 36] and in the
pentalayer-graphene/hBN moiré superlattice [42]. These _moiré fractional Chern
insulator_ (mFCI) states provide the best experimental platform to date for
the physics proposed in Refs. [1, 2, 3].

In this work, we will focus on $t$MoTe2, where mFCIs are observed at
fractional fillings $\nu=-2/3,-3/5$, as well as a Chern insulator (CI) at
$\nu=-1$. Throughout, $\nu$ is the electron filling measured from the charge
neutrality point. FCIs are expected to appear [1, 2, 3] in the CI bands found
in the $t$MoTe2 model first proposed in Ref. [43]. However, an accurate model
of the material is required to correctly reproduce the FCI states and
competing non-topological states. Otherwise serious disagreement with
experiment can occur, as has been somewhat overlooked in the flurry of recent
literature. For example, the spin polarization at $\nu=-1/3,-4/3$ cannot be
reproduced theoretically unless band-mixing is taken into account [31],
showing that a bands beyond the lowest valence manifold play an integral part
in the physics. Recent theoretical work [20, 21, 23, 24, 26, 28, 29, 31, 32,
44], some of which predicts the existence of fractional states not seen in the
experiment, might then be subject to change when band-mixing is properly
included. Even interaction-driven band mixing beyond the two-band Hilbert
space considered in Ref. [31] could be important if the single-particle bands
are energetically close. Hence it is of the utmost importance to obtain a
continuum model capturing the ab-initio band structure over the range of
energies accessible by the Coulomb potential.

Most importantly, there are two different sets of ab-initio parameters Ref.
[20, 21] used [20, 21, 23, 24, 26, 28, 29, 31, 32, 44] for interacting
calculations in the continuum model of Ref. [43]. These sets of parameters
give rise to different single particle phase diagrams, many-body spin
polarizations, and stabilities of 2/3 and 1/3 mFCI states Ref. [31]. Settling
the single-particle model of $t$MoTe2 is the subject of this first work in the
mFCI series.

We perform large-scale density functional theory (DFT) calculations to study
twisted bilayer MoTe2 at various twist angles for AA stacking, which is the
stacking configuration relevant to the experiments [33, 34, 35, 36]. To
accurately capture the impact of crystal structure relaxation, we rigorously
test 19 different van der Waals (vdW) exchange-correlation functionals, and
determine that DFT-D2 yields the most reliable lattice parameters when
compared with experimental results. With DFT-D2, we further use a highly
efficient two-step method to obtain the structure. Our first step is to use
the DFT-D2 functional to train a machine learning algorithm that generates
relaxed structures for a set of twist angles. The second step is to use these
relaxed structures as the initial configurations, and further perform the full
DFT relaxation to obtain final, accurate relaxed structures. The moiré bands
are eventually calculated with the final relaxed structures. Our DFT results
show that the $\pm\text{K}$-valley valence bands dominate the low-energy
physics while the $\Gamma$-valley valence bands are about 80meV away from the
valence band maximum (VBM) for the experimentally relevant angle
$3.89^{\circ}$. Furthermore, the second and the third valence bands in K
valley (or $-\text{K}$ valley related by time-reversal symmetry) undergo a gap
closing around $4.41^{\circ}$, changing the Chern number of the second valence
band from $-1$ at smaller angles to $1$ at larger angles. Our ab-initio
results for the $\pm\text{K}$ valley bands are closer to those of Ref. [21]
than to those of Ref. [20], although our second and third band per valley are
closer to each other than those of Ref. [21] at $3.89^{\circ}$, which shifts
their gap closing to smaller angles. Our $\Gamma$-valley band is also lower in
energy than that of Ref. [21], and hence will not contribute to many-body
physics near $\nu=-1$. The band structure in Ref. [20] looks different from
ours; we can reproduce similar band width of the topmost valence moiré bands
in Ref. [20] by using ISMEAR=1 parameters in the settings of VASP
calculations. However, in our calculations all through this work we use
ISMEAR=0 [45], which leads to the narrow bandwidth.

We then use the moiré model to capture the low-energy DFT bands. If we only
keep the first harmonics (FH) as in Ref. [43], we can only manage to capture
the top two valence bands in each valley. To capture the top three valence
bands in each valley, we add two extra second harmonic (SH) terms, and obtain
a good match with the DFT band structure and symmetry representations (reps).
We propose that this more accurate model be used in many-body calculations. In
our model with SH terms, the gap closing between the second and third top
valence bands in one valley happens around $4.2^{\circ}$ at zero displacement
field, which is consistent with our DFT calculation. We find that adding a
displacement field can change the Chern number of the top valence band from
$1$ to $0$ in K valley, and can also achieve a variety of Chern numbers (from
-2 to 2) for the second and third top valence bands in one valley. Accessing
these bands provides another route to integer Chern physics seen in twisted
transition metal dichalcogenides[46, 47, 48, 49, 50].

In addition to the experimentally relevant AA-stacking configuration, we also
study the AB-stacking configuration. Our DFT results show that the
$\Gamma$-valley bands in the AB stacking case is closer to the VBM (only about
30meV away) and are extremely flat. We build a FH model capturing these
$\Gamma$-valley bands, and show that they are extremely localized atomic bands
whose flatness comes from zero hopping among atomic orbitals on the triangular
lattice. We also used the $\pm\text{K}$-valley model in Ref. [43] to match the
top two valence bands in each valley.

In the rest of this paper, we discuss the DFT calculations at a range of twist
angles and stacking configurations in Sec. II, and the continuum models we
employ to faithfully reproduce these calculations in Sec. III. We conclude the
paper in Sec. IV, and provide more details in a series of appendices.

##  II DFT Calculations

In this section, we discuss the large-scale DFT calculations on the $t$MoTe2
at various twist angles and different stacking configurations. We will mainly
discuss the results for 3.89∘, which is the closest one to the structure in
recent experiments [33, 34, 35, 36], as well as the topological phase
transition from large angle to 3.89∘. A complete discussion can be found in
Appendix. A.

###  II.1 3.89∘ $t$MoTe2

Utilizing the coincidence lattice method [51, 52], we construct twisted
bilayer crystal structures of MoTe2 at various commensurate angles: 13.2∘,
9.43∘, 7.34∘, 5.09∘, 4.41∘, 3.89∘ and 3.48∘, considering both AA and AB
stacking. Here, AA (AB) stacking implies that, without any twist, the Mo/Te
atoms of the top layer respectively align with the Mo/Te (Te/Mo) atoms of the
bottom layer. Subsequently, large-scale DFT calculations using the Vienna _ab
initio_ simulation package (VASP) [53, 54, 55, 56, 57] are performed on these
structures.

Figure 1: Relaxation and band structure of 3.89∘ $t$MoTe2. $(a)$ Relaxation of
the 3.89∘ AA stacking. $(b)$ Relaxation of the 3.89∘ AB stacking. $(c)$ and
$(d)$ show the band structure and density of states (DOS) of the structures in
$(a)$ and $(b)$ respectively, with red lines indicating bands from the
$\Gamma$ valley and the green line marking the Fermi level.The notations
$k_{v1}$, $k_{v2}$, and $k_{v3}$ represent the three highest pairs of valence
bands of the AA-stacked configuration in the $\pm\text{K}$ valleys. The DOS is
normalized per molybdenum atom.

To capture the van der Waals (vdW) interactions between top and bottom layers,
we test 19 exchange-correlation functionals using the experimental bulk
crystal structure [58] as a benchmark (see Tab. 3 in Appendix. A.1).
Ultimately, we find that DFT-D2 gives the lattice parameter closest to the
experimental value [58], and we use this functional throughout the work.
Secondly, different pseudo-potential (PP) combinations have also been tested.
It has been found that the energy differences between different PPs are
negligible (see Fig. 7 of Appendix. A.1). Considering the computational cost,
we choose the PAW pseudo-potential and PBE exchange-correlation functional.

Based on these funtionals, we develop a highly efficient two-step relaxation
scheme which combines machine learning and DFT. First, we construct the
machine-learned force field (MLFF) using the relaxation data generated by
DFT-D2. The MLFF method has a much lower computational cost than the direct
relaxation with DFT since the moiré unit cells are very large with 1302 aotms.
MLFF can produce a relaxed structure quickly. Using this structure as an
initial guess, we further perform the full DFT relaxation. It turns out that
the MLFF is reasonably good, making the DFT relaxation quite fast. To achieve
total-force convergence with an accuracy of 5meV/Å and energy convergence of
$1\times 10^{-5}$eV, the DFT+MLFF method requires 17 ionic steps which take
7.5 hours using 16 NVIDIA A100 GPUs. In contrast, direct relaxation using DFT
requires 178 ionic steps and take 55 hours when using 8 NVIDIA A100 GPUs,
demanding approximately four times more resources than the DFT+MLFF method.

The faster convergence of the DFT+MLFF method (compared to the DFT direct
relaxation) does not sacrifice precision in the band structure. To show it, we
present the band structures obtained from three different approaches: MLFF
relaxation, DFT+MLFF relaxation, and DFT direct relaxation, as shown in Fig.
17 of Appendix. A.3.4. The band structures obtained from the DFT+MLFF and DFT
relaxed methods have qualitatively the same shapes, which is consistent with
the fact that the DFT relaxation upon the MLFF structures is quite fast. In
particular, the bands from the DFT+MLFF and the direct DFT relaxation methods
are extremely similar, verifying the validity of the DFT+MLFF method.
Therefore, we use the DFT+MLFF method in this work instead of the DFT direct
relaxation. The stacking-dependent corrugated moiré structures of AA and AB
configuration generated by the DFT+MLFF method are shown in Fig. 1(a,b).

AA-stacked $t$MoTe2 has a two-fold rotational symmetry axis along $y$ axis,
$C_{2y}$, and the $C_{3}$ symmetry. The relaxed structure exhibits a maximum
interlayer distance $d=7.62\mbox{\text{\AA}}$ in the AA region, where the
metal atoms in top layer is aligned with metal atoms in bottom layer, while a
minimum interlayer distance $d=7.0\mbox{\text{\AA}}$ in the MX (XM) region,
where the top layer metal (chalcogen) atoms are aligned with chalcogen (metal)
atoms of bottom layer, as shown in Fig. 6 of Appendix. A.1.

For AA $t$MoTe2 at $3.89$ degree, the VBM is located at the $\text{K}_{M}$
point in the moiré Brillouin zone (BZ), which is folded from the K point in
the untwisted bilayer structure. The top three pairs of valence bands, labeled
as $k_{v1}$, $k_{v2}$ and $k_{v3}$, originate from the $\pm\text{K}$ valleys,
exhibiting bandwidths of 12.8 meV, 16.2 meV, and 16.5 meV, respectively. The
combined effects of lattice relaxation and SOC lead to the $\Gamma$-valley
bands shifting downward by about 80 meV from the VBM. The $\Gamma$-valley
bands, illustrated in Fig. 1 (c) and marked by red dashed lines, contribute to
two distinct peaks in the density of states (DOS).

Figure 2: Irreducible representations at high-symmetry points $\Gamma$ and $K$
of valence bands of monolayer and twist MoTe2 with twist angle 5.09∘, 4.41∘,
and 3.89∘ for relaxed AA stacking configuration with SOC. The eigenvalues of
the $C_{3}$ are denoted as $\omega=e^{i\pi/3}$, $\omega^{*}=e^{-i\pi/3}$ and
$\bar{1}=e^{-i\pi}$. The critical point of band inversion between the 2nd and
3rd valence bands at $\gamma$ point of moire BZ happens at around twist angle
4.41∘.

The AB stacking configuration has a two-fold rotational symmetry $C_{2x}$ with
the axis along $x$ as well as the three-fold rotation symmetry $C_{3}$ with
axis along $z$. For the twist angle $3.89^{\circ}$, the top two pairs of
valence bands also come from $\pm\text{K}$ valley with band width of 16 meV,
while the first pair of ultra flat bands from $\Gamma$ valley is only 14.6 meV
below the VBM. The charge density of these ultra flat bands is highly
localized in the AB region of the moiré lattice (see Fig.20 of Appendix.
A.3.4).

###  II.2 Topological phase of valence bands from 5.09∘ to 3.89∘

In this part, we will discuss the topological phase transition around
$4.41^{\circ}$ of the valence bands in AA-stacked $t$MoTe2. Since the full DFT
calculations are performed in huge unit cells, it is impractical to calculate
the Chern number directly, and we turn to symmetry eigenvalues to efficiently
deduce the topology[59].

Fig. 2 shows the band structure and the $C_{3}$ eigenvalues at the high-
symmetry points. Our convention is to use the spinful $C_{3}$ eigenvalues,
which are labelled by $\omega=e^{i\pi/3}$, $\omega^{*}=e^{-i\pi/3}$ and
$\bar{1}=e^{-i\pi}$. Since $\pm\text{K}$-valley bands are valley-spin locked,
we distinguish two valleys by their spins (spin up/down for K/$-\text{K}$
valley).

For 5.09∘, we observe that there is a gap of 2.8 meV between $k_{v2}$ and
$k_{v3}$ bands at the moiré $\Gamma_{M}$ point. This gap is closed around
4.41∘, and reopens at angle 3.89∘. From the $C_{3}$ eigenvalues labeled in
Fig. 2, we see that this band crossing exchanges the symmetry representations
at $\Gamma_{M}$ for $k_{v2}$ and $k_{v3}$ bands. Explicitly, for 5.09∘, the
$k_{v2}$ band at $\Gamma_{M}$ has two spin-polarized states:
$\left|{\uparrow,\omega^{*}}\right>$ and $\left|{\downarrow,\omega}\right>$,
while the $k_{v3}$ band has $\left|{\uparrow,\omega}\right>$ and
$\left|{\downarrow,\omega^{*}}\right>$. When energy gap reopening occurs, the
symmetry representations of $k_{v2}$ and $k_{v3}$ at 3.89∘ switch in
comparison to those at 5.09∘.

The exchange of symmetry eigenvalues is proof of a band inversion and causes a
change in the spin/valley Chern number. Specifically, recall that the Chern
number can be determined from symmetry via
$e^{i\frac{2\pi}{3}C}=-\xi_{\Gamma_{M}}\xi_{\text{K}_{M}}\xi_{\text{K}^{\prime}_{M}}$
where $\xi_{\mathbf{k}}$ is the spinful $C_{3}$ eigenvalue at high-symmetry
point $\mathbf{k}$ [59]. Since $C_{2y}\mathcal{T}$ relates the moiré
$\text{K}_{M}$ and $\text{K}_{M}^{\prime}$ points and is anti-unitary,
$\xi_{\text{K}_{M}}=\xi_{\text{K}_{M}^{\prime}}$. Thus for $\theta$ slightly
larger than $4.41^{\circ}$, we find $e^{i\frac{2\pi}{3}C}=-(-1)(-1)\omega^{*}$
so that $C=1$ mod 3 for the second top spin-$\uparrow$ band, but for $\theta$
slightly smaller than $4.41^{\circ}$, we find
$e^{i\frac{2\pi}{3}C}=-(-1)(-1)\omega$ so that $C=-1$ mod 3. This topological
phase transition is matched by the continuum model as we will show in Sec.
III.

###  II.3 Comparison with Ref. [20, 21]

Two recent works [20, 21] have also studied the relaxation and band structure
of $t$MoTe2. The band structures in Ref. [20, 21] are different from our
result, mainly due to the different relaxed structures. Ref. [20] employed the
SCAN density functional with dDsC dispersion correction to perform crystal
structure relaxation using VASP, while Ref. [21] used SIESTA with DFT-D2
functional to perform the relaxation. The comparison of the relaxation results
of AA-$t$MoTe2 between this work and Ref. [20, 21] is shown in Fig. 3. As
shown in Fig. 3, the interlayer distance in our relaxed structure has the
qualitatively the same shape as those in Ref. [20, 21]—largest interlayer
distance at MM and the smallest interlayer distance at MX/XM. However, the
interlayer distance in our relaxed structure has smaller spatial fluctuations
than that in Ref. [20] as shown in Fig. 3(a,b), while our relaxed structure
has larger interlayer distance than that of Ref. [21] (Fig. 3(c)). The maximum
interlayer distances of this work, Ref. [20], and Ref. [21] are about 7.6Å,
7.8Å and 7.4Å, respectively, and the minimum interlayer distances are
respectively about 7.0Å, 7.0Å and 6.9Å.

Our relaxation result is consistent with the AA and AB stacking untwisted
bilayer structure. In the AA region, the stacking configuration is close to
that of AA untwisted bilayer structure, and thus the maximum interlayer
distance should be close to but slightly smaller than (due to the corrugation
effect due to the connection to other stacking configurations in the moiré
structure) the interlayer distance of AA stacking untwisted bilayer (7.7Å),
which is consistent with our results but not with Ref. [20]. Furthermore, the
MX region has the stacking configuration akin to that of AB untwisted bilayer
structure. As a result, the smallest interlayer distance in $t$MoTe2
structures should be close to but slightly larger than (due to the
corrugation) 7Å, which is consistent with our results but not with Ref. [21].
More details about the relaxation and its influence on band structure are
discussed in Appendix. A.3.5.

Figure 3: Relaxation results of AA twisted structure. (a) Relaxation result of
4.41∘ AA MoTe2 from Fig.2 (a) of Ref. [20]. (b) Our relaxation result of 4.41
∘ AA MoTe2. (c) The interlayer distance of 3.89∘ AA MoTe2 along the black
arrow in Fig. 1(a). The red curve is the result from Fig.1(d) from Ref. [21],
and the blue curve is our result. Figure 4:  The comparison between the
valence bands from the DFT calculation (black dots) and from the continuum
model (red/orange line) for twist angle $3.89^{\circ}$ for (a,b) AA-stacking
and (c) AB-stacking along the high-symmetry line. The red line comes from the
$\pm\text{K}$-valley model, while the orange line comes from the
$\Gamma$-valley model. Each red line is doubly degenerate (except along
$\Gamma_{M}-M_{M}$ in (c) where it is 4-fold degenerate). The orange line in
(c) has double degeneracy around $-30$meV and 4-fold degeneracy around
$-86$meV. The fitting in (a) is done with only FH terms, while SH terms are
added in (b).

##  III Continuum Models

In this section, we use the continuum model to fit the DFT results. Up to now,
FCI states were only found for hole doping experimentally [33, 34, 35, 36];
thus, we will focus on the model for the valence bands. We will first discuss
the AA-stacking and then discuss the AB stacking.

###  III.1 AA-Stacking

According to the DFT results, the low-energy valence bands mainly originate
from the $\pm\text{K}$ valleys in the monolayer. The symmetry group of AA-
stacking $t$MoTe2 (AA-$t$MoTe2) is generated by $C_{3}$, $C_{2y}$, and
$\mathcal{T}$, in addition to the moiré lattice translations [43]. We pick a
convention where the top layer is rotated by $-\theta/2$ and the bottom layer
by $\theta/2$. The moiré lattice constant is

$a_{M}=\frac{a_{0}}{2\sin\left(\frac{\theta}{2}\right)}\ ,$ (1)

where $a_{0}=3.52\text{\AA}$ is the lattice constant of monolayer MoTe2.

The continuum model for the monolayer $\pm\text{K}$ valleys in general reads

$H_{\eta,0}^{AA}=\int
d^{2}r\left(c^{\dagger}_{\eta,b,\bm{r}},c^{\dagger}_{\eta,t,\bm{r}}\right)h_{\eta,0}^{AA}(\bm{r})\left(\begin{matrix}c_{\eta,b,\bm{r}}\\\
c_{\eta,t,\bm{r}}\end{matrix}\right)\ ,$ (2)

where $c^{\dagger}_{\eta,l,\bm{r}}$ labels the basis of the continuum model,
$\eta=\pm$ labels the $\pm\text{K}$ valleys (or equivalently spins), $l=t,b$
labels the layer, and $\bm{r}$ labels the position. Ref. [43] proposed a model
with only the first harmonics (FH); however, here we add the certain second
harmonics (SH) terms in order to accurately match the higher bands. (See
Appendix. B.1 for details.) Expanding the potential terms and using symmetry,
we find the form

$\displaystyle\begin{split}&V_{\eta,l}(\bm{r})=Ve^{-(-)^{l}\mathrm{i}\psi}\sum_{i=1,2,3}e^{\mathrm{i}\bm{g}_{i}\cdot\bm{r}}+Ve^{(-)^{l}\mathrm{i}\psi}\sum_{i=1,2,3}e^{-\mathrm{i}\bm{g}_{i}\cdot\bm{r}}\\\
&\qquad+2V_{2}\sum_{i=1}^{3}\cos(\bm{g}_{2i}\cdot\bm{r})\\\
&t_{\eta}(\bm{r})=w\sum_{i=1,2,3}e^{-\eta\mathrm{i}\bm{q}_{i}\cdot\bm{r}}+w_{2}\sum_{i=1,2,3}e^{-\eta\mathrm{i}\bm{q}_{2i}\cdot\bm{r}}\
,\end{split}$ (3)

where we take $(-)^{t}=1$, $(-)^{b}=-1$, $\bm{g}_{i}=C_{3}^{i-1}\bm{b}_{M,1}$,
$\bm{b}_{M,1}=\frac{4\pi}{\sqrt{3}a_{M}}(1,0)^{T}$ and
$\bm{b}_{M,2}=\frac{4\pi}{\sqrt{3}a_{M}}$ are the basis moiré reciprocal
lattice vectors,
$\bm{q}_{1}=\frac{4\pi}{3a_{0}}2\sin\left(\frac{\theta}{2}\right)(0,1)^{T}$,
$\bm{q}_{2}=C_{3}\bm{q}_{1}$, $\bm{q}_{3}=C_{3}^{2}\bm{q}_{1}$,
$\bm{g}_{21}=\bm{b}_{M,1}+\bm{b}_{M,2}$, $\bm{g}_{2i}=C_{3}^{i-1}\bm{g}_{21}$,
$\bm{q}_{21}=\bm{b}_{M,1}+\bm{q}_{{}_{1}}$, and
$\bm{q}_{2i}=C_{3}^{i-1}\bm{q}_{21}$. $V$, $\psi$ and $w$ characterize the FH
terms, while $V_{2}$ and $w_{2}$ belong to the SH. The AA-stacking FH
$\pm\text{K}$-valley model (Eq. (2)) has effective inversion symmetry that
makes the two bands from the two valleys identical (in accord with the DFT
results which show small splitting about 1.2 meV The effective inversion
symmetry is natural with only FH terms, and we only include the SH terms that
preserve the effective inversion symmetry. In total, the model (Eq. (2)) has
six real parameters $m^{*}$, $V$, $\psi$, $w$, $V_{2}$ and $w_{2}$.

We fit to the DFT band structure at $\theta=3.89^{\circ}$ in two ways. (See
Appendix. B.1 for details.) First, we set $V_{2}=w_{2}=0$, which corresponds
to the FH model. In this case, we manage to fit the top 4 valence bands (2 in
each valley) with the corresponding FH parameters in Tab. 1, as shown in Fig.
4(a). Then, we allow nonzero $V_{2}$ and $w_{2}$, _i.e._ , adding the SH
terms. We are now able to fit the top 6 valence bands (3 in each valley) with
the corresponding FH+SH parameters in Tab. 1, as shown in Fig. 4 (b). The
match is not only good along the high-symmetry line but also good in the full
BZ as shown in Fig. 21 (a,b) in Appendix. B.3.

As shown in Fig. 22 (a,b) in Appendix. B.3, we can see that the $C_{3}$
eigenvalues for the top 6 valence bands match the DFT calculation in both FH
and FH+SH cases. Furthermore, the Chern numbers of the top 3 bands (in
decreasing order of energy) in K valley are (1,-1,0) in both cases, which are
consistent with the $C_{3}$ eigenvalues [59]. The $C_{3}$ eigenvalues and the
Chern numbers of the top two bands per valley are the same as those in Ref.
[21]. It is clear that adding the SH terms improves the reliability of the
model across a wider range of energies. We expect more remote bands to be
accessible in future experiments, and hence the our FH+SH model is an
essential improvement.

Model |  $m^{*}$ ($m_{e}$) | $V$ | $\psi(^{\circ})$ | $w$ | $V_{2}$ | $w_{2}$
---|---|---|---|---|---|---
FH | 0.60 | 16.5 | -105.9 | -18.8 | 0 | 0
FH+SH | 0.62 | 7.94 | -88.43 | -10.77 | 20.00 | 10.21
Table 1: Values of the parameters in the $\pm\text{K}$-valley continuum model
(Eq. (LABEL:main_eq:Kvalley_V_t_FH_SH)) for the AA-stacking $t$MoTe2.
$V,w,V_{2},w_{2}$ are in meV. “FH” means we only include the first harmonics,
whereas “FH+SH” means that we include both the first harmonics and the
effective-inversion-invariant second harmonics.

At last, we discuss the evolution of the bands of the AA-stacking moiré model
in Eq. (2) with the FH+SH parameters values in Tab. 1 as a funciton of the
twist angle $\theta$ and the displacement field $\varepsilon$. We will focus
on the K valley. As shown in Fig. 4 (a), the gap between second top and third
top bands closes around $4.2^{\circ}$ for zero displacement field, which is
close the the DFT’s $4.41^{\circ}$ in Fig. 2; the gap closing will change the
Chern numbers of the second and third top bands in K valley from (-1,0) to
(1,-2). Further increasing the angle at zero displacement field will cause a
band inversion between the third and fourth top valence bands, which changes
the Chern number of the third top band from -2 to -1 as shown in Fig. 4 (d).
Fig. 4 (b) shows that increasing the displacement field can trivialize the top
valence band, while a variety of Chern numbers (ranging from -2 to 2) can
arise for nonzero displacement field for the second and third top bands as
shown in Fig. 4 (c,d).

Figure 5:  Phase diagram of AA-stacking moiré model in Eq. (2) with the FH+SH
parameters values in Tab. 1. In this plot, band 1,2,3 refer to the top, second
top and third top valence bands in the K valley. In (a), we show the Chern
numbers of the three bands (blue for 1, red for 2 and orange for 3) for zero
displacement field, which shows a band inversion between the second and third
top bands around $4.2^{\circ}$. In (b,c,d), we show the direct gaps and the
Chern numbers of the top (a), second top (b) and third top (c) bands in the
$\pm\text{K}$ valley. Explicitly, the direct gaps are illustrated by the color
(the color bar is the same across plots) and the Chern numbers are labeled in
white.

###  III.2 AB-Stacking

The generators of the symmetry group of AB-stacking $t$MoTe2 (AB-$t$MoTe2),
which can be thought of as twisting the top layer of AA-$t$MoTe2 by another
$180^{\circ}$, are $C_{3},C_{2x}$, and $\mathcal{T}$. Note that $C_{2x}$ is
local to the monolayer K point (unlike $C_{2y}$ in AA-$t$MoTe2), and thus it
preserves the valley quantum number in the moiré model. This difference in the
valley symmetry group is important and, as we now show, leads to different
behavior with the potential for interesting many-body spin and Hubbard physics
[60, 61, 62, 63, 64, 65, 66].

The DFT results show that the low-energy valence bands of AB-$t$MoTe2 come
from both the $\pm\text{K}$ valleys and the $\Gamma$ valleys in the monolayer.
The $\pm\text{K}$-valley model was proposed in Ref. [43], which reads

$\displaystyle\begin{split}H_{\eta,0}^{AB}&=\sum_{l}\int
d^{2}rc^{\dagger}_{\eta,l,\bm{r}}\left[\frac{\hbar^{2}\nabla^{2}}{2m^{*}}+V_{\eta,l}(\bm{r})+(-)^{l}\frac{\varepsilon}{2}\right]c_{\eta,l,\bm{r}}\\\
&\ +\left[\int
d^{2}rc^{\dagger}_{\eta,b,\bm{r}}t_{\eta}(\bm{r})c_{\eta,t,\bm{r}}+h.c.\right]\
,\end{split}$ (4)

where $c^{\dagger}_{\eta,l,\bm{r}}$ labels the basis of the continuum model
with $\eta=\pm$ labeling the $\pm\text{K}$ valleys (or equivalently spins),
$l=t,b$ labels the layer, and $\bm{r}$ labels the position. Since AB-$t$MoTe2
is given by rotating the top layer of AA-$t$MoTe2 by an extra $180^{\circ}$,
the two layers in K or $-\text{K}$ valley now have opposite spin. With first
harmonics, the forms of $V_{\eta,l}(\bm{r})$ and $t_{\eta}(\bm{r})$ derived
from the symmetries read

$\displaystyle\begin{split}&V_{\eta,l}(\bm{r})=Ve^{-\mathrm{i}\psi}\sum_{i=1,2,3}e^{\mathrm{i}\bm{g}_{i}\cdot\bm{r}}+Ve^{\mathrm{i}\psi}\sum_{i=1,2,3}e^{-\mathrm{i}\bm{g}_{i}\cdot\bm{r}}\\\
&t_{\eta}(\bm{r})=w\sum_{i=1,2,3}e^{\eta\mathrm{i}(i-1)\frac{2\pi}{3}}e^{-\eta\mathrm{i}\bm{q}_{i}\cdot\bm{r}}\
,\end{split}$ (5)

where $V$ is real, and $w$ can made real by choosing the relative phase
between the two layers. For $\varepsilon=0$, the minimal continuum model has
the effective TR symmetry within each valley that flips layer (or equivalently
spin), making the bands from the two valleys identical.

On the other hand, the $\Gamma$-valley model has the following form

$H_{\Gamma}=\int
d^{2}r\left(\begin{matrix}\psi^{\dagger}_{\bm{r},b}&\psi^{\dagger}_{\bm{r},t}\end{matrix}\right)h_{\Gamma}^{AB}(\bm{r})\left(\begin{matrix}\psi^{\dagger}_{\bm{r},b}\\\
\psi^{\dagger}_{\bm{r},t}\end{matrix}\right)\ ,$ (6)

where we choose the kinetic term in the continuum model as the intralayer
spin-independent $\nabla^{2}$ term,

$\displaystyle\begin{split}&h_{\Gamma}^{AB}(\bm{r})=\frac{\hbar^{2}\nabla^{2}}{2m_{\Gamma}^{*}}+E_{\Gamma}+\left(\begin{matrix}V_{\Gamma,b}(\bm{r})-\epsilon/2&t_{\Gamma}(\bm{r})\\\
t_{\Gamma}^{\dagger}(\bm{r})&V_{\Gamma,t}(\bm{r})+\epsilon/2\end{matrix}\right)\
,\end{split}$ (7)
$\psi^{\dagger}_{\bm{r},l}=(\psi^{\dagger}_{\bm{r},l,\uparrow},\psi^{\dagger}_{\bm{r},l,\downarrow})\
,$ (8)

$t$ and $b$ correspond to the top and bottom layers, respectively,
$E_{\Gamma}$ accounts for the energy difference between the $\Gamma$-valley
and $\pm\text{K}$-valley bands, and $V_{\Gamma,l}(\bm{r})$ and
$t_{\Gamma}(\bm{r})$ are $2\times 2$ matrix functions. For zero displacement
field, the $\Gamma$-valley model has effective TR symmetry with in each spin-
subspace, which makes the bands from the two spins identical (resulting in at
least double degeneracy of each band).

We fit the low-energy bands at $\theta=3.89^{\circ}$, and the resulting
parameter values are summarized in Tab. 2. As the illustration, we show the
good match between the DFT bands and those from the models along the high-
symmetry line in Fig. 4 (b). The bands match well also in the full BZ, and the
$C_{3}$ eigenvalues also match the DFT calculation as discussed in Appendix.
C.3.

For the $\pm\text{K}$-valley model, we notice that the interlayer coupling can
be set to zero $w=0$ in the $\pm\text{K}$-valley model while keeping the match
of the bands good, which can be understood as the follows. The two layers in
one valley now have opposite spins; owing to the spin U(1) symmetry for the
low-energy states near $\pm\text{K}$ valleys in monolayer MoTe2, we expect the
spin $U(1)$ symmetry is approximately preserved in $t$MoTe2, which means the
interlayer coupling is very small for the AB-stacking. The zero interlayer
coupling makes the eigenstates have well-defined valley and layer/spin. As a
result, the two states with the same spin in the $\pm\text{K}$-valley model
are degenerate at $\text{K}_{M}$, since the combination of the effective TR
symmetry and the TR symmetry leaves the spin invariant; this is consistent
with DFT calcualtions as shown in Fig. 19 of Appendix. C.3. The top four
valence bands is in the $A_{1}@1a$ atomic limit, and the band structure can be
thought of as arising from hopping on this moiré triangular lattice.

The $\Gamma$-valley bands are extremely flat as shown in Fig. 4 (b), owing to
the large $m^{*}_{\Gamma}=10m_{e}$
($-\frac{\hbar^{2}|\bm{q}_{1}|^{2}}{2m^{*}_{\Gamma}}=-2.49$meV) compared to
the potential $V_{\Gamma}=72$meV. As a result, those flat bands are extremely
localized atomic states localized around the minima of the intralayer
potential (_i.e._ , 1a positions due to $\psi=0$), and their energies can be
approximately calculated from an array of decoupled harmonic oscillators
related by moiré translations. (See Appendix. C.3.)

At last, we discuss the effect of the displacement field. As shown in Fig. 25
in Appendix. C.3, the effect of the displacement field on the
$\pm\text{K}$-valley bands are just to shift the bands from different layers
relative to each other, since the valley is a good quantum number due to the
zero interlayer coupling (see Tab. 2). On the other hand, the effect of the
displacement field on the low-energy $\Gamma$-valley bands is negligible,
which is consistent with the fact that the very large interlayer coupling
makes the eigenstates equally distributed between the two layers (see Tab. 2).

$m^{*}$ ($m_{e}$) |  $V$ (meV) | $\psi(^{\circ})$ |  $w$ (meV)
---|---|---|---
0.62 | 53 | -56 | 0
$m_{\Gamma}^{*}$ ($m_{e}$) |  $V_{\Gamma}$ (meV) | $\psi_{\Gamma}(^{\circ})$ |  $w_{\Gamma}$ (meV)
10 | 72 | 0 | 300
Table 2: Values of the parameters in the $\pm\text{K}$-valley continuum model
in Eq. (LABEL:main_eq:Kvalley_model_AB) (first and second rows) and the
$\Gamma$-valley continuum model in Eq. (LABEL:main_eq:Gammavalley_model_AB)
(third and fourth rows) for the AB-stacking $t$MoTe2.

##  IV Conclusions

Our extensive DFT study, accelerated by machine learning, has confirmed that,
and our continuum model analysis show that the lowest two bands in the K
valley have Chern numbers $+1$ and $-1$ respectively, at a twist angle of
$3.89^{\circ}$ for the experimentally relevant AA stacking. These Chern
numbers mod 3 are accessible directly from the DFT data, and we compute the
exact value from the continuum model fit to the DFT bands. The Chern numbers
are consistent with Ref. [21] but do not agree with those computed by Ref.
[20] at larger twist angles and extrapolated to $3.89^{\circ}$.

Unlike the fitting in Ref. [20, 21] which captures the top two valence bands
per valley, we succeeded in matching the top three valence bands (as well as
much of the fourth band) per valley by adding only two higher harmonic terms
to the moiré Hamiltonian (see Eq. (2)). This Hamiltonian still preserves the
effective intra-valley inversion symmetry, and will serve as a faithful model
of the dispersion and topology for a wide range of electron fillings.
Excitingly, our model predicts a rich topological phase diagram accessible
through displacement fields in the remote bands. Our forthcoming work will
study the many-body physics of this model, contributing to the study of
correlations and topology [67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 50] now
accessible in experiment.

##  V Acknowledgments

The authors thank Allan H. MacDonald, Kin-Fai Mak, Jie Shan, Zhijun Wang and
Xiaodong Xu for helpful discussions, with special thanks to D. Xiao, Y. Zhang
and Jianpeng Liu for their helpful cross-checking of the details in the DFT
calculations. J.Y. and J. H.-A. are grateful for conversations with Pok Man
Tam. This work was supported by the Ministry of Science and Technology of
China (Grant No. 2022YFA1403800), the Science Center of the National Natural
Science Foundation of China (Grant No. 12188101) and the National Natural
Science Foundation of China (Grant No.12274436). H.W. acknowledge support from
the Informatization Plan of the Chinese Academy of Sciences
(CASWX2021SF-0102). B. A. B.’s work was primarily supported by the DOE Grant
No. DE-SC0016239 and the Simons Investigator Grant No. 404513. N.R. also
acknowledges support from the QuantERA II Programme that has received funding
from the European Union’s Horizon 2020 research and innovation programme under
Grant Agreement No 101017733 and from the European Research Council (ERC)
under the European Union’s Horizon 2020 Research and Innovation Programme
(Grant Agreement No. 101020833). . O.V. was funded by the Gordon and Betty
Moore Foundation’s EPiQS Initiative Grant GBMF11070, National High Magnetic
Field Laboratory through NSF Grant No. DMR-1157490 and the State of Florida.
J. H.-A. is supported by a Hertz Fellowship, with additional support from DOE
Grant No. DE-SC0016239 by the Gordon and Betty Moore Foundation through Grant
No. GBMF8685 towards the Princeton theory program, the Gordon and Betty Moore
Foundation’s EPiQS Initiative (Grant No. GBMF11070), Office of Naval Research
(ONR Grant No. N00014-20-1-2303), BSF Israel US foundation No. 2018226 and
NSF-MERSEC DMR. J. Y. is supported by the Gordon and Betty Moore Foundation
through Grant No. GBMF8685 towards the Princeton theory program.

##  Appendix A DFT Results

###  A.1 Atomic structures of $t$MoTe2

Figure 6: Rigid atomic structures of 13.2∘ $t$MoTe2. The purple one is Mo atom
and yellow-brown one is Te atom. (a) Bilayer AA stacking configuration. (b)
Bilayer AB stacking configuration. (c) 13.2∘ $t$MoTe2 AA configuration with
in-plane 2-fold rotational symmetry axis along the y-axis. (b) 13.2∘ $t$MoTe2
AB configuration with in-plane 2-fold rotational symmetry axis along the
x-axis. Both AA and AB have the $C_{3}$ 3-fold rotation symmetry axis
perpendicular to the plane. Mo atom is labeled by M, while chalcogen atom is
labeled by X. MX represents the top layer Mo atom is aligned with bottom layer
chalcogen atom. The similar with XM and MM.

Bulk 2H-MoTe2, as sketched in Fig. 6 (b), has a hexagonal structure with space
group 194 generated by the inversion $\mathcal{P}$, a three-fold symmetry
$C_{3}$ with axis along z axis, a two-fold symmetry $C_{2x}$ with axis along x
axis, and lattice translations. Experimentally, the lattice constant is found
to be a=3.519Å and c=13.976Å [58].

Before performing DFT calculations for band structure, we need to choose the
suitable exchange-correlation functionals to describe the electron-electron
interaction and the suitable pseudopotentials to deal with the interactions
between electrons and nucleus. To pick the best exchange-correlation
functionals, we try 19 different ones to calculate the relaxed bulk crystal
structure, and show the results in Tab. 3. We find that the DFT-D2 functional
(IVDW=10) provides a lattice constant of $a=b=3.518\mbox{\text{\AA}}$ and
$c=13.976\mbox{\text{\AA}}$ for the bulk MoTe2 primitive crystal structure,
which is closest to the experimental values in Ref. [58]. Therefore, we
conclude that the DFT-D2 functional (IVDW=10) is the most appropriate choice
for the MoTe2 structural relaxation among the 19 functionals.

| experiment | optB86 | optB88 | vdw-DF | vdw-DF-cx
---|---|---|---|---|---
a,b(Å) | 3.519 | 3.527 | 3.567 | 3.631 | 3.502
c(Å) | 13.964 | 14.032 | 14.213 | 15.007 | 13.867
| PBE | optPBE-vdw | rVV10 | SCAN+rVV10 | r2SCAN+rVV10
a,b(Å) | 3.551 | 3.580 | 3.546 | 3.503 | 3.542
c(Å) | 15.095 | 14.475 | 13.949 | 14.223 | 14.187
| rev-vdw-DF2 | vdw-DF2 | IVDW=10 | IVDW=11 | IVDW=12
a,b(Å) | 3.529 | 3.711 | 3.519 | 3.512 | 3.490
c(Å) | 14.028 | 14.891 | 13.976 | 13.984 | 13.649
| IVDW=20 | IVDW=21 | IVDW=263 | IVDW=4 | IVDW=3
a,b(Å) | 3.514 | 3.516 | 3.490 | 3.515 | 3.531
c(Å) | 13.923 | 13.826 | 13.709 | 14.047 | 14.222
Table 3: Relaxed lattice constant of bulk MoTe2 using different vdw functionals. The IVDW number corresponds to different vdw functionals provided in VASP. Here IVDW=10 is the DFT-D2 method of Grimme (marked by gray). IVDW=11 is the DFT-D3 method of Grimme with zero-damping function. IVDW=13 is the DFT-D4 method. IVDW=20 is the Tkatchenko-Scheffler method. IVDW=21 is the Tkatchenko-Scheffler method with iterative Hirshfeld partitioning. IVDW=263 is the Many-body dispersion energy with fractionally ionic model for polarizability method. IVDW=4 is the dDsC dispersion correction method. IVDW=3 is the DFT-ulg method. Experimental data is from [58].  | experiment | Mo-Te | Mo_pv-Te | Mo_sv-Te
---|---|---|---|---
a,b(Å) | 3.519 | 3.519 | 3.521 | 3.523
c(Å) | 13.964 | 13.976 | 13.985 | 13.997
Table 4: Relaxed lattice constant of bulk MoTe2 using different
pseudopotentials. Three pseudopotentials combinations show little difference.
Experimental data is from [58].  Figure 7: Band structure of bulk, monolayer,
AA bilayer and AB bilayer MoTe2 calculated with different pseudopotentials.
The crystal structures are relaxed with corresponding pseudopotentials and
IVDW=10(DFT-D2) functional.

After selecting the exchange-correlation functionals, we test the different
pseudopotentials—Mo-Te, Mo_pv-Te, Mo_sv-Te—where the pseudopotentials without
suffices are general pseudopotential, and the suffix ‘_sv’ (‘_pv’) means that
the inner s (p) electrons are considered as valence electrons. These
combinations show little difference in the relaxed lattice parameters, and is
consistent with experimental data. In addition, the band structure calculated
by different pseudopotentials coincides with each other in different MoTe2
structures. Therefore, it is legitimate to choose any of them. Considering the
balance of accuracy and computational cost, here we choose the "Mo-Te"
combination.

Before studing the electronic structures of the $t$MoTe2 system, we first
generated their crystal structures. We start from both the AA and AB
stacking[77]. As illustrated in Fig. 6, AA stacking means that, when there is
no twist, the Mo/Te atoms of the top layer aligns with the Mo/Te atoms of the
bottom layer. AB stacking means that, when there is no twist, the Mo/Te atoms
of the top layer is directly above the Te/Mo atoms of the bottom layer (see
Fig. 6 (a)). The bilayer untwisted crystal structure relaxed by DFT-D2
functional gives lattice parameter $a=3.5228\mbox{\text{\AA}}$, and thus two
primitive lattice vectors are $a_{1}=a(1,0,0)$ and
$a_{2}=\frac{a}{2}(-\frac{1}{2},\frac{\sqrt{3}}{2}$,0). For twisted
homobilayer with both top and bottom layer being MoTe2, a commensurate
structure occurs when the moiré lattice vector of top layer and bottom layer
satisfy the commensurate lattice condition
$\bm{a}_{Mb}=n_{1}\bm{a}_{1}+n_{2}\bm{a}_{2}=m_{1}\bm{a}_{1}^{r}+m_{2}\bm{a}_{2}^{r}=\bm{a}_{Mt}$
for certain integers $n_{1},n_{2},m_{1},m_{2}$, where $\bm{a}_{1}^{r}$ is the
primitive lattice vectors rotated by an angle $\theta$. In this way, we can
obtain the rigid $t$MoTe2 structures with different angles. All the rigid
structures were generated using a homemade software 2DTwist.

The moiré structures twisted from AA stacking has a hexagonal structure in
space group 150 , with $C_{3}$, $C_{2y}$ (see Fig. 6 (a)) as well as the time-
reversal (TR) symmetry. Thus, the little group in the K valley is generated by
$C_{3},C_{2y}\mathcal{T}$, and $\pm\text{K}$ valleys are exchanged by $C_{2y}$
and $\mathcal{T}$. In contrast, Fig. 6 (b) shows the moiré structures twisted
from AB stacking which has a hexagonal structure of space group 149, with
$C_{3}$, $C_{2x}$ and $\mathcal{T}$. For the AB configuration, the little
group in the K valley is generated by $C_{3}$ and $C_{2x}$, and the two
valleys are exchanged by $\mathcal{T}$.

As the twist angle becomes smaller but nonzero, a $\mathrm{U}(1)$ valley
symmetry emerges due to the exponential suppression of inter-valley scattering
off the moiré potential. This will enable us to build continuum models around
the monolayer $\pm\text{K}$ points using the little group symmetries to
constrain the low-order terms.

###  A.2 Relaxation of $t$MoTe2

Since the relaxation will greatly affect the band structure, it is necessary
to perform relaxation on $t$MoTe2. However, there are 1302 atoms for 3.89∘ and
1626 atoms for 3.48∘ in the moire cell of twisted structures, making the
relaxation process difficult to converge. We decided to construct a Machine
Learning Force Field (MLFF) to get relaxed structures in an efficient way.
MLFF is a machine learning algorithm that will “learn” energies and forces of
atoms from ab-initio calculation and can be applied to predict forces and
energies for similar systems. We note that MLFF is not a relaxed structure—it
is a function that maps structures to forces/energies; it can be used to
efficiently generate the relaxed structures. We firstly constructed a MLFF and
applied it to get an MLFF-relaxed structure. Then, we performed DFT relaxation
on the MLFF-relaxed structures. It only took around 20 DFT steps in relaxing
the largest moire structure to converge in this strategy. In comparison, the
relaxation from rigid 3.89∘ AA-$t$MoTe2 takes 178 steps to converge.

####  A.2.1 Construction of Machine Learning Force field

During the construction of MLFF, two software packages are used. One is the
VASP together with its integrated MLFF module [78], and the second is
NequIP[79], which is an MLFF built on an E(3) equivariant neural network.

We use the VASP MLFF module to generate the _ab-initio_ data needed for
training a high-precision MLFF. The VASP MLFF module itself is a way to
generate MLFF and accelerate the Molecular Dynamics (MD) simulation. It
firstly runs several _ab initio_ MD steps, collects all the energies, forces
and structure data into a dataset, then trains a MLFF, and finally estimates
the error by doing MLFF-based MD calculations. Here the MLFF-based MD
calculation is based on Bayesian linear regression, and the error is directly
estimated by the spread of the Gaussian distribution. If the estimated error
of MLFF is large, an additional _ab initio_ MD step will be performed to
enlarge dataset, the MLFF is retrained with the new dataset, and estimate the
error again. The procedure will be repeated until the total number of MD steps
exceeds the preset.

However, the MLFF algorithm in VASP is lightweight. It helps accelerate the MD
simulation, but the MLFF generated by the algorithm is not accurate enough.
Therefore, we will not directly use the generated MLFF; instead, we just run
the VASP MLFF module for tens of thousands of steps to generate a set of _ab-
initio_ MD data.

With the set of _ab-initio_ MD data generated from the VASP MLFF module, we
use the NequIP software to train an accurate MLFF. NequIP is based on an E(3)
equivariant neural network, meaning that the input and output of each neural
network layer are equivariant (in other words, covariant) under the rotation,
reflection and translation in 3D space. The NequIP software is reported to
outperform several other MLFF algorithms in both data efficiency and accuracy
[79].

The training process is summarized as the follows. We started from small
supercells of untwisted AA and AB bilayer MoTe2 with different in-plane shift
between top and bottom layer. We ran MD simulation using VASP MLFF module on
those structures and collected the data from all ab-initio MD steps. We merged
all the collected data and trained a NequIP MLFF, which is then used in Atomic
Simulation Environment [80] to relax the $t$MoTe2 at various angles.

####  A.2.2 Relaxation results

The relaxation results of AA and AB MoTe2 in different twist angles are listed
in Fig. 8, Fig. 9 and Fig. 10. As we mentioned before, the DFT+MLFF-relaxed
results are obtained by further DFT relaxation based on the MLFF-relaxed
structures. In Fig. 8 and Fig. 9, the MLFF-relaxed structures are close to the
DFT+MLFF-relaxed structure. The extra DFT relaxation in the DFT+MLFF method
only modified the quantitative details, keeping the qualitative shape of the
MLFF-relaxed structures unchanged. It indicates that the MLFF can reproduce
the main part of DFT relaxation, making the DFT relaxation easier to converge.

Figure 8: Relaxation results of 7.34∘ AA and AB $t$MoTe2. (a), (b), (c) and
(d) are interlayer distances of MLFF-relaxed AA structure, DFT+MLFF-relaxed AA
structure, MLFF-relaxed AB structure and DFT+MLFF-relaxed AB structure
respectively. (e), (f), (g) and (h) are intralayer displacements of MLFF-
relaxed AA structure, DFT+MLFF-relaxed AA structure, MLFF-relaxed AB structure
and DFT+MLFF-relaxed AB structure. Interlayer distance is the distance between
the top and bottom layer, while the intralayer displacement indicates the in-
plane displacement from rigid positions to relaxed positions of a Mo atom in
top layer. In the interlayer distance plots, rigid structures are selected to
have the same lattice constant as relaxed structures.

In structures of smaller twist angles, the local conformation is more similar
to untwisted structures. Because our MLFF are constructed from untwisted
structures, more accurate results in smaller angles are expected. Compare Fig.
9 with Fig. 8, the MLFF-relaxed results are indeed better at 3.89∘ than that
at $7.34^{\circ}$.

Figure 9: Relaxation results of 3.98∘ AA and AB $t$MoTe2. (a), (b), (c) and
(d) are interlayer distances of MLFF-relaxed AA structure, DFT+MLFF-relaxed AA
structure, MLFF-relaxed AB structure and DFT+MLFF-relaxed AB structure
respectively. (e), (f), (g) and (h) are intralayer displacements of MLFF-
relaxed AA structure, DFT+MLFF-relaxed AA structure, MLFF-relaxed AB structure
and DFT+MLFF-relaxed AB structure respectively.

Fig. 10 shows the relaxation results of AA and AB $t$MoTe2 in different twist
angles. In both AA and AB stacking, interlayer distance comes to the lowest
point in MX/XM/AB region and becomes higher in MM/XX/AA region. In AA
stacking, the interlayer distance of MX region is the same as XM region since
MX and XM configurations are related by $C_{2}$ symmetry. In the case of AB
stacking, however, the MM region is different from XX region, as they are not
symmetry-related. The interlayer distance of XX region is always higher than
MM region, which results in the asymmetry in Fig. 10 (b) and Fig. 10 (d).

Figure 10: Relaxation results of different angles along black arrows in Fig. 9
(b, d). (a) and (b), interlayer distance of AA structures given by MLFF and
DFT+MLFF, respectively. (c) and (d), interlayer distance of AB structures
given by MLFF and DFT+MLFF.

###  A.3 Electronic structures of $t$MoTe2

In this part, we discuss the electron band structures of $t$MoTe2 obtained
from the DFT+MLFF relaxed structure.

####  A.3.1 Energy bands of monolayer, AA bilayer and AB bilayer MoTe2

Figure 11: Band structure of monolayer, AB bilayer and AA bilayer MoTe2 with
and without relaxation. (a), band structure of monolayer, the effective mass
of the monolayer MoTe2 of the VBM at the $\Gamma$ and K points, and the CBM at
the K point, along two different directions, are indicated in the orange boxes
with units $m_{e}$. (b), band structure of AB bilayer, where the rigid and
relaxed) bands of AA and AB bilayer MoTe2 are indicated by black and red
lines, respectively. (c, band structure of AA bilayer. All the effective
masses are calculated for relaxed structures with SOC.

Before diving into electronic structure of the $t$MoTe2, we first study the
monolayer and bilayer structures. For the monolayer structure, the inclusion
of SOC brings a clear splitting for top valence band at K point. Additionally,
At the $\Gamma$ point, SOC causes a downward shift of the top valence band by
70 meV in comparison to the scenario without SOC, as depicted in Fig. 11 (a).

Fig. 11 (b) shows that lattice relaxation has negligible influence on the band
structure of AB stacking. The reason is that the functional we used is
optimized for the bulk MoTe2 which is also AB stacked. Conversely, Fig. 11 (c)
demonstrates that lattice relaxation significantly alters the band structure
of AA stacking. This is due to the fact that the relaxed interlayer distance
measures approximately 7.7 Å, which is around 10% greater than the interlayer
distance of the rigid structure at 7 Å. Relaxation causes a downward shift of
the first valence band at the $\Gamma$ point, resulting in an energy about
380meV lower at $\Gamma$ compared to the K point. For conduction bands, upon
relaxation, the bottom bands along $\Gamma$-K path are elevated to a similar
energy level as those at the K point.

####  A.3.2 13.2∘ $t$MoTe2

We start from the $13.2^{\circ}$ $t$MoTe2 and gradually decrease the twist
angle to follow the evolution of the band structure. In Fig. 12 (a), when spin
orbit coupling (SOC) is not considered, there is an isolated narrow valence
band (NVB) located at the top of the valence bands for AA stacking
$13.2^{\circ}$ $t$MoTe2 without relaxation. A very similar band is also
present for the AB stacking with negligible difference (Fig. 12 (d)). After
considering the SOC, the NVB start to go down and entangle with other bands.
(See Fig. 12 (b,e).) Further including the relaxations, we obtain stacking-
dependent corrugated moiré structures. Relaxation pushes the NVB further below
as shown in Fig. 12 (c,f).

By projecting the moiré bands into atoms’ orbitals as shown in Fig. 12
(l)-(n), it is clear that the NVB is consist of Mo atoms $d_{z^{2}}$ orbitals
which is the same band character of the VBM at $\Gamma$ point of monolayer
MoTe2 (see Fig. 14 (a)). So we can say that the NVB is the $\Gamma$ valley
band. Using the orbital nature, we can see that upon considering both SOC and
lattice relaxation effects, the $\Gamma$-valley bands are pushed down by 250
meV below the valence band maximum (VBM) for AA configuration, while about 180
meV for AB configuration. Then, the state around the VBM are composed of
molybdenum (Mo) $d_{x^{2}-y^{2}}$ and $d_{xy}$ orbitals which have the same
orbital characteristics as the VBM at the $\pm\text{K}$ points in a monolayer
MoTe2 (see Fig. 14 (a)). Therefore, the VBM of the moiré energy bands
primarily originates from the $\pm\text{K}$ valley. We note that the downward
shift of the $\Gamma$ bands, caused by the SOC and lattice relaxation effects,
is consistent with the lower $\Gamma$-valley bands after including the
relaxation in the AA-stacking untwisted case in Fig. 11(c).

Figure 12: Band structure and band character analysis (fatband) of 13.2∘
$t$MoTe2. (a) is the band structure for rigid AA without SOC, showing
negligible difference with rigid AB without SOC in (d). Considering SOC, the
top isolated valence band tangles with other valence bands both for rigid AA
and AB. SOC brings larger splitting in $\Gamma-M$ k-path for rigid AA in (b)
while larger splitting around $K$ point for rigid AB in (e). After the
relaxation, the bands from $\pm\text{K}$ valley are lifted up to the top of
valence bands, both for relaxed AA in (c) and relaxed AB in (f).

####  A.3.3 AA stacking: Evolution of band structure from 9.43∘ to 3.89∘

We also calculate the band structures of several relaxed AA stacking $t$MoTe2
structure at 9.43∘, 7.34∘, 5.09∘, 4.41∘, 3.89∘ and 3.48∘ as shown in Fig. 13.

We discuss the valence bands first. The valence band maximum is at the
$\text{K}_{M}$ points of the moire BZ. The top pair of valence bands remains
quasi degenerate in $\Gamma_{M}-K_{M}-M_{M}$ high symmetry line, suggesting
the presence of an extra symmetry than TR symmetry, which would map $-\bm{k}$
to $\bm{k}$. As shown in Fig. 14 on the orbital natures of the bands at 9.43∘,
the valence bands near the VBM mainly come from $\pm\text{K}$ valleys, owing
to their $d_{x^{2}-y^{2}}$ and $d_{xy}$ nature, similar to the discussion of
Appendix. A.3.2; the expectation is the bands labeled by red dashed lines in
Fig. 13, which comes form the $\Gamma$ valley owing to its $d_{z^{2}}$ orbital
nature.

Analyzing the evolution of the valence bands in Fig. 13 from panel (g) to (l),
it becomes evident that the bandwidth of the top two pairs of valence bands
decreases. narrows as the twist angle decreases—–a characteristic commonly
observed in twisted systems. Moreover, there is a discernible indication of
band inversion between $k_{v2}$ and $k_{v3}$ at $\Gamma_{M}$, occurring as the
twist angle decreases from 5.09∘ and 3.48∘. Specifically, at the $\Gamma_{M}$
point for 5.09∘, there is a gap of approximately 2.8 meV separating the second
and third pairs of top valence bands. This gap closes at a twist angle of
4.41∘ and reopens to about 1.2 meV at 3.89∘. This cycle of gap closing and
reopening suggests the possibility of a band inversion. To investigate this,
we calculated the irreducible representations of the six highest valence bands
at the $\text{K}_{M}$ and $\Gamma_{M}$ points for twist angles ranging from
7.34∘ to 3.89∘ and have presented the findings in Fig. 2.

Figure 13: Evolution of band structures of AA stacking configuration $t$MoTe2
with twist angles ranging from 9.43∘ to 3.48∘ . The six band structures are
calculated from DFT+MLFF relaxed moiré structure with the consideration of SOC
effects. The dashed red lines show the flat bands coming from $\Gamma$ valley.
Figure 14: Orbital-projected band structure of monolayer layer MoTe2 and AA
(AB) $t$MoTe2 with twisted angle 9.43∘. Bands are projected on $d_{xy}$,
$d_{x^{2}-y^{2}}$ and $d_{z^{2}}$ of Mo atoms. Figure obtained with the open-
source code PYPROCAR[81].

As shown in Fig. 13 from panel (a) to (f), the conduction bands are quite
messy. Based on the decomposition of orbital content depicted in Fig. 14, we
determine that the conduction band minimum (CBM) at the $\Gamma_{M}$ point for
a 9.43° twist exhibits the same orbital characteristics—specifically Mo
$d_{x}y$ and $d_{x^{2}-y^{2}}$—as those found in the electron pocket along the
$\Gamma$-K path in the monolayer band structure, also illustrated in Fig. 14.
The conduction band minimum is around the $\Gamma_{M}$ point rather than
$\text{K}_{M}$ point can be understood from the AB and AA band structures in
Fig. 11.

The presence of the conduction band minimum around the $\Gamma_{M}$ point, as
opposed to the $\text{K}_{M}$ point, can be elucidated by examining the band
structures for AB and AA configurations in Fig. 11. In the AB bilayer band
structure, the dip along the $\Gamma-K$ path is lower than at the K point
itself. Conversely, in the AA bilayer, the $\Gamma-K$ path’s minimum aligns
with the level at the K point. Given that the moiré unit cell’s stacking
transitions between AA and AB types, the moiré conduction band minimum
consequently shifts away from the $\text{K}_{M}$ point.

We discuss the erratic evolution of the electron bands in Appendix. C.4, which
we can attribute to the lowest energy states coming from an electron-like
pocket. This pocket occurs at generic momentum $\mathbf{p}_{i}$ in the
untwisted BZ, so that the bands are folded around a generic point
$\mathbf{p}_{i}\,\mathrm{mod}\,\mathbf{G}_{M}$ in the moiré BZ that depends
extremely sensitively on the angle through $\mathbf{G}_{M}$.

####  A.3.4 3.89∘ AA stacking $t$MoTe2 electronic structures

Let us focus on the 3.89° twisted AA stacking. In Fig. 15, we present the band
structure calculations for both the rigid and relaxed configurations, with and
without SOC. Before relaxation, the $\Gamma$ valley ultra flat bands are
located at the top of the valence bands when SOC is omitted. The introduction
of SOC raises some dispersive bands; however, even with SOC, the rigid
structure maintains flat $\Gamma$ valley bands at the top. Relaxation
dependent on stacking has a significant impact on the band structure. When SOC
is absent, the relaxation raises the $\pm\text{K}$ bands. Further
incorporating SOC, $\pm\text{K}$ bands new becomes the VBM, and the red dashed
bands originating from the $\Gamma$ valley, as denoted in Fig. 15 (d), are
markedly diminished, sinking approximately 80 meV below the VBM. Our results
are different from Ref. [21], where the $\Gamma$ bands are close to the top
two valence bands from $\pm\text{K}$ valley, with only 30 meV away from the
VBM. The top valence band in Ref. [21] has the bandwidth of about 9 meV, while
the width of the same band in our results is about 12.8 meV.

Figure 15: AA stacking configuration band structure of 3.89∘ $t$MoTe2. (a) and
(b) show the band structures without SOC and with SOC bands for rigid
structure, respectively. (c) and (d) show the band structures without SOC and
with SOC bands for relaxed structure, respectively Figure 16: Conduction bands
and valence bands with their density of states (DOS) in unit of states per Mo
atom of 3.89∘ for relaxed AA stacking and AB stacking configuration with SOC.
a and b are the conduction bands of AA and AB stacking configuration with DOS.
c and d are the valence bands of AA and AB stacking configuration. The green
line represents the Fermi level.

The density of states (DOS) calculations depicted in Fig. 16 (c) reveal two
distinct peaks within the valence bands around -75meV and -105meV, which
correspond to the $\Gamma$ valley band embedded in the backdrop of dispersive
bands. Additionally, several minor peaks around 0meV, 20meV, are observed,
which originate from the valence bands of the $\pm\text{K}$ valleys.

In Fig. 17, we show the band structure given by three different relaxation
methods: a, only using MLFF, b, two-step MLFF+DFT relaxation and c, relaxation
directly from rigid structure using DFT. MLFF can an accuracy considerably
close to the DFT calculations, especially the top three pairs of
$\pm\text{K}$-valley valence bands and the highest $\Gamma$-valley bands as
shown in Fig. 17(a). However, we find the forces of MLFF relaxed structure is
not small enough (with a mean absolute force of $7.8\times
10^{-2}eV/\mbox{\text{\AA}}$), meaning that the MLFF result can be relaxed
further using the more accurate DFT relaxation method. After the further DFT
relaxation, the mean absolute force reached $2.1\times
10^{-3}eV/\mbox{\text{\AA}}$, and the band structure changes quantitatively
(see in Fig. 17 (a) and (b)). If the relaxation is performed from the rigid
structure using DFT only, we will get almost the same band structure as the
two-step relaxation band.

Figure 17: Comparision of different relaxation band structures. a, band
structure of MLff relaxation. b, band structure of MLFF+DFT relaxation. c,
band structure of direct DFT relaxation from rigid structure. Figure 18:
Evolution of band structures for various angles range from 9.43∘ to 3.48∘ AB
stacking configuration $t$MoTe2. Top(bottom) row is for conduction(valence)
bands. The lattice was fully relaxed with DFT-D2 vdW functional. SOC are
considered in the calculation. Figure 19: Irreducible representations at high-
symmetry points $\Gamma$ and $K$ of valence bands of 3.89∘ for relaxed AB
stacking configuration with SOC. The eigenvalues of the $C_{3}$ are denoted as
$\omega=e^{i\pi/3}$, $\omega^{*}=e^{-i\pi/3}$ and $\bar{1}=e^{-i\pi}$

The squared wavefunctions $|\psi_{\Gamma_{M}}|^{2}$ of the top four valence
bands at the $\Gamma_{M}$ point in the K valley is illustrated in Fig. 20(a).
The wavefunctions of the first and second topmost valence bands in the K
valley are predominantly localized in the XM and MX regions, collectively
manifesting a hexagonal lattice pattern in real space. In contrast, the
wavefunction of the third topmost valence band primarily occupies the AA
region, delineating a triangular lattice. Meanwhile, the fourth band’s
wavefunction is concentrated at the XM, MX, and AA regions, constituting
another hexagonal lattice configuration.

####  A.3.5 Comparision with Ref. [20, 21]

Ref. [20] and Ref. [21] have studied the relaxation and band structures of
$t$MoTe2. Ref. [20] employed the SCAN density functional with dDsC dispersion
correction to perform crystal structure relaxation. In our assessment of 19
different functionals, as shown in Tab. 3, we find that the SCAN functional
and dDsC dispersion correction (IVDW=4 in Tab. 3) yielded a larger c-axis
lattice parameter and smaller a(b)-axis lattice parameters, whereas the DFT-D2
functional that we use (IVDW=10 in Tab. 3) provide the lattice parameters
closest to the experimental results. After crystal structure relaxation, Ref.
[20] showed a larger interlayer distance compared to our relaxation.

In contrast to both our approach and Ref. [20], Ref. [21] used SIESTA with
DFT-D2 functional to perform the DFT calculations. After the relaxation, Ref.
[21] gets a smaller inter-layer distance (ILD) of which the minimum is about
6.9 Å. This Mo-Mo ILD is smaller than our relaxation result shown in Fig. 3,
even smaller than ILD of bulk MoTe2 crystal structure with 7.0 Å.

Our relaxation results shows that, in the AA region, where the stacking
configuration is close to that of AA untwisted bilayer structure, and thus the
maximum interlayer distance should be close to but slightly smaller than (due
to the drag of surrounding environment which is not aligned in AA form) the
interlayer distance of AA stacking untwisted bilayer (7.7Å), which is not
satisfied by Ref. [20]. Furthermore, the MX region has the stacking
configuration akin to that of AB untwisted bilayer structure. As a result, the
smallest interlayer distance in $t$MoTe2 structures should be close to but
slightly larger than (due to the surrounding environment’s influence) 7Å,
which is violated by Ref. [21].

From Fig. 15, it is clear that relaxation will push the $\Gamma$-valley bands
down as well as lift the bands from $\pm\text{K}$ valleys up. The relaxation
in Ref. [21] gets the smaller ILD, resulting that their $\Gamma$-valley bands
are not pushed down significantly and close to the $\pm\text{K}$-valley bands.
The relaxed $\Gamma$-valley bands in Ref. [21] are only 30 meV below the
valence band maximum, while our $\Gamma$-valley bands are 80 meV below the
VBM. Meanwhile, the relaxation in Ref. [20] gives the maximal ILD of 7.8
Åwhich is larger that us, due to the overestimated c-axis lattice parameter by
SCAN functional. Ref. [20] calculated the band structure of 4.4∘ moiré
structure (1014 atoms per unit cell). Due to Ref. [20] only calculated 4 k
point along the high-symmetry line, the band width is estimated from their
fitted continuum model. From the model, the top valence band from
$\pm\text{K}$ valleys has band width of about 30 meV at 4.41∘, while our 4.41∘
structure (1014 atoms per unit cell) $\pm\text{K}$-valley band has a smaller
bandwidth. For 3.89∘ structure in Ref. [21], the top pair of the valence bands
have band width of about 9 meV, which is smaller than our result of 12.8 meV.

####  A.3.6 AB stacking: Evolution of band structure from 9.43∘ to 3.89∘

The symmetry group of the AB stacking can be given by replacing $C_{2y}$ in
the group for $AA$ stacking by $C_{2x}$. At large twist angle, such as 13.2∘,
the electronic structures of the two stacking configurations appear similar,
as illustrated in Fig. 12. However, they exhibit notable differences at
smaller twist angles (most notable for angles smaller than $5^{\circ}$, as
seen in Fig. 18 and Fig. 13.

As listed in Tab. 5 , the bandwidth of the top two pairs of bands narrows with
decreasing twist angle, specifically to 16 meV for 3.89∘ and 9.8 meV for
3.48∘. As the twist angle decreases, a pair of ultra flat valence bands
gradually moves up, and become isolated for $\theta<4.41^{\circ}$ as shown in
Fig. 18. When the twist angle decreases to 3.89∘, the two set of ultra flat
valence bands are separated with about 56.0 meV, bringing two distinct peak in
the DOS as shown in Fig. 16.

Similar with that of AA, the top two pairs of valence bands of AB
configuration consist of $d_{x^{2}-y^{2}}$ and $d_{xy}$ orbitals of Mo atoms
as shown in Fig. 14, indicating these bands come from the $\pm\text{K}$
valleys. Besides, the flat valence bands of AB come from $\Gamma$ valley
because they consist of Mo $d_{z^{2}}$ orbitals.

Distinct from from AA stacking, there are isolated moiré bands on the
conduction band side for the AB stacking. As the twist angle decreases to
4.41∘ (see Fig. 18(d)), there are 12 bands at the bottom of the conduction
bands that are isolated from the higher-energy bands. The degeneracy of these
bands suggests they likely originate from the pockets along the $\Gamma$-K
line of the conduction bands, as depicted in Fig. 11(b) and (c). Compared with
AA stacking, the valence bands from the $\Gamma$ valley in AB stacking are
extremely flat (see Fig. 18(k) and (l)), resulting in pronounced peaks in the
density of states (DOS) (see Fig. 16). To faciliate further comprehensive
analysis in Appendix. B and Appendix. C, we have also calculated the
irreducible representations of the top 4 valence bands for AB stacking, which
can be found in Fig. 19. In Fig. 20(b), we plot the charge densities of first
four pairs of valence bands’ wavefunctions at $\Gamma_{M}$. The top
$\pm\text{K}$ valley bands wavefunctions are localized at MX region, similar
with the top $\Gamma$ valley bands. The third pair valence bands’
wavefunctions are localized at MM region, forming triangular lattice, while
the forth pair valence bands’ wavefunctions form triangular rings with minimum
at MX region.

Figure 20: Charge densities of the valence-band wavefunctions at $\gamma$ point in relaxed 3.89∘ AA-twist and AB-twist MoTe2, The isosurface value for these plots is $2\times 10^{-5}e/\text{\AA}^{3}$ Twist angle | 9.43∘ | 7.34∘ | 5.09∘ | 4.41∘ | 3.89∘ | 3.48∘
---|---|---|---|---|---|---
Bandwidth (meV) | 135.4 | 89.8 | 36.1 | 22.5 | 16.0 | 9.8
Table 5: Band width of the top two pairs of valence band of 3.89 AB stacking
$t$MoTe2.

##  Appendix B Single-Particle Continuum Model and Fitting: AA-Stacking

The symmetry group of AA-stacked $t$MoTe2 (AA-$t$MoTe2) is generated by
$C_{3}$, $C_{2y}$, and $\mathcal{T}$, in addition to the moiré lattice
translations $T_{\bm{R}_{M}}$ [43], where $\bm{R}_{M}$ labels the moiré
lattice vectors. We label the moiré lattice basis vectors as
$\bm{a}_{M,1}=a_{M}(\frac{\sqrt{3}}{2},-\frac{1}{2})^{T}$ and
$\bm{a}_{M,2}=C_{3}\bm{a}_{M,1}$, and the moiré reciprocal lattice basis
vectors as $\bm{b}_{M,1}=\frac{4\pi}{\sqrt{3}a_{M}}(1,0)^{T}$ and
$\bm{b}_{M,2}=\frac{4\pi}{\sqrt{3}a_{M}}(\frac{1}{2},\frac{\sqrt{3}}{2})^{T}$,
where

$a_{M}=\frac{a_{0}}{2\sin\left(\frac{\theta}{2}\right)}\ ,$ (9)

and $a_{0}=3.52\text{\AA}$ is the lattice constant of the monolayer MoTe2.

Up to now, FCI states were only found for hole doping experimentally [33, 34,
35, 36]. From Fig. 13, it is immediately apparent that isolated, nearly flat
bands (understood to be an important precursor to the FCI phase [3, 2, 1])
appear only in the valence bands, which are accessible through hole doping. In
contrast the conduction bands do not have well-separated bands. As such, we
focus on building a model for the valence bands in this section.

We will disucss the continuum models for both $\pm\text{K}$ valleys and the
$\Gamma$ valley.

###  B.1 AA-Stacking: $\pm\text{K}$ Valleys

####  B.1.1 Microscopic Basis, Symmetries, and Inter-layer coupling

We now derive the moiré states that make up the continuum model basis for the
AA-stacked twisted heterostructure. This derivation follows [82, 83] by
expanding the tight-binding states around the monolayer K point. (K can be
straightforwardly obtained via time-reversal symmetry.) We consider a two-
layer system where the $l=t,b$ layer is twisted via the linear transformation
$M_{\ell}=1-\mathrm{i}\frac{(-)^{l}\theta}{2}\sigma_{y}$ to leading order
$\theta\ll 1$, where $(-)^{t}=-(-)^{b}=1$. We will derive the symmetry
representations as well as the form of the interlayer coupling within the two-
center approximation.

From our first principles calculations in Fig. 14, we see the valence band
maximum around K is spanned by $d_{x^{2}-y^{2}}+id_{xy}$ orbitals on the Mo
atoms, which we will refer to as $d$ for brevity in this part. The $d_{z^{2}}$
Mo orbitals and Te orbitals do not contribute significantly to the density of
states near the active bands. We write $M_{\ell}\mathbf{R}$ as the positions
of the Mo atoms on the $\ell$th layer, where $\mathbf{R}$ is an untwisted
lattice vector. The states carried by these orbitals are

$\displaystyle\left|{\mathbf{p},\ell}\right>=\frac{1}{\sqrt{N}}\sum_{\mathbf{R}}e^{iM_{\ell}\mathbf{R}\cdot\mathbf{p}}\left|{M_{\ell}\mathbf{R}^{\prime},\ell}\right>$
(10)

for $\ell=0,1$ corresponding to top/bottom, and
$\mathbf{p}=\mathbf{K}+\delta\mathbf{p}$ is a momentum near the $\mathbf{K}$
point. The intra-layer $C_{3}$ symmetry acts as

$\displaystyle C_{3}\left|{\mathbf{K}+\delta\mathbf{p},\ell}\right>$ (11)
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{\mathbf{R}}e^{iM_{\ell}\mathbf{R}\cdot
C_{3}(\mathbf{K}+\delta\mathbf{p})}\left|{M_{\ell}\mathbf{R},\ell}\right>e^{i\lambda_{d}}$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{\mathbf{R}}e^{iM_{\ell}\mathbf{R}\cdot(\mathbf{K}-\mathbf{G}+C_{3}\delta\mathbf{p})}\left|{M_{\ell}\mathbf{R},\ell}\right>e^{i\lambda_{d}}$
$\displaystyle=\left|{\mathbf{K}+C_{3}\delta\mathbf{p},\ell}\right>e^{i\lambda_{d}}$

where $e^{i\lambda_{d}}=e^{i\frac{2\pi}{3}}$ is the (spin-less) $C_{3}$
eigenvalue of $d$ orbital, and $\mathbf{G}=C_{3}\mathbf{K}-\mathbf{K}$ is a
$t$MoTe2reciprocal lattice vector obeying
$\mathbf{G}\cdot\mathbf{R}=0\,\mathrm{mod}\,2\pi$. We also used that the
rotation matrix $C_{3}$ commutes with $M_{\ell}$, since both are rotations.
Note that $e^{i\lambda_{d}}$ is simply an overall phase of the rotation
representation. Secondly, the $C_{2y}\mathcal{T}$ representation is

$\displaystyle
C_{2y}\mathcal{T}\left|{\mathbf{K}+\delta\mathbf{p},\ell}\right>$ (12)
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{\mathbf{R}}e^{-iM_{\ell}\mathbf{R}\cdot(\mathbf{K}+\delta\mathbf{p})}\left|{C_{2y}M_{\ell}\mathbf{R},-\ell}\right>$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{\mathbf{R}}e^{-iM_{\ell}\mathbf{R}\cdot(\mathbf{K}+\delta\mathbf{p})}\left|{M_{-\ell}C_{2y}\mathbf{R},-\ell}\right>$
$\displaystyle=\frac{1}{\sqrt{N}}\sum_{\mathbf{R}}e^{-iM_{\ell}C_{2y}\mathbf{R}\cdot(\mathbf{K}+\delta\mathbf{p})}\left|{M_{-\ell}\mathbf{R},-\ell}\right>$
$\displaystyle=\sum_{\ell^{\prime}}\left|{\mathbf{K}+C_{2y}(-\delta\mathbf{p}),\ell^{\prime}}\right>[\sigma_{x}]_{\ell\ell^{\prime}}$

using the fact that $\mathcal{T}$ is anti-unitary and $C_{2y}$ flips layer,
the reflection property $C_{2y}M_{\ell}C_{2y}^{-1}=M_{-\ell}$, and the trivial
transformation of $d$ orbital under $C_{2y}\mathcal{T}$ [31]. In this
discussion, we have neglected the spin degree of freedom, which can be easily
re-introduced because of the spin-valley locking in MoTe2. Lastly, we note
that the emergent intra-valley inversion (which will be discussed in Appendix.
B.1.2) has no representation on the microscopic basis since it is not a true
symmetry of the model, much like the emergent particle-hole symmetry in
twisted bilayer graphene [83, 84].

We now derive the form of the interlayer Hamiltonian. Formally, we compute the
overlap

$\displaystyle H^{inter}_{l,-l}(\mathbf{p},\mathbf{p}^{\prime})$
$\displaystyle=\frac{1}{N}\sum_{\mathbf{R},\mathbf{R}^{\prime}}e^{-iM_{l}\mathbf{R}\cdot\mathbf{p}+iM_{-l}\mathbf{R}^{\prime}\cdot\mathbf{p}^{\prime}}$
(13)
$\displaystyle\qquad\left<{M_{l}\mathbf{R},l}\right|H\left|{M_{-l}\mathbf{R}^{\prime},-l}\right>\
.$

It is convenient to shift the sum so that the bottom layer is unrotated and
the top layer is rotated by $-\theta$:

$\displaystyle\begin{split}&H^{inter}_{l,-l}(\mathbf{p},\mathbf{p}^{\prime})\\\
&=\frac{1}{N}\sum_{\mathbf{R},\mathbf{R}^{\prime}}e^{-iM\mathbf{R}\cdot\mathbf{p}+i\mathbf{R}^{\prime}\cdot\mathbf{p}^{\prime}}\left<{M\mathbf{R},l}\right|H\left|{\mathbf{R}^{\prime},-l}\right>\\\
\end{split}$ (14)

where $M=1-i\theta\sigma_{2}$. To proceed, we assume that the matrix element
of the Hamiltonian is only dependent on the distance between orbitals (the
“two-center" approximation) leading to

$\displaystyle\begin{split}&\left<{M\mathbf{R},l}\right|H\left|{\mathbf{R}^{\prime},-l}\right>\\\
&=\frac{1}{N\Omega}\sum_{\mathbf{q}\in
BZ}\sum_{\mathbf{G}}t_{\mathbf{q}+\mathbf{G}}e^{i(\mathbf{q}+\mathbf{G})\cdot(M\mathbf{R}-\mathbf{R}^{\prime})}\
.\end{split}$ (15)

Plugging this expression into the interlayer Hamiltonian and using
$M\mathbf{r}\cdot\mathbf{k}=(M\mathbf{r})^{T}\mathbf{k}=\mathbf{r}\cdot
M^{T}\mathbf{k}$ gives gives

$\displaystyle H^{inter}_{l,-l}(\mathbf{p},\mathbf{p}^{\prime})$
$\displaystyle=\sum_{\mathbf{G}_{1},\mathbf{G}_{2}}\frac{t_{\mathbf{p}+\mathbf{G}_{1}}}{\Omega}\delta_{\mathbf{p}+\mathbf{G}_{1},\mathbf{p}^{\prime}+M\mathbf{G}_{2}}\
.$ (16)

Using $\mathbf{p}=\mathbf{K}+\delta\mathbf{p}$ (and similarly for
$\mathbf{p}^{\prime}$), we see that
$H^{inter}_{l,-l}(\mathbf{p},\mathbf{p}^{\prime})$ connects momenta apart by
$M\mathbf{K}-\mathbf{K}=\mathbf{q}_{1}$. We now keep only the lowest order
$\mathbf{G}_{1}$ since $t_{\mathbf{p}}$ is rapidly decaying. Thus in real
space, this Hamiltonian can be written

$\displaystyle
H^{inter}(\mathbf{K}+\delta\mathbf{p},M\mathbf{K}+\delta\mathbf{p}^{\prime})$
$\displaystyle=\sum_{j=1}^{3}w\delta_{\delta\mathbf{p},\delta\mathbf{p}^{\prime}+\mathbf{q}_{i}}$
(17)

where $w=t_{\bm{K}}$. We will find agreement with this term and the symmetry-
based approach in the following section.

The intra-layer Hamiltonian, in the two-center approximation, is given by
expanding the monolayer dispersion
$h_{mono}(\mathbf{K}+\delta\mathbf{p})=-\frac{\hbar^{2}\delta\mathbf{p}^{2}}{2m^{*}}$
to leading order. Note that the two-center approximation cannot capture the
intra-layer moiré potential which arises due to relaxation within the moiré
unit cell.

Figure 21: The comparison between the valence bands from the DFT calculation
(black dots) and from the continuum model (red/orange line) for twist angle
$3.89^{\circ}$ in the fundamental zone for (a,b) AA-stacking and (c) AB-
stacking . The red line comes from the $\pm\text{K}$-valley model, while the
orange line comes fro mthe $\Gamma$-valley model. The horizontal axies labels
the 31 momenta in the $k$-mesh of the fundamental zone that we choose. Each
red line is at least doubly degenerate. The orange line in (c) has double
degeneracy around $-30$meV and 4-fold degeneracy around $-86$meV. The fitting
in (a) is done with FH terms, while SH terms are added in (b).

####  B.1.2 Continuum Model

In this part, we discuss the continuum model for the low-energy states around
the $\pm\text{K}$ valleys based on symmetries. Ref. [43] proposed a model with
first harmonics, and we will introduce more terms into it. In monolayer MoTe2,
the strong spin-orbit coupling locks the spin degree of freedom to the valley
degree of freedom. Explicitly, the highest valence band around the K valley is
made up of $d_{x^{2}-y^{2}}+\mathrm{i}d_{xy}$ Mo orbitals with spin $\uparrow$
[85]. The highest and second highest valence bands are seperated by the large
energy 200meV [43], and thus we only consider the highest electron valence
band around K and $-\text{K}$ (in both layers) to construct the moiré model.
The basis of the continuum model is labelled by $c^{\dagger}_{\eta,l,\bm{r}}$
with $\eta=\pm$ labels the $\pm\text{K}$ valleys (or equivalently spins),
$l=t,b$ labels the layer, and $\bm{r}$ labels the position. The wavefunction
of $c^{\dagger}_{\eta,l,\bm{r}}$ is a continuous approximation to the
microscopic basis in Eq. (10). Using Eqs. LABEL:eq:c2ytmicro and
LABEL:eq:c3micro and after accounting for the electron spin, the _spinful_
symmetry representations furnished by $c^{\dagger}_{\eta,l,\bm{r}}$ read

$\displaystyle\begin{split}&C_{3}c^{\dagger}_{\eta,l,\bm{r}}C_{3}^{-1}=c^{\dagger}_{\eta,l,C_{3}\bm{r}}e^{\mathrm{i}\eta\frac{\pi}{3}}\\\
&(C_{2y}\mathcal{T})c^{\dagger}_{\eta,l,\bm{r}}(C_{2y}\mathcal{T})^{-1}=c^{\dagger}_{\eta,\bar{l},C_{2y}\bm{r}}\\\
&\mathcal{T}c^{\dagger}_{\eta,l,\bm{r}}\mathcal{T}^{-1}=c^{\dagger}_{-\eta,l,\bm{r}}(-\eta)\\\
&T_{\bm{R}_{M}}c^{\dagger}_{\eta,l,\bm{r}}T_{\bm{R}_{M}}^{-1}=c^{\dagger}_{\eta,l,\bm{r}+\bm{R}_{M}}e^{-\mathrm{i}\eta\bm{R}_{M}\cdot\text{K}_{l}}\
,\end{split}$ (18)

where $\bar{l}=b,t$ for $l=t,b$ respectively, and $\text{K}_{l}$ is the
rotated K point of the $l$th layer. For AA-$t$MoTe2, we rotate the top layer
by $-\theta/2$ and the bottom layer by $\theta/2$, and thus

$\displaystyle\begin{split}&\text{K}_{b}=\frac{4\pi}{3a_{0}}\left(\cos\left(\frac{\theta}{2}\right),\sin\left(\frac{\theta}{2}\right)\right)^{T}\\\
&\text{K}_{t}=\frac{4\pi}{3a_{0}}\left(\cos\left(\frac{\theta}{2}\right),-\sin\left(\frac{\theta}{2}\right)\right)^{T}\
.\end{split}$ (19)

For the convenience of the discussions in the rest of this part, we define

$\displaystyle\begin{split}&\bm{q}_{1}=\text{K}_{b}-\text{K}_{t}=\frac{4\pi}{3a_{0}}2\sin\left(\frac{\theta}{2}\right)\left(\begin{matrix}0\\\
1\end{matrix}\right)\\\ &\bm{q}_{2}=C_{3}\bm{q}_{1}\ ,\
\bm{q}_{3}=C_{3}^{2}\bm{q}_{1}\ .\end{split}$ (20)

The general form of the single-particle Hamiltonian in the $\eta$ valley reads

$H_{\eta,0}^{AA}=\int
d^{2}r\left(c^{\dagger}_{\eta,b,\bm{r}},c^{\dagger}_{\eta,t,\bm{r}}\right)\left(\begin{matrix}h_{\eta,b}(\bm{r})&t_{\eta}(\bm{r})\\\
t_{\eta}^{*}(\bm{r})&h_{\eta,t}(\bm{r})\end{matrix}\right)\left(\begin{matrix}c_{\eta,b,\bm{r}}\\\
c_{\eta,t,\bm{r}}\end{matrix}\right)\ ,$ (21)

where $h_{\eta,l}(\bm{r})$ is hermitian. The symmetry requirements on the
terms in the Hamiltonian read

$\displaystyle\begin{split}&C_{3}:\
h_{\eta,l}(C_{3}\bm{r})=h_{\eta,l}(\bm{r})\ ,\
t_{\eta}(C_{3}\bm{r})=t_{\eta}(\bm{r})\\\ &C_{2y}\mathcal{T}:\
h_{\eta,\bar{l}}(C_{2y}\bm{r})=h_{\eta,l}^{*}(\bm{r})\ ,\
t_{\eta}(C_{2y}\bm{r})=t_{\eta}(\bm{r})\\\ &\mathcal{T}:\
h_{-\eta,l}(\bm{r})=h_{\eta,l}^{*}(\bm{r})\ ,\
t_{-\eta}(\bm{r})=t_{\eta}^{*}(\bm{r})\\\
&T_{\bm{R}_{M}}:h_{\eta,l}(\bm{r}+\bm{R}_{M})=h_{\eta,l}(\bm{r})\ ,\\\
&\qquad\ \
t_{\eta}(\bm{r}+\bm{R}_{M})=t_{\eta}(\bm{r})e^{-\mathrm{i}\eta\bm{q}_{1}\cdot\bm{R}_{M}}\end{split}$
(22)

Owing to $C_{3}$ symmetry, the dispersion around $\pm\text{K}$ in the
monolayer MoTe2 is quadratic. Thus, the spatial derivatives are kept to the
second order in the intra-layer term:

$h_{\eta,l}(\bm{r})=\frac{\hbar^{2}\nabla^{2}}{2m^{*}}+V_{\eta,l}(\bm{r})\ ,$
(23)

where $m^{*}$ is real and positive as we are considering the valence electron
bands, the $C_{2y}\mathcal{T}$ and $\mathcal{T}$ symmetries require $m^{*}$ to
be the same for all values of $\eta,l$, and hermiticity requires
$V_{\eta,l}(\bm{r})$ to be real. The kinetic term follows from the two-center
approximation (see the discussion after Eq. 17). Owing to the constraints
imposed by moiré lattice translations on the potential terms (see Eq.
(LABEL:eq:symmetry_requirements_Kvalley)), $V_{\eta,l}(\bm{r})$ and
$e^{\mathrm{i}\eta\bm{q}_{1}\cdot\bm{r}}t_{\eta}(\bm{r})$ can be expanded in
series of the moiré lattice vector $\bm{G}_{M}$:

$\displaystyle\begin{split}&V_{\eta,l}(\bm{r})=\sum_{\bm{G}_{M}}e^{-\eta\mathrm{i}\bm{G}_{M}\cdot\bm{r}}V_{\eta,l,\bm{G}_{M}}\\\
&t_{\eta}(\bm{r})=\sum_{\bm{G}_{M}}e^{-\eta\mathrm{i}(\bm{q}_{1}+\bm{G}_{M})\cdot\bm{r}}t_{\eta,\bm{q}_{1}+\bm{G}_{M}}\
,\end{split}$ (24)

where the TR symmetry requires
$V_{\eta,l,\bm{G}_{M}}=V_{-\eta,l,\bm{G}_{M}}^{*}$ and
$t_{\eta,l,\bm{q}_{1}+\bm{G}_{M}}=t_{-\eta,l,\bm{q}_{1}+\bm{G}_{M}}^{*}$.

In Ref. [43], only the first harmonics of $V_{\eta,l}(\bm{r})$ and
$t_{\eta}(\bm{r})$ is kept, _i.e._ , only including
$\bm{G}_{M}\in\\{\pm\bm{g}_{i}|i=1,2,3\\}$ with
$\bm{g}_{i}=C_{3}^{i-1}\bm{b}_{M,1}$ for the intralayer terms, and only
include $\bm{q}_{1}+\bm{G}_{M}\in\\{-\bm{q}_{i}|i=1,2,3\\}$ for the interlayer
terms. Combined with $C_{3}$, $C_{2y}\mathcal{T}$ and TR symmetries,
$V_{+,l}(\bm{r})$ and $t_{+}(\bm{r})$ take the forms of

$\displaystyle\begin{split}&V_{\eta,l}^{FH}(\bm{r})=Ve^{-(-)^{l}\mathrm{i}\psi}\sum_{i=1,2,3}e^{\mathrm{i}\bm{g}_{i}\cdot\bm{r}}+Ve^{(-)^{l}\mathrm{i}\psi}\sum_{i=1,2,3}e^{-\mathrm{i}\bm{g}_{i}\cdot\bm{r}}\\\
&t_{\eta}^{FH}(\bm{r})=w\sum_{i=1,2,3}e^{-\eta\mathrm{i}\bm{q}_{i}\cdot\bm{r}}\
,\end{split}$ (25)

where “$FH$” labels the first harmonics, $(-)^{t}=1$, $(-)^{b}=-1$, $V$ and
$\psi$ are real, and $w$ is chosen to be real non-positive by tuning the
relative phase between the two layers. The lowest-harmonics model has
effective inversion symmetry

$\mathcal{I}c^{\dagger}_{\eta,l,\bm{r}}\mathcal{I}^{-1}=c^{\dagger}_{\eta,\bar{l},-\bm{r}}\
.$ (26)

Adding displacement field would induce an energy difference between the two
layers, _i.e._ ,

$H_{\eta,\varepsilon}=\int
d^{2}r\sum_{l}c^{\dagger}_{\eta,l,\bm{r}}c_{\eta,l,\bm{r}}(-)^{l}\frac{\varepsilon}{2}\
,$ (27)

which breaks the effective inversion symmetry, as well as the $C_{2y}$
symmetry.

As discussed in the previous work (_e.g._ , Ref. [21]) and in Appendix. B.3,
the lowest-harmonics model can only well describe the highest two valence
bands per valley and the gap between 2nd and 3rd valence bands at $\Gamma$. In
our work, we want to capture the 3rd band in each valley by including one more
harmonics. In general, including higher harmonics may break the effective
inversion symmetry. However, according to the DFT band structure in Fig. 13,
the 3rd pairs of valence bands are still approximately degenerate for angles
below 5.09 degrees; therefore, we will keep the effective inversion symmetry
when adding extra harmonics. The second harmonics (SH) have
$|\bm{G}_{M}|=|\bm{b}_{M,1}+\bm{b}_{M,2}|$ for the intra-layer potential and
have $|\bm{q}_{1}+\bm{G}_{M}|$ = $\bm{q}_{1}+\bm{b}_{M,1}$ for the inter-layer
potential; combined with $C_{3}$, $C_{2y}\mathcal{T}$, TR and the effective
inversion symmetries, the form of the SH for the intra-layer potential read

$\displaystyle\begin{split}&V_{\eta,l}^{SH}(\bm{r})=2V_{2}\sum_{i=1}^{3}\cos(\bm{g}_{2i}\cdot\bm{r})\\\
&t_{\eta}^{SH}(\bm{r})=w_{2}\sum_{i=1,2,3}e^{-\eta\mathrm{i}\bm{q}_{2i}\cdot\bm{r}}\
,\end{split}$ (28)

where $\bm{g}_{21}=\bm{b}_{M,1}+\bm{b}_{M,2}$,
$\bm{g}_{2i}=C_{3}^{i-1}\bm{g}_{21}$,
$\bm{q}_{21}=\bm{b}_{M,1}+\bm{q}_{{}_{1}}$,
$\bm{q}_{2i}=C_{3}^{i-1}\bm{q}_{21}$, and we have enforced the effective
inversion symmetry, which makes $V_{2}$ and $w_{2}$ real. With these extra
terms, the final form of $V_{\eta}(\bm{r})$ and $t_{\eta}(\bm{r})$ reads

$\displaystyle\begin{split}&V_{\eta,l}(\bm{r})=Ve^{-(-)^{l}\mathrm{i}\psi}\sum_{i=1,2,3}e^{\mathrm{i}\bm{g}_{i}\cdot\bm{r}}+Ve^{(-)^{l}\mathrm{i}\psi}\sum_{i=1,2,3}e^{-\mathrm{i}\bm{g}_{i}\cdot\bm{r}}\\\
&\qquad+2V_{2}\sum_{i=1}^{3}\cos(\bm{g}_{2i}\cdot\bm{r})\\\
&t_{\eta}(\bm{r})=w\sum_{i=1,2,3}e^{-\eta\mathrm{i}\bm{q}_{i}\cdot\bm{r}}+w_{2}\sum_{i=1,2,3}e^{-\eta\mathrm{i}\bm{q}_{2i}\cdot\bm{r}}\
,\end{split}$ (29)

The moiré translations allow us to express the Hamiltonian in the momentum
space. The Fourier transformation of the basis reads

$c^{\dagger}_{\eta,l,\bm{r}}=\frac{1}{\sqrt{\mathcal{V}}}\sum_{\bm{k}}\sum_{\bm{Q}\in\mathcal{Q}_{l}^{\eta}}e^{-\mathrm{i}(\bm{k}-\bm{Q})\cdot\bm{r}}c^{\dagger}_{\eta,l,\bm{k}-\bm{Q}}\
,$ (30)

where

$\mathcal{Q}_{l}^{\eta}=\\{\bm{G}_{M}+\eta(-)^{l}\bm{q}_{1}\\}\ .$ (31)

As a result, $H_{\eta,0}$ reads

$H_{\eta,0}^{AA}=\sum_{\bm{k}}\sum_{\bm{Q},\bm{Q}^{\prime}\in\mathcal{Q}}c^{\dagger}_{\eta,\bm{k},\bm{Q}}h_{\eta,\bm{Q}\bm{Q}^{\prime}}^{AA}(\bm{k})c_{\eta,\bm{k},\bm{Q}^{\prime}}\
,$ (32)

where we have included the displacement field,
$\mathcal{Q}=\\{\bm{G}_{M}+\bm{q}_{1}\\}\cup\\{\bm{G}_{M}-\bm{q}_{1}\\}$,

$c^{\dagger}_{\eta,\bm{k},\bm{Q}}=c^{\dagger}_{\eta,l_{\bm{Q}}^{\eta},\bm{k}-\bm{Q}}\
,$ (33)

and $l_{\bm{Q}}^{\eta}=l$ for $\bm{Q}\in\mathcal{Q}_{l}^{\eta}$.The general
form of $h_{\eta,\bm{Q}\bm{Q}^{\prime}}^{AA}(\bm{k})$ reads

$\displaystyle\begin{split}&h_{\eta,\bm{Q}\bm{Q}^{\prime}}^{AA}(\bm{k})=\delta_{\bm{Q}\bm{Q}^{\prime}}\left(\frac{-\hbar^{2}\left(\bm{k}-\bm{Q}\right)^{2}}{2m^{*}}+\frac{\varepsilon}{2}\eta(-)^{\bm{Q}}\right)\\\
&\quad+\sum_{\bm{G}_{M}}\delta_{\bm{Q}_{1},\bm{Q}_{2}+\eta\bm{G}_{M}}V_{\eta,l_{\bm{Q}_{1}}^{\eta},\bm{G}_{M}}\\\
&\quad+\sum_{\bm{G}_{M}}\delta_{\bm{Q}_{1},\bm{Q}_{2}+\eta(\bm{q}_{1}+\bm{G}_{M})}t_{\eta,l_{\bm{Q}_{1}}^{\eta},\bm{q}_{1}+\bm{G}_{M}}\\\
&\quad+\sum_{\bm{G}_{M}}\delta_{\bm{Q}_{1},\bm{Q}_{2}-\eta(\bm{q}_{1}+\bm{G}_{M})}t_{\eta,l_{\bm{Q}_{1}}^{\eta},\bm{q}_{1}+\bm{G}_{M}}^{*}\
.\end{split}$ (34)

With the FH and SH terms in Eq. (LABEL:eq:Kvalley_V_t_FH_SH), the matrix
Hamiltonian $h_{\eta,\bm{Q}\bm{Q}^{\prime}}^{AA}(\bm{k})$ has the form

$\displaystyle\begin{split}&h_{\eta,\bm{Q}\bm{Q}^{\prime}}^{AA}(\bm{k})=\delta_{\bm{Q}\bm{Q}^{\prime}}\left(\frac{-\hbar^{2}\left(\bm{k}-\bm{Q}\right)^{2}}{2m^{*}}+\frac{\varepsilon}{2}\eta(-)^{\bm{Q}}\right)\\\
&\quad+V\sum_{i=1}^{3}\left[e^{-\eta(-)^{\bm{Q}}\mathrm{i}\psi}\delta_{\bm{Q}+\bm{g}_{i},\bm{Q}^{\prime}}+e^{\eta(-)^{\bm{Q}}\mathrm{i}\psi}\delta_{\bm{Q}-\bm{g}_{i},\bm{Q}^{\prime}}\right]\\\
&\quad+V_{2}\sum_{i=1}^{3}\left[\delta_{\bm{Q}+\bm{g}_{2i},\bm{Q}^{\prime}}+\delta_{\bm{Q},\bm{Q}^{\prime}+\bm{g}_{2i}}\right]\\\
&\quad+w\sum_{i=1}^{3}\left[\delta_{\bm{Q}+\bm{q}_{i},\bm{Q}^{\prime}}+\delta_{\bm{Q},\bm{Q}^{\prime}+\bm{q}_{i}}\right]\\\
&\quad+w_{2}\sum_{i=1}^{3}\left[\delta_{\bm{Q}+\bm{q}_{2i},\bm{Q}^{\prime}}+\delta_{\bm{Q},\bm{Q}^{\prime}+\bm{q}_{2i}}\right]\
,\end{split}$ (35)

where $(-)^{\bm{Q}}=\pm 1$ for $\bm{Q}\in\\{\bm{G}_{M}\pm\bm{q}_{1}\\}$, and
we used $\eta(-)^{\bm{Q}}=(-)^{l^{\eta}_{\bm{Q}}}$

###  B.2 AA-Stacking: $\Gamma$ Valley

In this part, we discuss the continuum model for the $\Gamma$ valley. In the
monolayer MoTe2, the bands around $\Gamma$ mainly come from the $d_{z^{2}}$
orbital, and time-reversal symmetry at $\Gamma$ leads to a Kramers degeneracy
between spin up and spin down states. Thus, the basis of the moiré model at
$\Gamma$ valley is labelled as $\psi^{\dagger}_{\bm{r},l,s}$, where $l$ is the
layer index, and $s=\uparrow/\downarrow$ labels the spin. Then, the continuum
model in general has the following form

$H_{\Gamma}^{AA}=\int
d^{2}r\left(\begin{matrix}\psi^{\dagger}_{\bm{r},b}&\psi^{\dagger}_{\bm{r},t}\end{matrix}\right)\left(\begin{matrix}h_{\Gamma,b}(\bm{r})&t_{\Gamma}(\bm{r})\\\
t_{\Gamma}^{\dagger}(\bm{r})&h_{\Gamma,t}(\bm{r})\end{matrix}\right)\left(\begin{matrix}\psi^{\dagger}_{\bm{r},b}\\\
\psi^{\dagger}_{\bm{r},t}\end{matrix}\right)\ ,$ (36)

where

$\psi^{\dagger}_{\bm{r},l}=(\psi^{\dagger}_{\bm{r},l,\uparrow},\psi^{\dagger}_{\bm{r},l,\downarrow})\
,$ (37)

$t$ and $b$ correspond to the top and bottom layers, respectively, and
$h_{\Gamma,l}(\bm{r})$ and $t_{\Gamma}(\bm{r})$ are $2\times 2$ moiré-periodic
matrix functions.

The spinful representations of the $d_{z^{2}}$ orbital symmetries of the moiré
model are

$\displaystyle\begin{split}&C_{3}\psi^{\dagger}_{\bm{r},l}C_{3}^{-1}=\psi^{\dagger}_{C_{3}\bm{r},l}e^{-\mathrm{i}s_{z}\frac{\pi}{3}}\\\
&C_{2y}\mathcal{T}\psi^{\dagger}_{\bm{r},l}(C_{2y}\mathcal{T})^{-1}=\psi^{\dagger}_{C_{2y}\bm{r},\bar{l}}\\\
&\mathcal{T}\psi^{\dagger}_{\bm{r},l}\mathcal{T}^{-1}=\psi^{\dagger}_{\bm{r},l}\mathrm{i}s_{y}\\\
&T_{\bm{R}_{M}}\psi^{\dagger}_{\bm{r},l}T_{\bm{R}_{M}}^{-1}=\psi^{\dagger}_{\bm{r}+\bm{R}_{M},l}\
,\end{split}$ (38)

where $s_{x},s_{y},s_{z}$ are the Pauli matrices for the spin index, and
$\bar{l}=t,b$ if $l=b,t$. Based on symmetries, we know

$\displaystyle\begin{split}&e^{-\mathrm{i}s_{z}\frac{\pi}{3}}h_{\Gamma,l}(\bm{r})e^{\mathrm{i}s_{z}\frac{\pi}{3}}=h_{\Gamma,l}(C_{3}\bm{r})\\\
&e^{-\mathrm{i}s_{z}\frac{\pi}{3}}t_{\Gamma}(\bm{r})e^{\mathrm{i}s_{z}\frac{\pi}{3}}=t_{\Gamma}(C_{3}\bm{r})\\\
&h_{\Gamma,l}^{*}(\bm{r})=h_{\Gamma,\bar{l}}(C_{2y}\bm{r})\ ,\
t_{\Gamma}^{*}(\bm{r})=t_{\Gamma}^{\dagger}(C_{2y}\bm{r})\\\
&\sigma_{y}h_{\Gamma,l}^{*}(\bm{r})\sigma_{y}=h_{\Gamma,l}(\bm{r})\ ,\
\sigma_{y}t_{\Gamma}^{*}(\bm{r})\sigma_{y}=t_{\Gamma}(\bm{r})\\\
&h_{\Gamma,l}(\bm{r})=h_{\Gamma,l}(\bm{r}+\bm{R}_{M})\ ,\
t_{\Gamma}(\bm{r})=t_{\Gamma}(\bm{r}+\bm{R}_{M})\ .\end{split}$ (39)

The intra-layer kinetic energy term must take the form

$\displaystyle\begin{split}h_{\Gamma,l}(\bm{r})=\frac{\hbar^{2}\nabla^{2}}{2m_{\Gamma}^{*}}s_{0}+V_{\Gamma,l}(\bm{r})\end{split}$
(40)

since a linear Rashba-like kinetic term is forbidden by $C_{3}$,
$\mathcal{T}$, and $M_{z}$ symmetry in the monolayer. The effective mass
$m_{\Gamma}^{*}$ is the same in both layers because of the $C_{2y}$ symmetry.
We note that the value of $m_{\Gamma}^{*}$ in Eq. (40) might not be equal to
that for monolayer MoTe2 due to the inter-layer coupling renormalizing the
bands, due to the fact that the bands near $\Gamma$ in the monolayer is quite
flat as shown in Fig. 11 (a). The potentials are periodic and can be expanded
as

$\displaystyle\begin{split}&V_{\Gamma,l}(\bm{r})=\sum_{\bm{G}}V_{\Gamma,l,\bm{G}}e^{\mathrm{i}\bm{G}\cdot\bm{r}}\\\
&t_{\Gamma}(\bm{r})=\sum_{\bm{G}}t_{\Gamma,\bm{G}}e^{\mathrm{i}\bm{G}\cdot\bm{r}}\end{split}$
(41)

whose components are restricted by symmetry to obey

$\displaystyle\begin{split}&e^{-\mathrm{i}s_{z}\frac{\pi}{3}}V_{\Gamma,l,\bm{G}}e^{\mathrm{i}s_{z}\frac{\pi}{3}}=V_{\Gamma,l,C_{3}\bm{G}}\
,\
e^{-\mathrm{i}s_{z}\frac{\pi}{3}}t_{\Gamma,\bm{G}}e^{\mathrm{i}s_{z}\frac{\pi}{3}}=t_{\Gamma,C_{3}\bm{G}}\\\
&V_{\Gamma,l,\bm{G}}^{*}=V_{\Gamma,\bar{l},-C_{2y}\bm{G}}\ ,\
t_{\Gamma,\bm{G}}^{*}=t^{\dagger}_{\Gamma,-C_{2y}\bm{G}}\\\
&s_{y}V_{\Gamma,l,\bm{G}}^{*}s_{y}=V_{\Gamma,l,-\bm{G}}\ ,\
s_{y}t_{\Gamma,\bm{G}}^{*}s_{y}=t_{\Gamma,-\bm{G}}\\\
&V_{\Gamma,l,\bm{G}}^{\dagger}=V_{\Gamma,l,-\bm{G}}\ .\end{split}$ (42)

To capture the dominant contribution, we only include the first harmonics in
$V_{\Gamma,l}(\bm{r})$, and the zeroth harmonic in the $t_{\Gamma}(\bm{r})$,
_i.e._ ,

$\displaystyle\begin{split}&V_{\Gamma,l}(\bm{r})=\sum_{\bm{G}\in\\{\pm\bm{g}_{i}|i=1,2,3\\}}V_{\Gamma,l,\bm{G}}e^{\mathrm{i}\bm{G}\cdot\bm{r}}\\\
&t_{\Gamma}(\bm{r})=t_{\Gamma,0}\ .\end{split}$ (43)

We neglect the zeroth harmonic in $V_{\Gamma,l}(\bm{r})$ since it is just a
total shift of energy (although it must be included to compare the energies of
the in-active $\Gamma$ valley and active $\pm\text{K}$ valleys). According to
Eq. (LABEL:eq:sym_rep_AA_Gamma_Valley), the first harmonic terms in
$V_{\Gamma,l}(\bm{r})$ can be written

$\displaystyle\begin{split}&V_{\Gamma,b,\pm\bm{g}_{i}}=V_{\Gamma}e^{\pm\mathrm{i}\psi_{\Gamma}}s_{0}\\\
&V_{\Gamma,t,\pm\bm{g}_{i}}=V_{\Gamma}e^{\mp\mathrm{i}\psi_{\Gamma}}s_{0}\
,\end{split}$ (44)

and similarly, $t_{\Gamma,0}$ has the form

$t_{\Gamma,0}=w_{\Gamma}e^{\mathrm{i}s_{z}\phi_{\Gamma}}\ .$ (45)

Here $V_{\Gamma}$, $\psi_{\Gamma}$, $w_{\Gamma}$ and $\phi_{\Gamma}$ are real.
We set $\phi_{\Gamma}=0$ as a gauge choice of the relative phase between the
layers. As a result, Eq. (LABEL:eq:Vt_AA_Gamma_Valley) is further simplified
to

$\displaystyle\begin{split}&V_{\Gamma,l}(\bm{r})=\sum_{i=1}^{3}V_{\Gamma}e^{\mathrm{i}(-)^{l}\psi_{\Gamma}}e^{\mathrm{i}\bm{g}_{i}\cdot\bm{r}}+\sum_{i=1}^{3}V_{\Gamma}e^{-\mathrm{i}(-)^{l}\psi_{\Gamma}}e^{-\mathrm{i}\bm{g}_{i}\cdot\bm{r}}\\\
&t_{\Gamma}(\bm{r})=w_{\Gamma}s_{0}\ .\end{split}$ (46)

With the simplification in Eq. (40) and Eq. (LABEL:eq:Vt_AA_Gamma_Valley_sim),
the $\Gamma$-valley continuum model in Eq. (36) has effective inversion
symmetry, _i.e._ ,

$\mathcal{I}\psi^{\dagger}_{\bm{r},l}\mathcal{I}^{-1}=\psi^{\dagger}_{-\bm{r},\bar{l}}\
.$ (47)

Similar to the $\pm\text{K}$-valley case, adding displacement breaks the
effective inversion symmetry, as well as the $C_{2y}$ symmetry, where the
displacement field term reads

$H_{\Gamma,\varepsilon}=\int
d^{2}r\sum_{l}\psi^{\dagger}_{\bm{r},l}\psi_{\bm{r},l}(-)^{l}\frac{\varepsilon}{2}\
.$ (48)

Including the first harmonics in $t_{\Gamma}(\bm{r})$ and the higher harmonics
in $V_{\Gamma,l}(\bm{r})$ may also break the effective inversion symmetry.

Figure 22: The $C_{3}$ eigenvalues of the bands from the continuum model for
twist angle $3.89^{\circ}$ for (a,b) AA-stacking and (c) AB-stacking . The
plot in (a) is done with FH terms, while SH terms are added in (b). In (a,b),
we only include $C_{3}$ eigenvalues for the K-valley bands (red) at
$\Gamma_{M}$ and $\text{K}_{M}$, since $K_{M}^{\prime}$ point is related to
the $K_{M}$ point by $C_{2y}\mathcal{T}$ and the $-\text{K}$ valley can be
obtained from the TR symmetry. In (c), we only include $C_{3}$ eigenvalues for
the K-valley bands (red) and the $\Gamma$-valley spin-up bands at $\Gamma_{M}$
and $\text{K}_{M}$, since $K_{M}^{\prime}$ point is related to the $K_{M}$
point by the effective TR symmetry, and those for the $-\text{K}$ valley and
the $\Gamma$-valley spin down can be obtained from the TR symmetry. Like in
Fig. 2, $\omega=e^{\mathrm{i}\pi/3}$ and $\bar{1}=-1$.

###  B.3 AA-Stacking: Fitting to the DFT Data

As shown in Fig. 13, the $\Gamma$-valley valence bands are below the 6th
highest valence bands for $\theta\leq 5.09^{\circ}$ and zero displacement
field $\varepsilon=0$. Therefore, we only use the $\pm\text{K}$-valley model
(Eq. (LABEL:eq:Kvalley_model_AA)) in the fitting.

We fit the DFT bands at $3.89^{\circ}$ in two ways. First, we set
$V_{2}=w_{2}=0$, which corresponds to the FH model. In this case, we manage to
fit the top 4 valence bands (2 in each valley) with the corresponding FH
parameters in Tab. 1, as shown in Fig. 4 (a) along the high-symmetry lines and
in Fig. 21 (a) for the $\bm{k}$ points in the fundamental zone. Then, we allow
nonzero $V_{2}$ and $w_{2}$, which means we add the effective-inversion-
symmetric SH terms. We are now able to fit the top 6 valence bands (3 in each
valley) with the corresponding FH+SH parameters in Tab. 1, as shown in Fig. 4
(b) along the high-symmetry lines and in Fig. 21 (b) for the $\bm{k}$ points
in the fundamental zone.

As shown in Tab. 1, the value of the effective mass $m^{*}$ is similar to the
monolayer and untwisted AA-stacking bilayer structures masses, as shown in
Fig. 11. We note that all bands from the $-\text{K}$-valley continuum model
are the same as those from the K-valley continuum model, due to the
combination of the effective inversion and the TR symmetry. By comparing Fig.
22 (a,b) to Fig. 2 (c), we can see that the $C_{3}$ eigenvalues for the top 6
valence bands from the model are the same as those from the DFT calculation in
both FH and FH+SH cases at $3.89^{\circ}$. At $3.89^{\circ}$, The Chern
numbers of the top 3 bands in K valley are (1,-1,0) for the highest, the
second highest and the third highest valence bands in both FH and FH+SH cases,
respectively, which are consistent with the $C_{3}$ eigenvalues. The Chern
numbers and the $C_{3}$ eigenvalues of the top two bands per valley are the
same as those in Ref. [21] at $3.89^{\circ}$. The quantum geometry on the
moiré Brillouin zone is shown in Fig. 23-24, with Chern number $C$ and
integrated Fubini-Study metric $G/2\pi$ (see Ref. [31]). We observe that the
Berry curvature and Fubini-Study metric are more uniform in the first band
than in the second (remote) band, with the second remote band showing a peak
around $\Gamma$ where the gap to the third band is smallest. Note also that
the effect of second harmonic terms is weak on the first ($C=-1$) valence
band, whereas the second valence band becomes significantly more strongly
peaked at the $\Gamma$ point and flatter elsewhere on the BZ. This can be
understood from the decrease in the gap at $\Gamma$ between the second and
third valence bands due to the improved accuracy of the higher-harmonic model.

Figure 23: Quantum geometry of the AA-stacked continuum FH model
($\theta=3.89^{\circ}$) in (LABEL:eq:Kvalley_model_AA) (without higher
harmonics). $(a,b)$ show the Berry curvature and Fubini-Study metric of the
highest valence band, respectively, and $(c,d)$ show the corresponding plots
for the second highest valence band. The Chern numbers are $C=-1,+$ for $(a)$
and $(c)$ respectively, and integrated Fubini-Study metrics are marked.
Figure 24: Quantum geometry of the AA-stacked continuum model
($\theta=3.89^{\circ}$) including second-shell harmonics in
(LABEL:eq:Kvalley_model_AA). $(a,b)$ show the Berry curvature and Fubini-Study
metric of the highest valence band, respectively, and $(c,d)$ show the
corresponding plots for the second highest valence band. The Chern numbers are
$C=-1,+$ for $(a)$ and $(c)$ respectively.

Since the FH+SH model captures more bands than the FH model, we use the FH+SH
model to discuss the phase diagram under the twist angle $\theta$ and the
displacement field $\varepsilon$. Fig. 5 in the Main Text shows the single-
particle phase diagram for the top three valence bands in K valley. The
highest valence band – most relevant to the filling factors where FCIs have
been observed. The top valence K-valley band displays a Chern number $C=1$
throughout the phase diagram at $\varepsilon=0$, and shows a phase transition
into a trivial insulator as $\varepsilon$ is turned on. The phase diagram of
the remote bands is richer, showing all Chern numbers $C=-2,\dots,2$
throughout the $\theta,\varepsilon$ plane. In particular, at $\varepsilon=0$,
the second and third top bands have a gap closing around $\theta=4.2^{\circ}$,
which is close to the gap closing around $4.41^{\circ}$ shown in Fig. 2.

##  Appendix C Single-Particle Continuum Model and Fitting: AB-Stacking
(2H-Stacking)

Replacing $C_{2y}$ by $C_{2x}$ in the generators of the symmetry group of
AA-$t$MoTe2 gives the generators of the symmetry group of AB-stacking $t$MoTe2
(AB-$t$MoTe2), which can be thought of as twisting the top layer of
AA-$t$MoTe2 by another $180^{\circ}$. The change of symmetry has a dramatic
effect on the model. Since $C_{2x}$ is local to the monolayer K point (unlike
$C_{2y}$), it preserves the valley quantum in the moiré model.

Different from the AA-$t$MoTe2, the DFT results show that the low-energy
valence bands of AB-$t$MoTe2 come from both the $\pm\text{K}$ valleys and the
$\Gamma$ valleys. Therefore, in the following, we will review the continuum
model for the $\pm\text{K}$ valleys proposed in Ref. [43], and will also
discuss the continuum model for the $\Gamma$ valley.

###  C.1 AB-Stacking: $\pm\text{K}$ Valleys

In this part, we review the continuum model of AB-$t$MoTe2 for the
$\pm\text{K}$ valleys proposed in Ref. [43]. The basis of the continuum model
is still labelled by $c^{\dagger}_{\eta,l,\bm{r}}$ with $\eta=\pm$ labeling
the $\pm\text{K}$ valleys (or equivalently spins), $l=t,b$ labels the layer,
and $\bm{r}$ labels the position. Since AB-$t$MoTe2 is given by rotating the
top layer of AA-$t$MoTe2 by an extra $180^{\circ}$, the two layers in one
valley now have opposite spin. Then, the symmetry reps furnished by
$c^{\dagger}_{\eta,l,\bm{r}}$ now read

$\displaystyle\begin{split}&C_{3}c^{\dagger}_{\eta,l,\bm{r}}C_{3}^{-1}=c^{\dagger}_{\eta,l,C_{3}\bm{r}}e^{-\mathrm{i}\eta(-)^{l}\frac{\pi}{3}}\\\
&C_{2x}c^{\dagger}_{\eta,l,\bm{r}}C_{2x}^{-1}=c^{\dagger}_{\eta,\bar{l},C_{2x}\bm{r}}(-\mathrm{i})\\\
&\mathcal{T}c^{\dagger}_{\eta,l,\bm{r}}\mathcal{T}^{-1}=c^{\dagger}_{-\eta,l,\bm{r}}(-\eta)\\\
&T_{\bm{R}_{M}}c^{\dagger}_{\eta,l,\bm{r}}T_{\bm{R}_{M}}^{-1}=c^{\dagger}_{\eta,l,\bm{r}+\bm{R}_{M}}e^{-\mathrm{i}\eta\bm{R}_{M}\cdot\text{K}_{l}}\
,\end{split}$ (49)

where $\bar{l}=b,t$ for $l=t,b$ respectively, $\text{K}_{l}$ is the rotated K
point of the $l$th layer, and $(-)^{t}=-(-)^{b}=1$.

The general form of the single-particle Hamiltonian in the $\eta$ valley reads

$\displaystyle\begin{split}H_{\eta,0}^{AB}&=\sum_{l}\int
d^{2}rc^{\dagger}_{\eta,l,\bm{r}}\left[\frac{\hbar^{2}\nabla^{2}}{2m^{*}}+V_{\eta,l}(\bm{r})\right]c_{\eta,l,\bm{r}}\\\
&\ +\left[\int
d^{2}rc^{\dagger}_{\eta,b,\bm{r}}t_{\eta}(\bm{r})c_{\eta,t,\bm{r}}+h.c.\right]\
,\end{split}$ (50)

where $V_{\eta,l}(\bm{r})$ is real, and we have used the fact that the
dispersion around $\pm\text{K}$ in the monolayer MoTe2 is quadratic and the
fact that $m^{*}$ is the same for two valleys and two layers owing to $C_{2x}$
and $\mathcal{T}$ symmetries. The symmetry requirements on
$V_{\eta,l}(\bm{r})$ and $t_{\eta}(\bm{r})$ read

$\displaystyle\begin{split}&C_{3}:\
V_{\eta,l}(C_{3}\bm{r})=V_{\eta,l}(\bm{r})\ ,\
t_{\eta}(C_{3}\bm{r})=t_{\eta}(\bm{r})e^{-\mathrm{i}\eta\frac{2\pi}{3}}\\\
&C_{2x}:\ V_{\eta,l}(C_{2x}\bm{r})=V_{\eta,\bar{l}}(\bm{r})\ ,\
t_{\eta}^{\dagger}(C_{2x}\bm{r})=t_{\eta}(\bm{r})\\\ &\mathcal{T}:\
V_{-\eta,l}(\bm{r})=V_{\eta,l}^{*}(\bm{r})\ ,\
t_{-\eta}(\bm{r})=t_{\eta}^{*}(\bm{r})\\\
&T_{\bm{R}_{M}}:V_{\eta,l}(\bm{r}+\bm{R}_{M})=V_{\eta,l}(\bm{r})\\\ &\qquad\ \
t_{\eta}(\bm{r}+\bm{R}_{M})=t_{\eta}(\bm{r})e^{\mathrm{i}\eta\bm{q}_{1}\cdot\bm{R}_{M}}\end{split}$
(51)

By keeping the first harmonics, the symmetry requirements lead to

$\displaystyle
V_{\eta,l}(\bm{r})=Ve^{-\mathrm{i}\psi}\sum_{i=1,2,3}e^{\mathrm{i}\bm{g}_{i}\cdot\bm{r}}+Ve^{\mathrm{i}\psi}\sum_{i=1,2,3}e^{-\mathrm{i}\bm{g}_{i}\cdot\bm{r}}$
(52) $\displaystyle
t_{\eta}(\bm{r})=w\sum_{i=1,2,3}e^{\eta\mathrm{i}(i-1)\frac{2\pi}{3}}e^{-\eta\mathrm{i}\bm{q}_{i}\cdot\bm{r}}\
.$

As can be checked from the expressions above, this lowest-harmonics model has
effective anti-unitary symmetry

$\mathcal{K}c^{\dagger}_{\eta,l,\bm{r}}\mathcal{K}^{-1}=c^{\dagger}_{\eta,\bar{l},\bm{r}}$
(53)

which flips the layers and acts like time-reversal. This symmetry can be
anticipated since the opposite layers have opposite spins in the twisted AB
stacking. Adding displacement field would induce an energy difference between
the two layers, _i.e._ ,

$H_{\eta,\varepsilon}=\int
d^{2}r\sum_{l}c^{\dagger}_{\eta,l,\bm{r}}c_{\eta,l,\bm{r}}(-)^{l}\frac{\varepsilon}{2}\
,$ (54)

which breaks the effective TR symmetry, as well as the $C_{2x}$ symmetry.
Including higher harmonics for the intralayer and interlatyer terms may break
the effective TR symmetry.

Similar to the discussion in Appendix. B.1, we can express the Hamiltonian in
the momentum space based on the Fourier transformation in Eq. (30), leading to

$H_{\eta,0}^{AB}=\sum_{\bm{k}}\sum_{\bm{Q},\bm{Q}^{\prime}\in\mathcal{Q}}c^{\dagger}_{\eta,\bm{k},\bm{Q}}h_{\eta,\bm{Q}\bm{Q}^{\prime}}^{AB}(\bm{k})c_{\eta,\bm{k},\bm{Q}^{\prime}}\
,$ (55)

where $c^{\dagger}_{\eta,\bm{k},\bm{Q}}$ is defined in Eq. (33), and
$\mathcal{Q}$ is defined right below Eq. (32). Specifically,
$h_{\eta,\bm{Q}\bm{Q}^{\prime}}^{AB}(\bm{k})$ reads

$\displaystyle\begin{split}h_{\eta,\bm{Q}\bm{Q}^{\prime}}^{AB}(\bm{k})&=\delta_{\bm{Q}\bm{Q}^{\prime}}\left(\frac{-\hbar^{2}\left(\bm{k}-\bm{Q}\right)^{2}}{2m^{*}}+\frac{\varepsilon}{2}\eta(-)^{\bm{Q}}\right)\\\
&\quad+V\sum_{i=1}^{3}\sum_{\alpha=\pm}\left[e^{-\alpha\mathrm{i}\psi}\delta_{\bm{Q}+\alpha\bm{g}_{i},\bm{Q}^{\prime}}\right]\\\
&\quad+w\sum_{i=1}^{3}\sum_{\alpha=\pm}\left[e^{\mathrm{i}\alpha\frac{2\pi(i-1)}{3}}\delta_{\bm{Q},\bm{Q}^{\prime}+\alpha\bm{q}_{i}}\right]\
,\end{split}$ (56)

where we have included the displacement field, and $w$ is chosen to be real by
picking the relative phase between the two layers, and $(-)^{\bm{Q}}=\pm 1$
for $\bm{Q}\in\\{\bm{G}_{M}\pm\bm{q}_{1}\\}$. We can see that without the
displacement field, $h_{\eta,\bm{Q}\bm{Q}^{\prime}}^{AB}(\bm{k})$ is the same
for two valleys due to the effective time-reversal symmetry.

###  C.2 AB-Stacking: $\Gamma$ Valley

In this part, we discuss the continuum model for the $\Gamma$ valley. Thus,
the basis of the moire model at $\Gamma$ valley for AB-$t$MoTe2 is the same as
that for AA-$t$MoTe2 which is labelled as $\psi^{\dagger}_{\bm{r},l,s}$ with
$l$ being the layer index and $s=\uparrow/\downarrow$ labeling the spin. The
important difference is that, in the $\Gamma$ valley, the momentum space
lattice is triangular instead of a honeycomb at $\pm\text{K}$ points. The reps
of the symmetries of the moiré model read

$\displaystyle\begin{split}&C_{3}\psi^{\dagger}_{\bm{r},l}C_{3}^{-1}=\psi^{\dagger}_{C_{3}\bm{r},l}e^{-\mathrm{i}s_{z}\frac{\pi}{3}}\\\
&C_{2x}\psi^{\dagger}_{\bm{r},l}(C_{2x})^{-1}=\psi^{\dagger}_{C_{2x}\bm{r},\bar{l}}(-\mathrm{i}s_{x})\\\
&\mathcal{T}\psi^{\dagger}_{\bm{r},l}\mathcal{T}^{-1}=\psi^{\dagger}_{\bm{r},l}\mathrm{i}s_{y}\\\
&T_{\bm{R}_{M}}\psi^{\dagger}_{\bm{r},l}T_{\bm{R}_{M}}^{-1}=\psi^{\dagger}_{\bm{r}+\bm{R}_{M},l}\
,\end{split}$ (57)

where $s_{x},s_{y},s_{z}$ are the Pauli matrices for the spin index, and
$\bar{l}=t,b$ if $l=b,t$.

After choosing the kinetic term in the continuum model as the intralayer spin-
independent $\nabla^{2}$ term similar to the AA-stacking case, the form of the
continuum model reads

$H_{\Gamma}=\int
d^{2}r\left(\begin{matrix}\psi^{\dagger}_{\bm{r},b}&\psi^{\dagger}_{\bm{r},t}\end{matrix}\right)h_{\Gamma}^{AB}(\bm{r})\left(\begin{matrix}\psi^{\dagger}_{\bm{r},b}\\\
\psi^{\dagger}_{\bm{r},t}\end{matrix}\right)\ ,$ (58)

where

$h_{\Gamma}^{AB}(\bm{r})=\left(\begin{matrix}\frac{\hbar^{2}\nabla^{2}}{2m_{\Gamma}^{*}}s_{0}+V_{\Gamma,b}(\bm{r})&t_{\Gamma}(\bm{r})\\\
t_{\Gamma}^{\dagger}(\bm{r})&\frac{\hbar^{2}\nabla^{2}}{2m_{\Gamma}^{*}}s_{0}+V_{\Gamma,t}(\bm{r})\end{matrix}\right)+E_{\Gamma}\
,$ (59)
$\psi^{\dagger}_{\bm{r},l}=(\psi^{\dagger}_{\bm{r},l,\uparrow},\psi^{\dagger}_{\bm{r},l,\downarrow})\
,$ (60)

$t$ and $b$ correspond to the top and bottom layers, respectively,
$E_{\Gamma}$ accounts for the energy difference between the $\Gamma$-valley
and $\pm\text{K}$-valley bands, and $V_{\Gamma,l}(\bm{r})$ and
$t_{\Gamma}(\bm{r})$ are $2\times 2$ matrix functions. Similar to the
discussion for AA-$t$MoTe2, we only include the terms up to the first
harmonics in $V_{\Gamma,l}(\bm{r})$, and the zeroth harmonics in the
$t_{\Gamma}(\bm{r})$. Then, combined with the symmetry reps in Eq.
(LABEL:eq:sym_rep_Gamma_AB), we arrive at

$\displaystyle\begin{split}&V_{\Gamma,l}(\bm{r})=\left[\sum_{i=1}^{3}V_{\Gamma}e^{\mathrm{i}\psi_{\Gamma}}e^{\mathrm{i}\bm{g}_{i}\cdot\bm{r}}+\sum_{i=1}^{3}V_{\Gamma}e^{-\mathrm{i}\psi_{\Gamma}}e^{-\mathrm{i}\bm{g}_{i}\cdot\bm{r}}\right]s_{0}\\\
&t_{\Gamma}(\bm{r})=w_{\Gamma}s_{0}\ .\end{split}$ (61)

With the simplification in Eq. (59) and Eq. (LABEL:eq:Vt_AB_Gamma_Valley_sim),
the $\Gamma$-valley continuum model in Eq. (58) has an effective TR symmetry

$\mathcal{K}\psi^{\dagger}_{\bm{r},l}\mathcal{K}^{-1}=\psi^{\dagger}_{\bm{r},\bar{l}}$
(62)

arising because of the spin-layer locking in the AB stacked structure. Similar
to the $\pm\text{K}$-valley case, adding displacement breaks the effective
inversion symmetry, as well as the $C_{2x}$ symmetry, where the displacement
field term reads

$H_{\Gamma,\varepsilon}=\int
d^{2}r\sum_{l}\psi^{\dagger}_{\bm{r},l}\psi_{\bm{r},l}(-)^{l}\frac{\varepsilon}{2}\
.$ (63)

Including the first harmonics in $t_{\Gamma}(\bm{r})$ and the higher harmonics
in $V_{\Gamma,l}(\bm{r})$ can also break the effective TR symmetry.

To express the Hamiltonian in the momentum space, we use the Fourier
transformation of the basis, which reads

$\psi^{\dagger}_{\bm{r},l,s}=\frac{1}{\sqrt{\mathcal{V}}}\sum_{\bm{k}}\sum_{\bm{G}_{M}}e^{-\mathrm{i}(\bm{k}-\bm{G}_{M})\cdot\bm{r}}\psi^{\dagger}_{\bm{k}-\bm{G}_{M},l,s}\
.$ (64)

As shown in Eq. (LABEL:eq:Vt_AB_Gamma_Valley_sim), the Hamiltonian has the
same form in the two spin subspaces. As a result, $H_{\eta,0}$ reads

$H_{\Gamma,0}^{AB}=\sum_{\bm{k},\bm{G}_{M},\bm{G}_{M}^{\prime}}\sum_{l,l^{\prime},s}\psi^{\dagger}_{\bm{k}-\bm{G}_{M},l,s}h_{ll^{\prime},\bm{G}_{M}\bm{G}_{M}^{\prime}}^{AB}(\bm{k})\psi_{\bm{k}-\bm{G}_{M}^{\prime},l^{\prime},s}\
,$ (65)

where

$\displaystyle\begin{split}&h_{ll,\bm{G}_{M}\bm{G}_{M}^{\prime}}^{AB}(\bm{k})=\left[-\frac{\hbar^{2}(\bm{k}-\bm{G})^{2}}{2m_{\Gamma}^{*}}+(-)^{l}\frac{\varepsilon}{2}+E_{\Gamma}\right]\delta_{\bm{G}_{M},\bm{G}_{M}^{\prime}}\\\
&\qquad+V_{\Gamma}\sum_{i=1}^{3}\left[e^{\mathrm{i}\psi_{\Gamma}}\delta_{\bm{G}_{M},\bm{G}_{M}^{\prime}-\bm{g}_{i}}+e^{-\mathrm{i}\psi_{\Gamma}}\delta_{\bm{G}_{M},\bm{G}_{M}^{\prime}+\bm{g}_{i}}\right]\\\
&h_{l\bar{l},\bm{G}_{M}\bm{G}_{M}^{\prime}}^{AB}(\bm{k})=w_{\Gamma}\delta_{\bm{G}_{M},\bm{G}_{M}^{\prime}}\
,\end{split}$ (66)

and we have included the displacement field.

###  C.3 AB-Stacking: Fitting to the DFT Data

As shown in Fig. 18, both the $\pm\text{K}$-valley and the $\Gamma$-valley
valence bands are close to the charge neutrality. Therefore, we only use both
the $\pm\text{K}$-valley model (Eq. (56)) and the $\Gamma$-valley model (Eq.
(LABEL:eq:Gammavalley_model_AB)) in the fitting. We mainly fit to the top 4
valence bands (2 in each valley) for the $\pm\text{K}$-valley model (Eq. (56))
and top 6 valence bands in the $\Gamma$ valley; similar to the case of AA
stacking, we fit to the DFT bands at $\theta=3.89^{\circ}$. The good match
between those DFT bands and those from the models is shown in Fig. 4 (b) and
Fig. 21 (b), where the model parameter values obtained from the fitting are
shown in Tab. 2. In the following, we provide more details for the fitting for
the $\pm\text{K}$-valley model (Eq. (56)) and the $\Gamma$-valley model,
separately.

####  C.3.1 $\pm\text{K}$-valley

We choose the value of $m^{*}$ to be the effective mass ($\sim 0.62m_{e}$) of
the top valence band for the untwisted AB-stacking bilayer structure around K,
as shown in Fig. 11. The value of effective mass around K is similar for the
monolayer MoTe2 and the untwisted AB-stacking bilayer MoTe2.

Owing to the effective TR symmetry Eq. (53) and the true microscopic TR
symmetry, the $-\text{K}$ valley has the same bands as the $+\text{K}$ valley,
leading to at least double degeneracy of each band from the
$\pm\text{K}$-valley continuum model. Along the $\Gamma_{M}-M_{M}$ line, the
bands from the $\pm\text{K}$-valley continuum model have 4-fold degeneracy
because we set the interlayer coupling to be zero ($w=0$), which is justified
by the nearly 4-fold degeneracy for the top 4 DFT valence bands along
$\Gamma_{M}-M_{M}$, as shown in Fig. 4. The small $w=0$ can be understood as
the follows. The two layers in one valley now have opposite spin; owing to the
spin U(1) symmetry for the low-energy states near $\pm\text{K}$ valleys in
monolayer MoTe2, we expect the spin $U(1)$ symmetry is approximately preserved
in $t$MoTe2, which means the interlayer coupling is very small for the AB-
stacking. The zero interlayer coupling makes the eigenstates have well-defined
valley and layer/spin. As a result, the two states with the same spin are
degenerate at $\text{K}_{M}$, since the combination of the effective TR
symmetry and the TR symmetry leaves the spin invariant. This is consistent
with the DFT result in Fig. 19. In Fig. 22 (b), we show the $C_{3}$
eigenvalues of the bands from the $\pm\text{K}$-valley model, which is
consistent with Fig. 19 and which shows that each of the top four valence
bands is in the $A_{1}@1a$ atomic limit, representing an atomic $s$ orbital at
the origin of the center of moiré Wigner-Seitz unit cell.

####  C.3.2 $\Gamma$-valley

In our fitting, we choose the values of $m_{\Gamma}^{*}$ to be the effective
mass of the top valence band for the untwisted AB-stacking bilayer structure
around $\Gamma$ in Fig. 11. (There is some anisotropy, so in practice we round
the effective mass to $10m_{e}$. The resulting bands show little dependence on
the precise value of the mass.) Although the value of effective mass around K
is similar for the monolayer MoTe2 and the untwisted AB-stacking bilayer
MoTe2, it differs a lot around $\Gamma$ as shown in Fig. 11. The change of the
$\Gamma$ effective mass from the monolayer MoTe2 to the untwisted AB-stacking
bilayer MoTe2 comes from the interlayer coupling between the monolayer
conduction and valence bands. This effect can only be taken into account in
the $\Gamma$-valley continuum model by changing the value of $m^{*}$, since
the the monolayer conduction bands are not explicitly included in the
$\Gamma$-valley continuum model. Therefore, we do not take the values of
$m_{\Gamma}^{*}$ from the monolayer; instead, we round it to the $\Gamma$
effective mass in the untwisted AB-stacking bilayer MoTe2. Furthermore, we
determine the interlayer coupling $w_{\Gamma}$ for the $\Gamma$-valley model
based on the interlayer coupling at $\Gamma$ for the untwisted AB-stacking
bilayer structure, as shown in Fig. 11.

The $\Gamma$-valley continuum model in Eq. (58) preserves spin
$\mathrm{U}(1)$. Owing to the effective TR symmetry Eq. (53) and the TR
symmetry, the spin-down bands are the same as the spin-up bands, leading to at
least double degeneracy of each band from the $\Gamma$-valley continuum model.
As shown in Fig. 22 (c), the $C_{3}$ eigenvalues of the bands from the
$\pm\text{K}$-valley model are consistent with the DFT calculation shown in
Fig. 19.

The $\Gamma$-valley bands are extremely flat as shown in Fig. 4. To understand
it, first we note that the interlayer coupling $w_{\Gamma}$ is very large
(shown in Tab. 2), allowing us to project the model into
$\psi^{\dagger}_{\bm{r},s}=\frac{1}{\sqrt{2}}(\psi^{\dagger}_{t,\bm{r},s}+\psi^{\dagger}_{b,\bm{r},s})$
basis. (Note that we have rotated the relative phase between
$\psi^{\dagger}_{t,\bm{r},s}$ and $\psi^{\dagger}_{b,\bm{r},s}$ to make
$w_{\Gamma}$ positive.) Then, the $\Gamma$-valley continuum model in Eq. (59)
with zero displacement field becomes

$\tilde{h}^{AB}(\bm{r})=\frac{\hbar^{2}\nabla^{2}}{2m_{\Gamma}^{*}}+w_{\Gamma}+E_{\Gamma}+V_{\Gamma}(\bm{r})\
,$ (67)

where the spin degeneracy is implicit, and

$\displaystyle\begin{split}V_{\Gamma}(\bm{r})&=\sum_{i=1}^{3}V_{\Gamma}e^{\mathrm{i}\psi_{\Gamma}}e^{\mathrm{i}\bm{g}_{i}\cdot\bm{r}}+\sum_{i=1}^{3}V_{\Gamma}e^{-\mathrm{i}\psi_{\Gamma}}e^{-\mathrm{i}\bm{g}_{i}\cdot\bm{r}}\\\
&=2V_{\Gamma}\sum_{i=1}^{3}\cos(\bm{g}_{i}\cdot\bm{r}+\psi_{\Gamma})\
.\end{split}$ (68)

For convenience of the discussion, let us consider the hole Hamiltonian

$-\tilde{h}^{AB}(\bm{r})=-\frac{\hbar^{2}\nabla^{2}}{2m_{\Gamma}^{*}}-w_{\Gamma}-E_{\Gamma}-V_{\Gamma}(\bm{r})\
;$ (69)

the highest bands of $\tilde{h}^{AB}(\bm{r})$ are lowest bands of
$-\tilde{h}^{AB}(\bm{r})$. With $\psi=0$ and $V_{\Gamma}>0$ in Tab. 2, the
minimum of the potential $-V_{\Gamma}(\bm{r})$ is at the moiré lattice points
$\bm{r}=\bm{R}_{M}$, which has value $-V_{\Gamma}(\bm{R}_{M})=-6V_{\Gamma}$.
To go from one moiré lattice point $\bm{R}_{M}$ to its nearest neighbors
$\bm{R}_{M}$ to $\bm{R}_{M}+\bm{a}_{M,1}$, the minimum energy barrier is
$8V_{\Gamma}=576$meV. The energy barrier is much larger than the kinetic
energy term estimated for $m^{*}_{\Gamma}=10m_{e}$:
$\frac{\hbar^{2}|\bm{q}_{1}|^{2}}{2m^{*}_{\Gamma}}=2.49$meV. As a result, the
lowest-energy states of $-\tilde{h}^{AB}(\bm{r})$ are bounded around
$\bm{R}_{M}$; then we can expand $-V_{\Gamma}(\bm{r})$ around $\bm{R}_{M}$ to
the second order of $\bm{r}-\bm{R}_{M}$ (the first order of
$\bm{r}-\bm{R}_{M}$ is forbidden by the $C_{3}$ symmetry):

$\displaystyle\begin{split}-V_{\Gamma}(\bm{r})&=V_{\Gamma}\sum_{i}(\bm{g}_{i}\cdot\bm{r})^{2}+O(|\bm{r}|^{4})\\\
&=\frac{3}{2}V_{\Gamma}|b_{M,1}|^{2}\left|\bm{r}\right|^{2}+O(|\bm{r}|^{4})\
.\end{split}$ (70)

Then, the effective Hamiltonian for those lowest-energy states can be
approximated by an array of decoupled harmonic oscillators centered around the
$\bm{R}_{M}$. Harmonic oscillators centered around different $\bm{R}_{M}$’s
are related by the moiré lattice translations, and thus have the same
spectrum. As a result, we have a set of flat bands with energies given by
solving

$h_{\Gamma,harmonic}^{AB}(\bm{r})=-\frac{\hbar^{2}\nabla^{2}}{2m_{\Gamma}^{*}}+w_{\Gamma}-E_{\Gamma}+\frac{3}{2}V_{\Gamma}|b_{M,1}|^{2}\left|\bm{r}\right|^{2}$
(71)

which yields

$E_{\Gamma,harmonic}^{AB}(n_{x},n_{y})=w_{\Gamma}E_{\Gamma}-\hbar\omega-(n_{x}+n_{y})\frac{\hbar\omega}{2}\
,$ (72)

where $n_{x},n_{y}=0,1,2,3,...$ and

$\hbar\omega=\hbar\sqrt{\frac{3V_{\Gamma}|b_{M,1}|^{2}}{m_{\Gamma}^{*}}}=56.76\text{meV}\
.$ (73)

This spectrum tells us (i) after including spin, the highest valence band
$(n_{x},n_{y})=(0,0)$ is doubly degenerate and the second highest valence band
$(n_{x},n_{y})=(1,0),(0,1)$ is 4-fold degenerate, (ii) the highest and second
valence bands have spinless $C_{3}$ eigenvalues $1$ and
$e^{\pm\mathrm{i}2\pi/3}$, respectively, and (iii) the energy difference
between the highest valence band and the second highest valence is
$\hbar\omega=56.76\text{meV}$. All of the features are consistent with the DFT
calculation. In Fig. 21 (b), the flat bands from the $\Gamma$ valley around
$-30$meV and $-86$meV have near 2-fold degeneracy (spinless $C_{3}$ eigenvalue
$1$) and near 4-fold degeneracy (spinless $C_{3}$ eigenvalue
$e^{\pm\mathrm{i}2\pi/3}$), and the energy difference between them is
$56.74$meV.

At last, we discuss the effect of the displacement field. As shown in Fig. 25,
the effect of the displacement field on the $\pm K$-valley bands are just to
shift the bands from different layers relative to each other, since the valley
is a good quantum number due to the zero interlayer coupling (see Tab. 2). On
the other hand, the effect of the displacement field on the low-energy
$\Gamma$-valley bands is negligible, which is consistent with the fact that
the very large interlayer (see Tab. 2) coupling make the eigenstates equally
distributed between the two layers.

Figure 25: The evolution of the bands of the AB-stacking $t$MoTe2 based on the
continuum model with parameter values in Tab. 2. The red and orange lines come
from the $\pm\text{K}$ valley and the $\Gamma$ valley, respectively.

###  C.4 AB-Stacking Electron Bands

For completeness, we also briefly discuss the conduction bands although they
are not the focus of this work. The Fermi surface of the electron bands in the
untwisted structure is more complex than the valence bands (see Fig. 11),
exhibiting an electron-like pocket along the $\Gamma\text{K}$ line which is
below the $\pm\text{K}$-valley states in both bilayer configurations. Let us
denote the center of the pocket by $s\mathbf{K}$, $s\sim.7$, along with its
$C_{3}$\- and $\mathcal{T}$-related partners. Thus in total there are 6
pockets in the untwisted BZ. Fig. 11 shows that there are only two low-energy
states (spin $\uparrow$ and spin $\downarrow$) at $s\mathbf{K}$, unlike the
four total states (spin $\uparrow$ and spin $\downarrow$ in both layers) at
$\mathbf{K}$. This can be understood from the simple model
$\frac{\mathbf{k}^{2}}{2m_{p}}\tau_{0}-w_{p}\tau_{1}$ where $\tau$ is a layer
Pauli matrix, $\mathbf{k}$ is the momentum measured from the Fermi pockets
with mass $m_{p}$, and $w_{p}$ is the effective inter-layer coupling, which
splits the two layers into bonding/anti-bonding states with energy
$\mathbf{k}^{2}/2m_{p}\mp w_{p}$.

Comparing the monolayer and bilayer band structures in Fig. 11, we estimate
$w_{p}\sim 200$meV. Since this is a large hybridization that removes the anti-
bonding state from the low-energy bands, an effective moiré model on the
bonding state is

$\displaystyle H_{\mathbf{p}_{i}}$
$\displaystyle=\frac{(-i\nabla+\mathbf{p}_{i})^{2}}{2m_{p}}+V_{p}(\mathbf{r})$
(74)

where $V_{p}(\mathbf{r})$ is a moiré-periodic potential and
$\mathbf{p}_{j}=R(\frac{2\pi j}{6})s\mathbf{K}$ is the center of the pocket.
Note that the momentum-space origin of $H_{\mathbf{p}_{i}}$, $\mathbf{p}_{i}$
mod $\mathbf{G}_{M}$, depends sensitively on the value of $s$ and $\theta$ and
is generically off the moiré high-symmetry points. The resulting shifting and
folding of the $\mathbf{p}_{i}$ pocket onto the moiré BZ explains the non-
uniform evolution of the conduction bands in Fig. 13 and 18 as a function of
twist angle. In comparison, the valence bands which originate at the K points
converge smoothly.

We show band structures at $\theta=4.41^{\circ},3.89^{\circ}$ in Fig. 26 using
$s=.68$ as an example to illustrate the strong effect of $\mathbf{p}_{i}$.
Here we chose similar parameters to the $K$-valley model in Tab. 2, namely
$m_{p}=.62m_{e}$, $V=26.5$meV (half the $K$ valley value), and
$\psi=-52^{\circ}$. We see that Fig. 26 qualitatively matches the lowest bands
in Fig. 18. Since the only symmetry preserved at $\mathbf{p}_{i}$ is $C_{2x}$,
symmetry does not strongly constrain the form of Eq. 74. Nevertheless, the
minimal model here successfully reproduces many features of the bands.

Figure 26: Conduction bands in the AB-stacked structure using $s=.68$ for
$\theta=4.41^{\circ}$ (a) and $\theta=3.89^{\circ}$ (b). Note that there are 6
low energy bands (per valley) in each case, although in (a) they are nearly
degenerate. The features of the band structures are well-matched by the
_lowest_ 12 bands in ab-initio in Fig. 18 between both angles. One should note
that the ab-initio band structures also show bands from the $K$ valley at
slightly higher energies.

## References

  * Neupert  _et al._ [2011] T. Neupert, L. Santos, C. Chamon, and C. Mudry, Phys. Rev. Lett. 106, 236804 (2011).
  * Sheng  _et al._ [2011] D. N. Sheng, Z.-C. Gu, K. Sun, and L. Sheng, Nature Communications 2, 389 (2011), arXiv:1102.2658 [cond-mat.str-el] .
  * Regnault and Bernevig [2011] N. Regnault and B. A. Bernevig, Phys. Rev. X 1, 021014 (2011).
  * Tang  _et al._ [2011] E. Tang, J.-W. Mei, and X.-G. Wen, Phys. Rev. Lett. 106, 236802 (2011).
  * Sun  _et al._ [2011] K. Sun, Z. Gu, H. Katsura, and S. Das Sarma, Phys. Rev. Lett. 106, 236803 (2011).
  * Bergholtz and Liu [2013] E. J. Bergholtz and Z. Liu, International Journal of Modern Physics B 27, 1330017 (2013), https://doi.org/10.1142/S021797921330017X .
  * Parameswaran  _et al._ [2013] S. A. Parameswaran, R. Roy, and S. L. Sondhi, Comptes Rendus Physique 14, 816 (2013), topological insulators / Isolants topologiques.
  * Abouelkomsan  _et al._ [2020] A. Abouelkomsan, Z. Liu, and E. J. Bergholtz, Phys. Rev. Lett. 124, 106803 (2020).
  * Ledwith  _et al._ [2020] P. J. Ledwith, G. Tarnopolsky, E. Khalaf, and A. Vishwanath, Phys. Rev. Res. 2, 023237 (2020).
  * Repellin and Senthil [2020] C. Repellin and T. Senthil, Phys. Rev. Res. 2, 023238 (2020).
  * Parker  _et al._ [2021] D. Parker, P. Ledwith, E. Khalaf, T. Soejima, J. Hauschild, Y. Xie, A. Pierce, M. P. Zaletel, A. Yacoby, and A. Vishwanath, “Field-tuned and zero-field fractional chern insulators in magic angle graphene,” (2021), arXiv:2112.13837 [cond-mat.str-el] .
  * Wilhelm  _et al._ [2021] P. Wilhelm, T. C. Lang, and A. M. Läuchli, Phys. Rev. B 103, 125406 (2021).