$V_{ISP}$, $D[\mathcal{I}\mathcal{T}]$ is broken by the potential difference
$V$, and both $D[\mathcal{I}\mathcal{T}]$ and $M_{1}\mathcal{T}$ are broken by
the moiré potential term proportional to $V_{1}$ in general. If $\theta=0$ so
that the hBN axis is aligned with the R$n$G axis, then the $M_{1}\mathcal{T}$
symmetry is restored. However, we focus on $\theta=0.76715^{\circ}$ throughout
this work.

Taking into account the lowest order moiré potential (parametrized by
$V_{0},V_{1}$ and $\psi_{\xi}$), as well as the effect of the out-of-plane
displacement field (parametrized by $V$ the inter-layer potential energy
difference), the R$n$G/hBN effective model acquires the following form:

$H^{\text{eff}}_{n,\xi}=h^{\text{eff}}_{n}(-i\nabla)+[V_{0}+2V_{1}\sum_{i}\cos(\mathbf{g}_{j}\cdot\mathbf{r}+\psi_{\xi})]\begin{pmatrix}1&0\\\
0&0\\\ \end{pmatrix}\ ,$ (85)

where the second term is obtained from $\mathcal{U}^{\dagger}(0)V_{\text{hBN,
bottom}}(\mathbf{r})\mathcal{U}(0)$, with $\mathcal{U}(0)$ the chiral basis at
$\mathbf{k}=0$ (see Eq. (45)), $V_{\text{hBN, bottom}}$ defined in Eq. (72)
and $\mathbf{g}_{j}\equiv
R(\frac{2\pi}{3}(j-1))(\mathbf{q}_{2}-\mathbf{q}_{3})$. This projection
preserves the locality of the moiré potential. For the case of hBN/R$n$G/hBN
structures, we have

$\begin{split}H^{\text{eff}}_{n,\xi_{t},\xi_{b}}=h^{\text{eff}}_{n}(-i\nabla)+&\begin{pmatrix}V_{b0}+2V_{b1}\sum_{i}\cos(\mathbf{g}_{i}\cdot\mathbf{r}+\psi_{\xi_{b}})&0\\\
0&V_{t0}+2V_{t1}\sum_{i}\cos(\mathbf{g}_{i}\cdot\mathbf{r}+\psi_{\xi_{t}}+\frac{2\pi
n}{3})\\\ \end{pmatrix}\ .\end{split}$ (86)

For the $(\xi_{t},\xi_{b})=(1,0)$ configuration, we impose
$\mathcal{I}\mathcal{T}$ symmetry so that $V_{b0}=V_{t0}$ and
$-\psi_{b}=\psi_{t}+\frac{2\pi n}{3}$.

Figure 18: 3D schematic of the pentalayer graphene (R5G) from the side $(a)$
and the top $(b)$. The unit cell consists of 10 carbon atoms and is connected
by a line to guide the eye. Red and blue are used to color the A and B
sublattices in each layer.

###  A.5 $C_{3}$ eigenvalues for the continuum model R$n$G/hBN

Now we take the effective model in Eq. (85) and analyze the $C_{3}$ symmetry
eigenvalues at $\tilde{\Gamma}_{M}$, $\tilde{K}_{M}$ and
$\tilde{K}^{\prime}_{M}$ points in the continuum model moiré Brillouin zone
for the graphene $K$-valley. This Brillouin zone, as well as the labeling of
high symmetry points, should be distinguished from the moiré Brillouin zone at
a commensurate twist, as depicted in Fig. 19. At different commensurate twist
angles, the graphene $\mathbf{K}_{G}$-point can be folded onto either
$\Gamma_{M}$, $K_{M}$ or $K^{\prime}_{M}$, as explained in Eq. (3), while for
simplicity of our discussion below, we will boost the BZ such that
$\mathbf{K}_{G}$ is consistently situated at the BZ center, which is denote as
$\tilde{\Gamma}_{M}$.

Figure 19: (a) Depiction of the moiré Brillouin zone (BZ) for three different
commensurate moiré configurations (see Eq. (3)). Depending on $m-n\text{ mod
}3$, with $(m,n)\in\mathbb{Z}^{2}$ labeling the commensurate configuration,
the graphene $K$-point (at $\mathbf{K}_{G}$) is either folded onto the moiré
$K_{M}$, $\Gamma_{M}$ or $K^{\prime}_{M}$ point. (b) Depiction of the
continuum model moiré BZ in which we focus on the degrees of freedom centered
around the K graphene valley. The high symmetry points are labeled as
$\tilde{\Gamma}_{M}$, $\tilde{K}_{M}$ and $\tilde{K}^{\prime}_{M}$ to signify
our convention of always boosting the graphene $\mathbf{K}_{G}$ onto
$\tilde{\Gamma}_{M}$. The three moiré reciprocal points circled in red are
considered in the tripod model analysis in App. A.5 for deducing the $C_{3}$
eigenvalue at $\tilde{K}_{M}$.

Let us just recall that the $D[C_{3}]$ operator in the $2n\times 2n$ model
(c.f. Eq. (43)) takes the form

$D[C_{3}]_{ll^{\prime}}=-\delta_{ll^{\prime}}e^{i\frac{2\pi}{3}(l-1-\lfloor\frac{n}{2}\rfloor)}e^{-i\frac{\pi}{3}\sigma_{3}},$
(87)

so by projecting onto the chiral basis we obtain
$D_{\text{eff}}[C_{3}]=\mathcal{U}^{\dagger}(0)D[C_{3}]\mathcal{U}(0)$ with
the form

$D_{\text{eff}}[C_{3}]=\begin{pmatrix}(\omega^{*})^{\frac{n}{2}-1}&0\\\
0&\omega^{\frac{n}{2}+1}\end{pmatrix}\quad\text{for even }n,\quad
D_{\text{eff}}[C_{3}]=\begin{pmatrix}(\omega^{*})^{\lfloor\frac{n}{2}\rfloor-1}&0\\\
0&\omega^{\lfloor\frac{n}{2}\rfloor-1}\end{pmatrix}\quad\text{for odd }n,$
(88)

where $\omega=e^{2i\pi/3}$. We shall focus on the lowest conduction band, and
denote the corresponding eigenvalue as $\eta_{\mathbf{k}}$. To proceed
analytically, we consider a reduced Hamiltonian that involves only the
reciprocal lattice points closest to the high symmetry point of interest.
Recall that the effective Hamiltonian in the full reciprocal space takes the
following form:

$H^{\text{eff}}_{n}(\mathbf{k})_{\mathbf{G},\mathbf{G}^{\prime}}=\delta_{\mathbf{G},\mathbf{G}^{\prime}}h^{\text{eff}}_{n}(\mathbf{k}-\mathbf{G})+\sigma_{11}(\delta_{\mathbf{G},\mathbf{G}^{\prime}}V_{0}+\delta_{\mathbf{G}-\mathbf{G}^{\prime},-\mathbf{g}_{j}}V_{1}e^{i\psi_{\xi}}+\delta_{\mathbf{G}-\mathbf{G}^{\prime},\mathbf{g}_{j}}V_{1}e^{-i\psi_{\xi}})$
(89)

where $\sigma_{11}\equiv(\sigma_{0}+\sigma_{3})/2$.

At $\tilde{\Gamma}_{M}$ (i.e. $\mathbf{k}=0$), there is one closest reciprocal
lattice point (i.e., $\mathbf{G}=0$), from which we obtain a 2 $\times$ 2
model describing the perturbation of the zero modes

$H^{\text{eff}}_{n}(\tilde{\Gamma}_{M})=(2V_{ISP}+\frac{V_{0}}{2})\mathds{1}+(\frac{V_{0}}{2}-2V)\sigma_{3}$
(90)

and we can neglect higher shells due to their large kinetic energy. Below we
consider $n=5$ for concreteness, though it is straightforward to generalize
our results to arbitrary $n$ by implementing the appropriate representation of
$C_{3}$. Following Eq. (88), the projected $C_{3}$ symmetry operator is
$D_{\text{eff}}[C_{3}]=\text{diag}(\omega^{*},\omega)$ for $n=5$. Thus,

$\eta_{\tilde{\Gamma}_{M}}=\begin{cases}\omega^{*},\quad\text{for}\quad\frac{V_{0}}{2}-2V>0,\\\
\omega,\quad\text{for}\quad\frac{V_{0}}{2}-2V<0.\end{cases}$ (91)

At $\tilde{K}_{M}$ (i.e., $\mathbf{k}=\mathbf{q}_{1}$), there are three
closest reciprocal lattice vectors, namely
$\mathbf{G}=0,\mathbf{q}_{1}-\mathbf{q}_{2},\mathbf{q}_{1}-\mathbf{q}_{3}$,
see Fig. 19(b). Other reciprocal points like $\mathbf{q}_{2}-\mathbf{q}_{1}$
and $\mathbf{q}_{3}-\mathbf{q}_{1}$ are at larger distances and will be
neglected in our first-order degenerate perturbation analysis. Notice that
$\tilde{\mathbf{K}}_{M}-(\mathbf{q}_{1}-\mathbf{q}_{2})=C_{3}\tilde{\mathbf{K}}_{M}$
and
$\tilde{\mathbf{K}}_{M}-(\mathbf{q}_{1}-\mathbf{q}_{3})=C^{2}_{3}\tilde{\mathbf{K}}_{M}$,
we obtain the following 6 $\times$ 6 Hamiltonian

$H^{\text{eff}}_{n}(\tilde{\mathbf{K}}_{M})=\begin{pmatrix}h^{\text{eff}}_{n}(\tilde{\mathbf{K}}_{M})+V_{0}\sigma_{11}&V_{1}e^{i\psi_{\xi}}\sigma_{11}&V_{1}e^{-i\psi_{\xi}}\sigma_{11}\\\
V_{1}e^{-i\psi_{\xi}}\sigma_{11}&h^{\text{eff}}_{n}(C_{3}\tilde{\mathbf{K}}_{M})+V_{0}\sigma_{11}&V_{1}e^{i\psi_{\xi}}\sigma_{11}\\\
V_{1}e^{i\psi_{\xi}}\sigma_{11}&V_{1}e^{-i\psi_{\xi}}\sigma_{11}&h^{\text{eff}}_{n}(C^{2}_{3}\tilde{\mathbf{K}}_{M})+V_{0}\sigma_{11}\end{pmatrix},\qquad\sigma_{11}\equiv(\sigma_{0}+\sigma_{3})/2\
.$ (92)

The $C_{3}$ operator that commutes with
$H^{\text{eff}}_{n}(\tilde{\mathbf{K}}_{M})$ takes the form

$\begin{pmatrix}0&0&1\\\ 1&0&0\\\ 0&1&0\end{pmatrix}\otimes
D_{\text{eff}}[C_{3}],$ (93)

which can be used to diagonalize $H^{\text{eff}}_{n}$ into three 2 $\times$ 2
blocks that correspond to symmetry eigenvalues $\omega^{j}$ ($j=0,1,2$). The
highest energy in each of these symmetry sectors are found to be

$\displaystyle\epsilon^{(c)}_{j}(\tilde{K}_{M})$
$\displaystyle=f_{0}+V_{1}\cos(\psi_{\xi}-(j+1)\frac{2\pi}{3})+\sqrt{f^{2}_{+}+\Big{(}\frac{V_{0}}{2}+V_{3}(|\tilde{K}_{M}|)+V_{1}\cos(\psi_{\xi}-(j+1)\frac{2\pi}{3})\Big{)}^{2}}$
(94a)

where $f_{0}=H_{0}(|\tilde{K}_{M}|)+V_{0}/2$, see Eq. (80), and we define

$f_{\pm}^{2}=\mathinner{\\!\left\lvert\tilde{K}_{M}\right\rvert}^{4}\big{[}\mathinner{\\!\left\lvert\tilde{K}_{M}\right\rvert}^{6}\beta^{2}+(|\tilde{K}_{M}|^{2}\gamma+\delta)^{2}\pm
2\mathinner{\\!\left\lvert\tilde{K}_{M}\right\rvert}^{3}\beta(|\tilde{K}_{M}|^{2}\gamma+\delta)\cos
3\theta_{K}\big{]},\quad\text{with}\quad\theta_{K}=\tan^{-1}\frac{\tilde{K}_{M,y}}{\tilde{K}_{M,x}}.$
(95)

Notice that $f_{\pm}$ is always non-negative. A similar analysis for the
$\tilde{K}_{M}^{\prime}$ point leads to

$\displaystyle\epsilon^{(c)}_{j}(\tilde{K}_{M}^{\prime})$
$\displaystyle=f_{0}+V_{1}\cos(\psi_{\xi}+(j+1)\frac{2\pi}{3})+\sqrt{f^{2}_{-}+\Big{(}\frac{V_{0}}{2}+V_{3}(|\tilde{K}_{M}^{\prime}|)+V_{1}\cos(\psi_{\xi}+(j+1)\frac{2\pi}{3})\Big{)}^{2}}$
(96a)

We consider certain parameter regime where it is easy to identify the energy
of the lowest conduction band at $\tilde{K}_{M}$ and $\tilde{K}_{M}^{\prime}$:
(I) with negative (or small positive) displacement field such that
$\frac{V_{0}}{2}+V_{3}(|\tilde{K}_{M}|)>0$, it is clear that the lowest
conduction state is identified by
$\min\\{\cos(\psi_{\xi}-2\pi/3),\cos(\psi_{\xi}+2\pi/3),\cos\psi_{\xi}\\}$;
(II) with very large and positive displacement field where
$|\frac{V_{0}}{2}+V_{3}(|\tilde{K}_{M}|)|\gg f_{\pm},V_{1}$, one can Taylor
expand the square root and realize that the lowest conduction state
corresponds to maximizing
$|\frac{V_{0}}{2}+V_{3}+V_{1}\cos(\psi_{\xi}+\frac{2\pi j}{3})|$. Since
$\frac{V_{0}}{2}+V_{3}$ acquires a large negative value in this case, the
lowest conduction state again corresponds to choosing
$\min\\{\cos(\psi_{\xi}-2\pi/3),\cos(\psi_{\xi}+2\pi/3),\cos\psi_{\xi}\\}$.
For reference, using SK values of pristine graphene parameters, we have
$f_{\pm}\sim 20$ meV for experimental twist angle $\theta=0.77^{\circ}$, while
$V_{0},V_{1}\sim 5$ meV as determined by our band-structure fitting analysis.

Eq. (76) implies $\psi_{\xi=0}=223.5^{\circ}$, while
$\psi_{\xi=1}=16.5^{\circ}$, which suggest

$\eta_{\text{K}}=\begin{cases}\omega^{*},\quad\text{for}\quad\xi=0\\\
\omega,\quad\text{for}\quad\xi=1\end{cases}\quad\text{and}\quad\quad\eta_{\text{K}^{\prime}}=\begin{cases}\omega^{*},\quad\text{for}\quad\xi=0\\\
1,\quad\text{for}\quad\xi=1\end{cases}.$ (97)

Irrespective of the value of $\psi_{\xi}$ (and hence the stacking
configuration), we consistently have
$\eta_{K}\eta_{\text{K}^{\prime}}=\omega$. Altogether, using the relationship
between point-group symmetry eigenvalues and Chern number [52], we find

$\exp(\frac{2i\pi
C_{c}}{3})=\begin{cases}1,\quad\text{for}\quad\frac{V_{0}}{2}-2V>0,\\\
\omega^{*},\quad\text{for}\quad\frac{V_{0}}{2}-2V<0.\end{cases}$ (98)

Furthermore, for very large positive displacement field where
$-(\frac{V_{0}}{2}+V_{3}(|\tilde{K}_{M}|))\gg f_{\pm}$ and $V_{1}$, the three
highest energies in the conduction band
$(\epsilon^{(c)}_{0},\epsilon^{(c)}_{1},\epsilon^{(c)}_{2})$ would stick
together, which corresponds to band-folding with an “empty” moiré superlattice
(empty lattice approximation). This can be understood from the fact that the
moiré potential provided by the bottom hBN has negligible effects on the
conduction band electrons localized on the top layer (see Fig. 5).

##  Appendix B Lattice structure and Structural relaxation

There are two types of structures R$n$G/hBN. The two structures are
distinguished by their stackings in the AA regions of the moiré structure.
(See the example of zero twist in Fig. 20(a),(b).) In one structure, the
carbon A/B in the bottom layer of graphene is aligned with nitrogen/boron in
the AA region, which we call $\xi=1$ configuration; the other $\xi=0$
configuration corresponds to that the carbon A/B in the bottom layer of
graphene align with boron/nitrogen. The symmetries of the $\xi=0,1$ structures
are reduced compared to the pristine R$n$G as hBN breaks inversion and nonzero
small twist angles breaks mirrors. Only the three-fold rotation $C_{3}$ (with
axis perpendicular to the sample) can be preserved for all twist angles. The
$\xi=0,1$ structures are related by a $C_{2}$ rotation (about the axis
perpendicular to the sample) on the hBN only.

We also consider hBN/R$n$G/hBN structures with hBN nearly aligned on top and
on bottom. We choose the alignments of two hBNs such that there is only one
moiré pattern and $C_{3}$ symmetry (with axis perpendicular to the sample) is
preserved. The stackings of the top and bottom hBNs (relative to graphene) are
labeled by $\xi_{t},\xi_{b}$, respectively. For simplicity, we in this work
only study $(\xi_{b},\xi_{t})=(1,1)$ (Fig. 20(c)) and
$(\xi_{b},\xi_{t})=(1,0)$ (Fig. 20(d)). We do not need to study
$(\xi_{b},\xi_{t})=(0,0)$, since it is related to $(\xi_{b},\xi_{t})=(1,1)$ by
a $C_{2}$ rotation. We note that $(\xi_{b},\xi_{t})=(1,0),(0,1)$ can restore
the inversion symmetry if the two hBNs are related by inversion and their
inversion center is aligned with that of the R$n$G; we choose
$(\xi_{b},\xi_{t})=(1,0)$ as a representative as this type of structure. In
our DFT calculation for $(\xi_{b},\xi_{t})=(1,0)$ structure, the two inversion
centers are nearly aligned with inversion symmetry preserved up to energy
error about 0.1meV, which will be discussed in App. C. Certainly, inversion
symmetry can be broken by the displacement field.

Figure 20:  Rigid R$5$G/hBN and hBN/R$5$G/hBN lattice structures with zero
twist. (a),(b) are pentalayer graphene stack on hBN, (a) with a carbon atom
strictly colinear with a boron atom in bottom layer along z-axis, (b) with a
carbon atom strictly colinear with a nitrogen atom in bottom layer along
z-axis. (c),(d) are pentalayer graphene between two hBN layers. (a) and (b)
are connected by applying $C_{2}$ rotation (with axis along z-axis) on hBN
only, (d) has an approximate inversion symmetry while (c) breaks the inversion
symmetry stronly.

The rigid lattice constant for Graphene is 0.246nm and for hBN is 0.25nm,
rigid layer distance for both Graphene-Graphene and Graphene-hBN are 0.336nm.
Our calculations mainly focus on the $0.77^{\circ}$ twisted angle, which the
moiré lattice constant is approximately 11.6nm, and the number of carbon atom
in the moiré superlattice is $4472$ per layer, the number of boron and
nitrogen atom in the moiré superlattice are both $2163$ per layer.

We perform classical structure relaxation implemented in LAMMPS[33]. During
the relaxation, we fix each hBN layer to simulate a thick substrate and keep
the moiré unit cell unchanged. We used empirical inter-atomic potentials in
LAMMPS to perform relaxation. In graphene and hBN systems, the inter-layer
interaction acts differently than intra-layer one due to the van der Waals
interaction. Therefore, they are usually treated with different empirical
potentials in classical molecular dynamics simulations. For intra-layer
interactions within graphene layers, we used the reactive empirical bond order
potential[34]. For inter-layer interactions, we used an inter-layer potential
developed for graphene and hBN systems[35]. This combination of empirical
potentials has been well tested in Ref. [35], and have good agreement with DFT
results.

Figure 21: Relaxation results of 0.767 ∘ R5G/hBN${}_{\xi_{b}=1}$. $(a)$ Rigid
structure of R5G/hBN${}_{\xi_{b}=1}$ for $\xi=1$. $(b)$-$(f)$ show the intra-
layer displacement of each layer. $(g)$-$(k)$ show the inter-layer distance of
each layer. The Gr0 layer is the lowest layer, which lies on the hBN
substrate. The intra-layer displacement is the in-plane displacement from
rigid position to relaxed position of an atom in the corresponding layer. The
inter-layer distance indicates the distance between the corresponding layer
and the layer under it. Figure 22: Relaxation results of 0.767 ∘
R5G/hBN${}_{\xi_{b}=0}$. (a), structure of R5G/hBN${}_{\xi_{b}=0}$. (b)-(f),
intra-layer displacement of each layer. (g)-(k), inter-layer distance of each
layer. The Gr0 layer is the lowest layer, which lies on the hBN substrate.

The relaxation results for R$n$G/hBN configurations are listed in Fig. 21 and
Fig. 22. The interlayer distance between the Gr0 layer and hBN becomes the
largest(3.52Å) in the AA region and the smallest(3.29Å) in the AB region. In
upper graphene layers, the inter-layer distance between graphene is around
3.33Å in all regions. On the other hand, the magnitude of intra-layer
displacement also decreases with the layer number, which is reasonable since
the hBN substrate has less effect on the upper layers. The atoms near the AB
region tend to rotate in a clockwise direction, against the global twist
(counterclockwise). This will enlarge the AB region, where the local stacking
energy reaches the minimum.

The relaxation results for hBN/R$5$G/hBN are also listed in Fig. 23 and Fig.
24. In these configurations, the distance between upper and lower hBN is fixed
at 22.66Å.

Figure 23: Relaxation results of 0.767 ∘
hBN${}_{\xi_{t}=1}$/R5G/hBN${}_{\xi_{b}=1}$. (a), structure of
hBN${}_{\xi_{t}=1}$/R5G/hBN${}_{\xi_{b}=1}$. (b)-(f), intra-layer displacement
of each layer. (g)-(k), inter-layer distance of each layer. The Gr0 layer is
the lowest, and Gr4 is the highest layer. The Gr0 and Gr4 layers are near hBN
substrates. Figure 24: Relaxation results of 0.767 ∘
hBN${}_{\xi_{t}=0}$/R5G/hBN${}_{\xi_{b}=1}$. (a), structure of
hBN${}_{\xi_{t}=0}$/R5G/hBN${}_{\xi_{b}=1}$. (b)-(f), intra-layer displacement
of each layer. (g)-(h), inter-layer distance of each layer. The Gr0 layer is
the lowest layer, and Gr4 is the highest layer. The Gr0 and Gr4 layers are
near hBN substrates.

##  Appendix C Slater-Koster method

The band structure are obtained within the SK tight-binding model using the
parameters listed in Eq. (39). The diagonalization is performed with an open-
source software WannierTools[110]. In our study, we consider an internal
symmetrical polarization (ISP) as 5mV/Å and $V_{\mathrm{pp}\pi}^{0}$ as -2.81
eV [111], a value exceeding the conventional -2.7 eV. This adjustment is made
to enhance the compatibility with our DFT calculations of ABC-stacked graphene
with different number layers in Fig. 25.

Our study delves into the distinctive characteristics of the low-energy
valence and conduction bands. By comparing the relaxed band structure and the
rigid band structure in Figs. 26, 27, 28, 29 and 30, we can find significant
changes in the band structure in some specific cases whose gap decrease at
charge neutrality. In order to delve deeper into this pattern, we
systematically investigated the flatter bandwidths across different numbers of
graphene layers, ranging from 3 to 7. The outcomes of our analysis are
succinctly presented in Figs. 32, 33, 34, 35 and 36 whose panels (a)-(t)
showcase the band structures of R$n$G/$\mathrm{hBN}_{\xi_{b}=1}$,
$\mathrm{hBN}_{\xi_{t}=1}$/R$n$G/$\mathrm{hBN}_{\xi_{b}=1}$,
R$n$G/$\mathrm{hBN}_{\xi_{b}=0}$ and
$\mathrm{hBN}_{\xi_{t}=0}$/R$n$G/$\mathrm{hBN}_{\xi_{b}=1}$ respectively.
Applying displacement field would help generate flat bands.

Figure 25:  Comparison of the tight-binding Hamiltonian (red circles) and ab
initio DFT (black lines) band structure calculations along the
$\Gamma-\text{K}-\text{M}$ path for ABC-stacked graphene: (a) Illustration of
the $\gamma_{2}$ bond correction in the SK tight-binding model. Panels (b) and
(c) for trilayer, (d) and (e) for tetralayer, (f) and (g) for pentalayer, (h)
and (i) for hexalayer, and (j) and (k) for septalayer graphene depict band
structures near the K point within two energy windows, from -0.1 eV to 0.1 eV
and from -1 eV to 1 eV, using a set value of $V_{\mathrm{pp}\pi}^{0}=-2.81$ eV
and an internal symmetrical polarization (ISP) of 5mV/Å.  Figure 26:
Comparison of band structures between rigid and relaxed structures of R3G and
0.76715∘ twisted angle setting ISP=5mV/Å. E is an applied electrical ﬁeld and
ISP is an internal symmetrical polarization due to the diﬀerent chemical
environment of outer and inner atoms in trilayer graphene. The positive
direction of E and ISP are shown in Fig. 2. The structures corresponding to
the first row to the fourth row are R3G/$\mathrm{hBN}_{\xi_{b}=1}$,
$\mathrm{hBN}_{\xi_{t}=1}$/R3G/$\mathrm{hBN}_{\xi_{b}=1}$,
R3G/$\mathrm{hBN}_{\xi_{b}=0}$ and
$\mathrm{hBN}_{\xi_{t}=0}$/R3G/$\mathrm{hBN}_{\xi_{b}=1}$ respectively, and
the applied electrical fields corresponding to the first column to the fifth
column are -5, -2 0, 2 and 5mV/Å respectively. The band structures are
depicted with orange dotted lines for rigid structures and black dotted lines
for relaxed structures. Here, both $K$ valley and $K^{\prime}$ valley bands
are included. Figure 27: Comparison of band structures between rigid and
relaxed structures of R4G and 0.76715∘ twisted angle setting ISP=5mV/Å. E is
an applied electrical ﬁeld and ISP is an internal symmetrical polarization due
to the diﬀerent chemical environment of outer and inner atoms in tetralayer
graphene. The positive direction of E and ISP are shown in Fig. 2. The
structures corresponding to the first row to the fourth row are
R4G/$\mathrm{hBN}_{\xi_{b}=1}$,
$\mathrm{hBN}_{\xi_{t}=1}$/R4G/$\mathrm{hBN}_{\xi_{b}=1}$,
R4G/$\mathrm{hBN}_{\xi_{b}=0}$ and
$\mathrm{hBN}_{\xi_{t}=0}$/R4G/$\mathrm{hBN}_{\xi_{b}=1}$ respectively, and
the applied electrical fields corresponding to the first column to the fifth
column are -5, -2 0, 2 and 5mV/Å respectively. The band structures are
depicted with orange dotted lines for rigid structures and black dotted lines
for relaxed structures. Here, both $K$ valley and $K^{\prime}$ valley bands
are included. Figure 28: Comparison of band structures between rigid and
relaxed structures of R5G and 0.76715∘ twisted angle setting ISP=5mV/Å. E is
an applied electrical ﬁeld and ISP is an internal symmetrical polarization due
to the diﬀerent chemical environment of outer and inner atoms in pentalayer
graphene. The positive direction of E and ISP are shown in Fig. 2. The
structures corresponding to the first row to the fourth row are
R5G/$\mathrm{hBN}_{\xi_{b}=1}$,
$\mathrm{hBN}_{\xi_{t}=1}$/R5G/$\mathrm{hBN}_{\xi_{b}=1}$,
R5G/$\mathrm{hBN}_{\xi_{b}=0}$ and
$\mathrm{hBN}_{\xi_{t}=0}$/R5G/$\mathrm{hBN}_{\xi_{b}=1}$ respectively, and
the applied electrical fields corresponding to the first column to the fifth
column are -5, -2 0, 2 and 5mV/Å respectively. The band structures are
depicted with orange dotted lines for rigid structures and black dotted lines
for relaxed structures. Here, both $K$ valley and $K^{\prime}$ valley bands
are included. Figure 29: Comparison of band structures between rigid and
relaxed structures of R6G and 0.76715∘ twisted angle setting ISP=5mV/Å. E is
an applied electrical ﬁeld and ISP is an internal symmetrical polarization due
to the diﬀerent chemical environment of outer and inner atoms in hexalayer
graphene. The positive direction of E and ISP are shown in Fig. 2. The
structures corresponding to the first row to the fourth row are
R6G/$\mathrm{hBN}_{\xi_{b}=1}$,
$\mathrm{hBN}_{\xi_{t}=1}$/R6G/$\mathrm{hBN}_{\xi_{b}=1}$,
R6G/$\mathrm{hBN}_{\xi_{b}=0}$ and
$\mathrm{hBN}_{\xi_{t}=0}$/R6G/$\mathrm{hBN}_{\xi_{b}=1}$ respectively, and
the applied electrical fields corresponding to the first column to the fifth
column are -5, -2 0, 2 and 5mV/Å respectively. The band structures are
depicted with orange dotted lines for rigid structures and black dotted lines
for relaxed structures. Here, both $K$ valley and $K^{\prime}$ valley bands
are included. Figure 30: Comparison of band structures between rigid and
relaxed structures of R7G and 0.76715∘ twisted angle setting ISP=5mV/Å. E is
an applied electrical ﬁeld and ISP is an internal symmetrical polarization due
to the diﬀerent chemical environment of outer and inner atoms in septalayer
graphene. The positive direction of E and ISP are shown in Fig. 2. The
structures corresponding to the first row to the fourth row are
R7G/$\mathrm{hBN}_{\xi_{b}=1}$,
$\mathrm{hBN}_{\xi_{t}=1}$/R7G/$\mathrm{hBN}_{\xi_{b}=1}$,
R7G/$\mathrm{hBN}_{\xi_{b}=0}$ and
$\mathrm{hBN}_{\xi_{t}=0}$/R7G/$\mathrm{hBN}_{\xi_{b}=1}$ respectively, and
the applied electrical fields corresponding to the first column to the fifth
column are -5, -2 0, 2 and 5mV/Å respectively. The band structures are
depicted with orange dotted lines for rigid structures and black dotted lines
for relaxed structures. Here, both $K$ valley and $K^{\prime}$ valley bands
are included.

##  Appendix D Valley-Resolved Band Structure

In superlattices form by R$n$G and hBN, the valleys of R$n$G emerge as a good
quantum number in the low energy part of bands. Inherited from the gapless
Dirac dispersion band structure at $\text{K},\text{K}^{\prime}$ valleys in
untwisted case, the low energy bands of R$n$G-hBN superlattice either belongs
to $K$ valley or $K^{\prime}$ valley. Thus, We can construct an operator to
distinguish the valley degree of freedom with eigenvalue $+1$ for the states
at one valley and eigenvalue $-1$ for the states at the other valley, which is
in the same spirit as the spin-z operator in basis
$\\{\ket{\uparrow},\ket{\downarrow}\\}$.

The valley operator for Dirac fermion model — single layer graphene(SLG) in
real space reads [38]

$\mathcal{V}_{z}^{Gr}=\frac{i}{3\sqrt{3}}\sum_{\ll
i,j\gg,s}\eta_{ij}\sigma_{z}^{i\,j}c^{\dagger}_{i,s}c_{j,s}$ (99)

where $\ll i,j\gg$ denotes the next-nearest-neighbor(NNN) hopping term, $s$
denote spin, $\eta_{ij}=\pm 1$ for clockwise or counterclockwise hopping and
$\sigma_{z}$ is in the sublattice degree of freedom, meaning that
$\sigma_{z}^{i\,j}=1$ if both $i,j$ are on the $A$ sublattice,
$\sigma_{z}^{i\,j}=-1$ if both $i,j$ are on the $B$ sublattice and
$\sigma_{z}^{i\,j}=0$ otherwise. The schematic diagram of valley operator in
real space for an unit cell is shown in Fig. 31(a). The factor
$\eta_{ij}\sigma_{z}^{ij}$ forms the valley flux, which distinguishes the NNN
hopping in $K$ direction and $K^{\prime}$ direction. To specify the valley
flux, we compare it with Haldane’s local magnetic flux(Fig. 31(b),(c)). While
the latter holds the same flux with both sublattice A and B, valley flux means
sublattice A/B holds opposite local magnetic flux. After performing Fourier
transformation, the valley operator of single layer graphene in sublattice is
proportional to the identity matrix in reciprocal space,
$\mathcal{V}_{z}^{Gr}(k)=\left(\begin{matrix}f(k)&\\\
&f(k)\end{matrix}\right)$. $f(k)$ is an odd function in Brillouin Zone
$f(k)=-f(-k)$, which satisfies
$\bra{\psi(K)}\mathcal{V}_{z}(K)\ket{\psi(K)}=1$ and
$\bra{\psi(K^{\prime})}\mathcal{V}_{z}(K^{\prime})\ket{\psi(K^{\prime})}=-1$(Fig.
31(d),(e)).

Figure 31:  Valley operator for single layer graphene. (a) shows the valley
operator in real space for one unit cell; (b) shows the Haldane’s NNN hopping
term with local magnetic flux and (c) shows the corresponding NNN hopping term
with valley flux; (d) shows the average value of valley operator
$\bra{\psi(k)}\mathcal{V}_{z}^{Gr}(k)\ket{\psi(k)}$ in high symmetry point
path of reciprocal space. (e) is the energy band of SLG follow the same path
as (d), weighted by $\bra{\psi(k)}\mathcal{V}_{z}^{Gr}(k)\ket{\psi(k)}$. (f)
is the flow diagram of valley resolved bands calculation, this flow is
applicable for both SLG and hBN-R$n$G.

Since valley operator in real space only contains intra-layer hopping term, it
can be easily generalised to R$n$G by sum up valley operator of each layer
together $\mathcal{V}_{z}=\sum_{l}\mathcal{V}_{z}^{l}$, where $l$ denotes the
layer index. And for twisted hBN-R$n$G system, valley operator of each layer
can be obtained numerically by considering all the NNN hopping term multiplied
by valley flux factor $\eta_{ij}\sigma_{z}^{ij}$ for the atoms in superlattice
and this layer. The process to obtain real space valley operator for moiré
system corresponds to valley band folding in momentum space.

To identify the valley quantum number of the energy bands in the momentum
space, we solve the valley operator of hBN-R$n$G $\mathcal{V}_{z}$ in
reciprocal space. Using the average value of $\mathcal{V}_{z}(k)$ under the
eigenvectors of $H(k)$, we can distinguish the valley components of each
state. And in the low energy region, it satisfies
$\bra{\psi(k)}\mathcal{V}_{z}(k)\ket{\psi(k)}\simeq 1$ for $K$ valley and
$\bra{\psi(k)}\mathcal{V}_{z}(k)\ket{\psi(k)}\simeq-1$ for $K^{\prime}$
valley. Based on the eigenvalue of each state, we can make a bipartition for
bands with different valley components and obtain the valley-resolved bands.
The calculation steps for valley-resolved bands are shown in Fig. 31(f). The
calculations are implemented in WannierTools package[110]. The valley-resolved
bands of hBN-R$n$G(n=3,4,5,6,7) are shown in Figs. 32, 33, 34, 35 and 36.

Figure 32: Tight-binding band structures of R3G and 0.76715∘ twist angle
setting ISP=$5$meV/Å. E is an applied electrical ﬁeld, and ISP is an internal
symmetrical polarization due to the diﬀerent chemical environment of outer and
inner atoms in trilayer graphene. The positive direction of E and ISP are
shown in Fig. 2. The structures corresponding to the first row to the fourth
row are R3G/$\mathrm{hBN}_{\xi_{b}=1}$,
$\mathrm{hBN}_{\xi_{t}=1}$/R3G/$\mathrm{hBN}_{\xi_{b}=1}$,
R3G/$\mathrm{hBN}_{\xi_{b}=0}$ and
$\mathrm{hBN}_{\xi_{t}=0}$/R3G/$\mathrm{hBN}_{\xi_{b}=1}$ respectively, and
the applied electrical fields corresponding to the first column to the fifth
column are -5, -2 0, 2 and 5mV/Å respectively. The band structures are
depicted with blue dotted lines for $K$ valley and red dotted lines for
$K^{\prime}$ valley.  Figure 33: Tight-binding band structures of R4G and
0.76715∘ twisted angle setting ISP=5mV/Å. E is an applied electrical ﬁeld and
ISP is an internal symmetrical polarization due to the diﬀerent chemical
environment of outer and inner atoms in tetralayer graphene. The positive
direction of E and ISP are shown in Fig. 2. The structures corresponding to
the first row to the fourth row are R4G/$\mathrm{hBN}_{\xi_{b}=1}$,
$\mathrm{hBN}_{\xi_{t}=1}$/R4G/$\mathrm{hBN}_{\xi_{b}=1}$,
R4G/$\mathrm{hBN}_{\xi_{b}=0}$ and
$\mathrm{hBN}_{\xi_{t}=0}$/R4G/$\mathrm{hBN}_{\xi_{b}=1}$ respectively, and
the applied electrical fields corresponding to the first column to the fifth
column are -5, -2, 0, 2 and 5mV/Å respectively. The band structures are
depicted with blue dotted lines for $K$ valley and red dotted lines for
$K^{\prime}$ valley.  Figure 34: Tight-binding band structures of R5G and
0.76715∘ twisted angle setting ISP=5mV/Å. E is an applied electrical ﬁeld and
ISP is an internal symmetrical polarization due to the diﬀerent chemical
environment of outer and inner atoms in pentalayer graphene. The positive
direction of E and ISP are shown in Fig. 2. The structures corresponding to
the first row to the fourth row are R5G/$\mathrm{hBN}_{\xi_{b}=1}$,
$\mathrm{hBN}_{\xi_{t}=1}$/R5G/$\mathrm{hBN}_{\xi_{b}=1}$,
R5G/$\mathrm{hBN}_{\xi_{b}=0}$ and
$\mathrm{hBN}_{\xi_{t}=0}$/R5G/$\mathrm{hBN}_{\xi_{b}=1}$ respectively and the
applied electrical fields corresponding to the first column to the fifth
column are -5, -2, 0, 2 and 5mV/Å respectively. The band structures are
depicted with blue dotted lines for $K$ valley and red dotted lines for
$K^{\prime}$ valley.  Figure 35: Tight-binding band structures of R6G and
0.76715∘ twisted angle setting ISP=5mV/Å. E is an applied electrical ﬁeld and
ISP is an internal symmetrical polarization due to the diﬀerent chemical
environment of outer and inner atoms in hexalayer graphene. The positive
direction of E and ISP are shown in Fig. 2. The structures corresponding to
the first row to the fourth row are R6G/$\mathrm{hBN}_{\xi_{b}=1}$,
$\mathrm{hBN}_{\xi_{t}=1}$/R6G/$\mathrm{hBN}_{\xi_{b}=1}$,
R6G/$\mathrm{hBN}_{\xi_{b}=0}$ and
$\mathrm{hBN}_{\xi_{t}=0}$/R6G/$\mathrm{hBN}_{\xi_{b}=1}$ respectively, and
the applied electrical fields corresponding to the first column to the fifth
column are -5, -2, 0, 2 and 5mV/Å respectively. The band structures are
depicted with blue dotted lines for $K$ valley and red dotted lines for
$K^{\prime}$ valley.  Figure 36: Tight-binding band structures of R7G and
0.76715∘ twisted angle setting ISP=5mV/Å. E is an applied electrical ﬁeld and
ISP is an internal symmetrical polarization due to the diﬀerent chemical
environment of outer and inner atoms in septalayer graphene. The positive
direction of E and ISP are shown in Fig. 2. The structures corresponding to
the first row to the fourth row are R7G/$\mathrm{hBN}_{\xi_{b}=1}$,
$\mathrm{hBN}_{\xi_{t}=1}$/R7G/$\mathrm{hBN}_{\xi_{b}=1}$,
R7G/$\mathrm{hBN}_{\xi_{b}=0}$ and
$\mathrm{hBN}_{\xi_{t}=0}$/R7G/$\mathrm{hBN}_{\xi_{b}=1}$ respectively, and
the applied electrical fields corresponding to the first column to the fifth
column are -5, -2, 0, 2 and 5mV/Å respectively. The band structures are
depicted with blue dotted lines for $K$ valley and red dotted lines for
$K^{\prime}$ valley.

##  Appendix E Parameter Fitting and Phase diagrams

In this part, we discuss the fitting of the $2n\times 2n$ continuum model (Eq.
(77)) and the $2\times 2$ effective continuum model (Eq. (86)). We will also
discuss the single-particle phase diagram given by the $2n\times 2n$ continuum
model. Without loss of generality, we focus on the K valley.

We discuss the $2n\times 2n$ continuum model first. From SK method, we can
directly derive the parameter values for R$n$G ($n\geq 3$), which read

$v_{F}=542.1\text{meV}\cdot\text{nm}\ ,\ t_{1}=355.16\text{meV}\ ,\
v_{3}=34.5\text{meV}\cdot\text{nm}\ ,\ t_{2}=-7\text{meV}\ .$ (100)

Regarding to the moiré potential parameters in Eq. (72) and Eq. (76), we adopt
the values of $\psi_{\xi_{b}}$ and $\psi_{\xi_{t}}$ based on the perturbative
expression in Eq. (18), while treating $V_{b0}$, $V_{b1}$, $V_{t0}$ and
$V_{t1}$ as independent parameters.

For R$n$G/hBN structures, we fix $V_{t0}=0$ and $V_{t1}=0$, and fit $V_{b0}$
and $V_{b1}$ to the DFT+SK band structure in the K valley; according to the
effective model in Eq. (17), $|V_{b0}|$ (or $\sqrt{V_{b0}^{2}+4t_{2}^{2}}$ for
$n=3$) is just the gap of the bands at $\text{K}_{M}$, which can be directly
read out form the DFT+SK bands, and we adopt $V_{b0}>0$ owing to Eq. (18). (If
the gap at $\text{K}_{M}$ is smaller than $|2t_{2}|=14$meV, we just choose
$V_{b0}=0$.) Therefore, for single hBN, we effectively only need to optimize
one parameter which is $V_{b1}$. The parameter values after the optimization
are listed in Tab. 1, where $V_{b0}$, $V_{b1}$ and $\psi_{\xi_{b}}$ are
labeled as $V_{0}$, $V_{1}$ and $\psi$, respectively, in short. The fitting is
remarkably good as shown in Figs. 37-41.

For hBN/R$n$G/hBN structures with configuration $(\xi_{b},\xi_{t})=(1,1)$, we
have four tuning parameters $V_{b0}$, $V_{b1}$, $V_{t0}$ and $V_{t1}$.
According to the effective model in Eq. (17), $(V_{b0}+V_{t0})/2$ is just a
total energy shift for low-energy bands, while $|V_{b0}-V_{t0}|$ (or
$\sqrt{|V_{b0}-V_{t0}|^{2}+4t_{2}^{2}}$ for $n=3$) corresponds to the gap at
$\text{K}_{M}$ in the K valley; thus, we just choose $V_{b0}$ to be the
corresponding R$n$G/hBN value, and determine $V_{t0}$ from the gap. (Here we
choose $V_{t0}>V_{b0}$ since the opposite gives considerably worse fitting,
and we choose $V_{t0}=V_{b0}$ for $n=3$ if the gap at $\text{K}_{M}$ in the K
valley is larger than $|2t_{2}|=14$meV. Therefore, we need to optimize two
parameters $V_{b1}$ and $V_{t1}$ in this case. For double hBN with
configuration $(\xi_{b},\xi_{t})=(1,0)$, the extra inversion symmetry requires
$V_{b0}=V_{t0}$ and $V_{b1}=V_{t1}$. According to the effective model Eq.
(17), $V_{b0}=V_{t0}$ is just effectively a total energy shift for low-energy
bands, which we can just choose $V_{b0}$ to be the corresponding R$n$G/hBN
value. As a result, we only have one parameter to optimize which is $V_{b1}$.
The parameter values for double hBN after the optimization are listed in Tab.
3, and the fitting is good as shown in Figs. 37-41.

Based on the $2n\times 2n$ continuum model, we can plot the single-particle
phase diagrams as a function of twisted angle $\theta$ and the inter-layer
potential energy difference $V$ for all the configurations, as shown in Figs.
42-46. In the phase diagrams, we only consider the Chern numbers of lowest
conduction band and the highest valence band. In sum, we find that the
R$n$G/hBN configuration always achieve Chern numbers $0$ or $n$ when the
direct gap of those bands is larger than 1meV, while $\pm 1$ states can happen
in the hBN/R$n$G/hBN configuration with considerable direct gap ($>$2meV).

We now turn to the $2\times 2$ effective continuum model in Eq. (86). In
principle, all parameters in the $2\times 2$ effective model can be derived
from the $2n\times 2n$ model in Eq. 77, as discussed in App. A.4. However, as
shown in Fig. 10 of the Main Text, directly using the values derived from the
perturbation theory can only partly well match the low-energy features, mainly
around $\text{K}_{M}$ point. (Recall that we focus on the K valley.) To
achieve a better match to the band structure, we treat
$\\{\alpha,\beta,\gamma,V_{1b}\\}$ as tuning parameters for R$n$G/hBN
structures ($\\{\alpha,\beta,\gamma,V_{1b},V_{1t}\\}$ for hBN/R$n$G/hBN
structures), and optimize them around their values derived from the
perturbation theory, while preserving the form of the model. Other parameters
take the values derived from the $2n\times 2n$ continuum model. With those
fitting parameter values listed in Tab. 1 and Tab. 4, we can improve the
matching away from $\text{K}_{M}$ point, especially at $\Gamma_{M}$, and
capture well the low-energy features of the DFT+SK bands, as shown in Figs.
47-41.

| $\alpha$ | $\beta$ | $\gamma$ | $\delta$ | $V_{b0}$ | $V_{b1}$ | $V_{t0}$ | $V_{t1}$
---|---|---|---|---|---|---|---|---
n=3, $(\xi_{b},\xi_{t})$=(1,1) | 59.61(103.79) | 864.58(1263.00) | 81.70(120.10) | -7.00 | 0 | 0.09(3.20) | 0 | 12.03(11.09)
n=4, $(\xi_{b},\xi_{t})$=(1,1) | 55.21(103.79) | -1208.20(-1927.70) | -197.86(-287.42) | 21.37 | 1.44 | 0.92(5.76) | 5.40 | 8.87(7.08)
n=5, $(\xi_{b},\xi_{t})$=(1,1) | 49.37(103.79) | 1879.20(2942.40) | 405.76(597.60) | -48.92 | 1.50 | 6.27(7.29) | 6.48 | 6.46(7.91)
n=6, $(\xi_{b},\xi_{t})$=(1,1) | 34.99(103.79) | -3000.00(-4491.10) | -650.25(-1154.70) | 99.57 | 1.56 | 5.70(6.02) | 7.52 | 5.53(7.78)
n=7, $(\xi_{b},\xi_{t})$=(1,1) | 36.69(103.79) | 4113.40(6855.10) | 15.67(2132.70) | -189.97 | 1.47 | 4.90(5.45) | 5.79 | 5.49(7.79)
n=3, $(\xi_{b},\xi_{t})$=(1,0) | 69.25(103.79) | 870.13(1263.00) | 92.72(120.10) | -7.00 | 0 | 1.76(6.72) | 0 | 1.76(6.72)
n=4, $(\xi_{b},\xi_{t})$=(1,0) | 58.73(103.79) | -1324.00(-1927.70) | -137.21(-287.42) | 21.37 | 1.44 | 3.78(7.65) | 1.44 | 3.78(7.65)
n=5, $(\xi_{b},\xi_{t})$=(1,0) | 52.95(103.79) | 1914.20(2942.40) | 414.48(597.60) | -48.92 | 1.50 | 4.14(5.43) | 1.50 | 4.14(5.43)
n=6, $(\xi_{b},\xi_{t})$=(1,0) | 53.64(103.79) | -3357.80(-4491.10) | -89.39(-1154.70) | 99.57 | 1.56 | 4.78(7.80) | 1.56 | 4.78(7.80)
n=7, $(\xi_{b},\xi_{t})$=(1,0) | 48.46(103.79) | 4619.20(6855.10) | 1876.40(2132.70) | -189.97 | 1.47 | 5.11(7.22) | 1.47 | 5.11(7.22)
Table 4: Parameter values of the $2\times 2$ effective moiré model (Eq. (86))
for the hBN/R$n$G/hBN structure of $n=3,4,5,6,7$ layers. Here
$\alpha,\beta,\gamma,\delta$ are reported in meV$\cdot$nm2,meV$\cdot$nmn,
meV$\cdot$nmn-1 and meV, respectively, while $V_{b0}$, $V_{b1}$, $V_{t0}$ and
$V_{t1}$ are in meV. If the parameter values are difference from those
directly derived from the perturbation theory of $2n$-band model, the latter
will be provided in the brackets. The values of $\psi_{t}$ and $\psi_{b}$ are
the same as those in Tab. 3; we omit those values here due to the length
limit.  Figure 37:  The comparison between the DFT+SK bands (black) and the
bands from the $2n\times 2n$ continuum model (blue lines) in Eq. (77) for
$n=3$. The parameter values are listed in Tab. 1 and Tab. 3 of Main Text.
Figure 38:  The comparison between the DFT+SK bands (black) and the bands from
the $2n\times 2n$ continuum model (blue lines) in Eq. (77) for $n=4$. The
parameter values are listed in Tab. 1 and Tab. 3 of Main Text.  Figure 39:
The comparison between the DFT+SK bands (black) and the bands from the
$2n\times 2n$ continuum model (blue lines) in Eq. (77) for $n=5$. The
parameter values are listed in Tab. 1 and Tab. 3 of Main Text.  Figure 40:
The comparison between the DFT+SK bands (black) and the bands from the
$2n\times 2n$ continuum model (blue lines) in Eq. (77) for $n=6$. The
parameter values are listed in Tab. 1 and Tab. 3 of Main Text.  Figure 41:
The comparison between the DFT+SK bands (black) and the bands from the
$2n\times 2n$ continuum model (blue lines) in Eq. (77) for $n=7$. The
parameter values are listed in Tab. 1 and Tab. 3 of Main Text.  Figure 42:
Phase diagrams of R3G-hBN superlattices of various stacking configurations:
(a) $\xi=1$ , (b) $\xi=0$, (c) $(\xi_{b},\xi_{t})=(1,1)$, and (d)
$(\xi_{b},\xi_{t})=(1,0)$. Panel (i) shows the single-particle indirect gap
$\Delta_{\text{ind},\nu=-4}$ at filling $\nu=-4$; (ii) shows the indirect gap
$\Delta_{\text{ind},\nu=0}$ at filling $\nu=0$; (iii) shows the indirect gap
$\Delta_{\text{ind},\nu=4}$ at filling $\nu=4$; (iv) shows the minimal direct
gap $\Delta_{\text{val}}$ around the highest valence band in one valley; (v)
shows the bandwidth $W_{\text{val}}$ of the highest valence band; (vi) shows
the minimal direct gap $\Delta_{\text{cond}}$ around the lowest conduction
band in one valley; (vii) shows the bandwidth $W_{\text{cond}}$ of the lowest
conduction band. Chern numbers of the highest valence band and the lowest
conduction band in K valley are indicated on panel (iv) and panel (vi),
respectively, where boundaries of topologically distinct phases can be seen as
the direct gap closes.  Figure 43: Phase diagrams of R4G-hBN superlattices of
various stacking configurations: (a) $\xi=1$ , (b) $\xi=0$, (c)
$(\xi_{b},\xi_{t})=(1,1)$, and (d) $(\xi_{b},\xi_{t})=(1,0)$. Panel (i) shows
the single-particle indirect gap $\Delta_{\text{ind},\nu=-4}$ at filling
$\nu=-4$; (ii) shows the indirect gap $\Delta_{\text{ind},\nu=0}$ at filling
$\nu=0$; (iii) shows the indirect gap $\Delta_{\text{ind},\nu=4}$ at filling
$\nu=4$; (iv) shows the minimal direct gap $\Delta_{\text{val}}$ around the
highest valence band in one valley; (v) shows the bandwidth $W_{\text{val}}$
of the highest valence band; (vi) shows the minimal direct gap
$\Delta_{\text{cond}}$ around the lowest conduction band in one valley; (vii)
shows the bandwidth $W_{\text{cond}}$ of the lowest conduction band. Chern
numbers of the highest valence band and the lowest conduction band in K valley
are indicated on panel (iv) and panel (vi), respectively, where boundaries of
topologically distinct phases can be seen as the direct gap closes.  Figure
44: Phase diagrams of R5G-hBN superlattices of various stacking
configurations: (a) $\xi=1$ , (b) $\xi=0$, (c) $(\xi_{b},\xi_{t})=(1,1)$, and
(d) $(\xi_{b},\xi_{t})=(1,0)$. Panel (i) shows the single-particle indirect
gap $\Delta_{\text{ind},\nu=-4}$ at filling $\nu=-4$; (ii) shows the indirect
gap $\Delta_{\text{ind},\nu=0}$ at filling $\nu=0$; (iii) shows the indirect
gap $\Delta_{\text{ind},\nu=4}$ at filling $\nu=4$; (iv) shows the minimal
direct gap $\Delta_{\text{val}}$ around the highest valence band in one
valley; (v) shows the bandwidth $W_{\text{val}}$ of the highest valence band;
(vi) shows the minimal direct gap $\Delta_{\text{cond}}$ around the lowest
conduction band in one valley; (vii) shows the bandwidth $W_{\text{cond}}$ of
the lowest conduction band. Chern numbers of the highest valence band and the
lowest conduction band in K valley are indicated on panel (iv) and panel (vi),
respectively, where boundaries of topologically distinct phases can be seen as
the direct gap closes.  Figure 45: Phase diagrams of R6G-hBN superlattices of
various stacking configurations: (a) $\xi=1$ , (b) $\xi=0$, (c)
$(\xi_{b},\xi_{t})=(1,1)$, and (d) $(\xi_{b},\xi_{t})=(1,0)$. Panel (i) shows
the single-particle indirect gap $\Delta_{\text{ind},\nu=-4}$ at filling
$\nu=-4$; (ii) shows the indirect gap $\Delta_{\text{ind},\nu=0}$ at filling
$\nu=0$; (iii) shows the indirect gap $\Delta_{\text{ind},\nu=4}$ at filling
$\nu=4$; (iv) shows the minimal direct gap $\Delta_{\text{val}}$ around the
highest valence band in one valley; (v) shows the bandwidth $W_{\text{val}}$
of the highest valence band; (vi) shows the minimal direct gap
$\Delta_{\text{cond}}$ around the lowest conduction band in one valley; (vii)
shows the bandwidth $W_{\text{cond}}$ of the lowest conduction band. Chern
numbers of the highest valence band and the lowest conduction band in K valley
are indicated on panel (iv) and panel (vi), respectively, where boundaries of
topologically distinct phases can be seen as the direct gap closes.  Figure
46: Phase diagrams of R7G-hBN superlattices of various stacking
configurations: (a) $\xi=1$ , (b) $\xi=0$, (c) $(\xi_{b},\xi_{t})=(1,1)$, and
(d) $(\xi_{b},\xi_{t})=(1,0)$. Panel (i) shows the single-particle indirect
gap $\Delta_{\text{ind},\nu=-4}$ at filling $\nu=-4$; (ii) shows the indirect
gap $\Delta_{\text{ind},\nu=0}$ at filling $\nu=0$; (iii) shows the indirect
gap $\Delta_{\text{ind},\nu=4}$ at filling $\nu=4$; (iv) shows the minimal
direct gap $\Delta_{\text{val}}$ around the highest valence band in one
valley; (v) shows the bandwidth $W_{\text{val}}$ of the highest valence band;
(vi) shows the minimal direct gap $\Delta_{\text{cond}}$ around the lowest
conduction band in one valley; (vii) shows the bandwidth $W_{\text{cond}}$ of
the lowest conduction band. Chern numbers of the highest valence band and the
lowest conduction band in K valley are indicated on panel (iv) and panel (vi),
respectively, where boundaries of topologically distinct phases can be seen as
the direct gap closes.  Figure 47:  The comparison between the DFT+SK bands
(black) and the bands from the $2\times 2$ effective continuum model (red
lines) in Eq. (86) for $n=3$. The parameter values are listed in Tab. 2 and
Tab. 4 of Main Text.  Figure 48:  The comparison between the DFT+SK bands
(black) and the bands from the $2\times 2$ effective continuum model (red
lines) in Eq. (86) for $n=4$. The parameter values are listed in Tab. 2 and
Tab. 4 of Main Text.  Figure 49:  The comparison between the DFT+SK bands
(black) and the bands from the $2\times 2$ effective continuum model (red
lines) in Eq. (86) for $n=5$. The parameter values are listed in Tab. 2 and
Tab. 4 of Main Text.  Figure 50:  The comparison between the DFT+SK bands
(black) and the bands from the $2\times 2$ effective continuum model (red
lines) in Eq. (86) for $n=6$. The parameter values are listed in Tab. 2 and
Tab. 4 of Main Text.  Figure 51:  The comparison between the DFT+SK bands
(black) and the bands from the $2\times 2$ effective continuum model (red
lines) in Eq. (86) for $n=7$. The parameter values are listed in Tab. 2 and
Tab. 4 of Main Text.