# Topological flat bands in rhombohedral tetralayer and multilayer graphene on
hexagonal boron nitride moire superlattices

Youngju Park  Yeonju Kim  Department of Physics, University of Seoul, Seoul,
02504, Korea  Bheema Lingam Chittari  Department of Physical Sciences, Indian
Institute of Science Education and Research Kolkata, Mohanpur 741246, West
Bengal, India  Jeil Jung<EMAIL_ADDRESS>Department of Physics, University
of Seoul, Seoul, 02504, Korea  Department of Smart Cities, University of
Seoul, Seoul 02504, Korea

###### Abstract

We show that rhombohedral four-layer graphene (4LG) nearly aligned with a
hexagonal boron nitride (hBN) substrate often develops nearly flat isolated
low energy bands with non-zero valley Chern numbers. The bandwidths of the
isolated flatbands are controllable through an electric field and twist angle,
becoming as narrow as $\sim$10 meV for interlayer potential differences
between top and bottom layers of $|\Delta|\approx 10\sim 15$ meV and
$\theta\sim 0.5^{\circ}$ at the graphene and boron nitride interface. The
local density of states (LDOS) analysis shows that the nearly flat band states
are associated to the non-dimer low energy sublattice sites at the top or
bottom graphene layers and their degree of localization in the moire
superlattice is strongly gate tunable, exhibiting at times large
delocalization despite of the narrow bandwidth. We verified that the first
valence bands’ valley Chern numbers are $C^{\nu=\pm 1}_{V1}=\pm n$,
proportional to layer number for $n$LG/BN systems up to $n=8$ rhombohedral
multilayers.

##  I Introduction

Nearly aligned van der Waals 2D layered heterostructures leads to moire
superlattices Hass et al. (2008); Miller et al. (2009, 2010); Sadowski et al.
(2006); De Heer et al. (2010); Brihuega et al. (2012); Ohta et al. (2012);
Lopes dos Santos et al. (2012, 2007); Shallcross et al. (2010, 2008); Landgraf
et al. (2013); Shallcross et al. (2013); Bistritzer and MacDonald (2011); Moon
and Koshino (2013, 2012); Jung et al. (2014); San-Jose et al. (2012); San-Jose
and Prada (2013); Stauber et al. (2013); Bistritzer and MacDonald (2010); Wang
et al. (2012); Schmidt et al. (2014); Carr et al. (2018); Koshino et al.
(2018); Kang and Vafek (2018, 2019); Tarnopolsky et al. (2019); Po et al.
(2019) with enlarged moire lattice constants ($\ell_{m}$) and reduced moire
Brillouin zones (MBZ) where the bandwidth suppression typically enhances the
Coulomb correlation effects. Twisted bilayer graphene (tBG) is a
representative system showing a variety of ordered phases near its magic angle
($\approx 1^{\circ}$) as a function of carrier doping Cao et al. (2018a); Kim
et al. (2017); Sharpe et al. (2019); Cao et al. (2018b); Yankowitz et al.
(2019); Cao et al. (2020). Other related graphene moire systems such as
twisted trilayer graphene (tTG) Lopez-Bezanilla and Lado (2020); Shin et al.
(2021); Park et al. (2021); Hao et al. (2021); Ramires and Lado (2021); Qin
and MacDonald (2021); Phong et al. (2021); Kim et al. (2022); Fischer et al.
(2022); Shen et al. (2023); Yankowitz (2023), twisted monolayer-bilayer
graphene (tMBG) Park et al. (2020); Rademaker et al. (2020); Ma et al. (2021);
Xu et al. (2021); Chen et al. (2021); He et al. (2021); Li et al. (2022); Tong
et al. (2022), and twisted double bilayer graphene (tDBG) Chebrolu et al.
(2019); Liu et al. (2020); Lee et al. (2019); Shen et al. (2020); Choi and
Choi (2019); Koshino (2019); Li et al. (2020); Sinha et al. (2020) are current
systems of interest. The gapped massive Dirac layer-based heterostructures
Song et al. (2013); Javvaji et al. (2020) including transition metal
dichalcogenides homo and hetero twisted bilayers Naik and Jain (2018); Naik et
al. (2020); Wu et al. (2019); Zhan et al. (2020); Wang et al. (2020); Zhang et
al. (2020); Shabani et al. (2021); Rademaker (2022) achieves reduced
bandwidths at relatively large twist angles allowing to enhance the Coulomb
interactions. Another interesting 2D layered moire system is formed by
graphene on hexagonal boron nitride (G/BN) moire superlattices Bistritzer and
MacDonald (2011); Lopes dos Santos et al. (2012); Jung et al. (2014); Wong et
al. (2015); Kerelsky et al. (2019); Choi et al. (2019); Moon and Koshino
(2014) with lattice mismatch Song et al. (2013); Yankowitz et al. (2012a).
When a graphene layer ($a_{\mathrm{G}}=2.46~{}\AA$) is aligned to the
hexagonal boron nitride ($a_{\mathrm{hBN}}=2.5025~{}\AA$), the lattice
mismatch ($\epsilon=a_{\mathrm{G}}/a_{\mathrm{hBN}}-1\approx 1.7\%$) generates
a moire pattern leading to secondary Dirac cone features Yankowitz et al.
(2012b); Dean et al. (2013); Ponomarenko et al. (2013) and a primary Dirac
point band gap formsAmet et al. (2013); Hunt et al. (2013) due to the moire
pattern strains Woods et al. (2014); Jung et al. (2015); San-Jose et al.
(2014). Those moire strain profiles are tunable with twist angle
$\theta=0^{\circ}\sim 1^{\circ}$ leading to superlattice lengths
${\ell_{m}}\approx{a_{G}}/{\sqrt{\epsilon^{2}+\theta^{2}}}\approx$ 14$\sim$10
nm that modifies band properties of the aligned single layer graphene and
multilayer graphene Leconte et al. (2017); Kim et al. (2018); Leconte and Jung
(2020); Chen et al. (2019a, b); Chittari et al. (2019); Zhang et al. (2019);
Chen et al. (2020); González et al. (2021); Yang et al. (2022); Calderón et
al. (2022). The aligned ABC stacked rhombohedral trilayer graphene on
hexagonal boron nitride (3LG/BN) Chen et al. (2019a, b); Chittari et al.
(2019); Zhang et al. (2019); Chen et al. (2020); González et al. (2021); Yang
et al. (2022); Calderón et al. (2022) is one of the important systems that
develop interface-interaction induced band isolation where a perpendicular
electric field can achieve narrow bandwidths comparable to the Coulomb energy
($U\approx 25~{}$meV) Chen et al. (2019a, b, 2020); Yang et al. (2022), making
it a powerful platform for studying flat band driven phenomena. In this
manuscript we demonstrate that in rhombohedral four-layer graphene on boron
nitride (4LG/BN) and $n$LG/BN systems up to $n=8$ we can achieve even narrower
nearly flat bands than those of 3LG/BN whose bandwidths are typically
$W\approx 30\sim 40~{}$meV in the absence of an electric field. While the low
energy bandwidths become progressively narrower for increasing number of
layers $n$ in rhombohedral multilayers, with 4LG narrowing between a factor
two to four over 3LG, comparison against 5LG, 6LG, 7LG, and 8LG systems
aligned on hexagonal boron nitride in the absence of perpendicular electric
fields reveals that 5LG systems can already host optimally narrowest bands
comparable to those of 6$\sim$8LG depending on the specific moire substrate
potential used. The manuscript is structured as follows. In Sec. II we present
the full-bands continuum model Hamiltonian, In Sec. III we discuss the results
on the bandwidths, Chern numbers and local density of states as a function of
twist angles and electric fields, and in Sec. IV we present the summary and
conclusions. Extended data for $n=5\sim 8$ systems are presented in the
appendix.

##  II Model Hamiltonian

We model the rhombohedral tetralayer graphene on hexagonal boron nitride
(4LG/BN) as

$\displaystyle{H}={H}_{4LG}+{H}_{M}$ (1)

as a sum of the four layers Hamiltonian $H_{4LG}$ and the effective intralayer
moire pattern potential $H_{M}$ in the graphene layer contacting BN. The
$H_{4LG}$ term is a full-bands tight binding bands model

$\displaystyle H_{4LG}=\begin{pmatrix}H_{11}&H_{12}&H_{13}&{\bm{0}}_{2\times
2}\\\ H_{12}^{\dagger}&H_{22}&H_{12}&H_{13}\\\
H_{13}^{\dagger}&H_{12}^{\dagger}&H_{33}&H_{12}\\\ {\bm{0}}_{2\times
2}&H_{13}^{\dagger}&H_{12}^{\dagger}&H_{44}\end{pmatrix},$ (2)

whose hopping parameters are chosen to match the LDA density functional theory
where the intra-$(H_{ll})_{2\times 2},(l=1,2,3,4)$ and inter-layer
$(H_{ij})_{2\times 2},(i\neq j=1,2,3)$ hamiltonian terms are given by

$\displaystyle H_{ll}({\bm{k}})$
$\displaystyle=\begin{pmatrix}u_{A_{l}}&\upsilon_{0}\pi^{\dagger}\\\
\upsilon_{0}\pi&u_{B_{l}}\end{pmatrix}+V_{ll}\mathbb{1},$ (3) $\displaystyle
H_{12}({\bm{k}})$
$\displaystyle=\begin{pmatrix}-\upsilon_{4}\pi^{\dagger}&-\upsilon_{3}\pi\\\
t_{1}&-\upsilon_{4}\pi^{\dagger}\end{pmatrix},$ $\displaystyle
H_{13}({\bm{k}})$ $\displaystyle=\begin{pmatrix}0&t_{2}\\\ 0&0\end{pmatrix}.$

where $\pi=(\nu p_{x}+ip_{y})$ is defined in terms of the valley index
$\nu=\pm 1$ using the momentum vector ${\bm{p}}=(p_{x},p_{y})$ measured from
the principal Dirac points
${\bm{K}}_{\nu}=\left(\nu\frac{4\pi}{3a_{G}},0\right)$. The intralayer nearest
neighbor sublattices ($A_{l}~{}\&~{}B_{l}$ where $l=1,2,3,4$ is the layer
index) are connected through the $t_{0}=-3.1$ eV hopping term whose magnitude
is slightly larger than $t_{0}\simeq-2.6~{}$eV obtained from density
functional theory (DFT) local density approximation (LDA) to partially account
for the Coulomb interaction-driven Fermi velocity enhancement. We use the
Fermi velocity parameters $\upsilon_{i}=(\sqrt{3}a/2\hbar)|t_{i}|$ in the
Hamiltonian. The term $V_{ll}\mathbb{1}$ is used to introduce the interlayer
potential difference ($\Delta$) between contiguous layers through a
perpendicular external electric field. For convenience we use equal magnitude
potential drops proportional to $\Delta$ given by

$V=\Delta\left(\frac{3}{2},\frac{1}{2},-\frac{1}{2},-\frac{3}{2}\right)$ (4)

such that the interlayer potential difference between top and bottom graphene
layers is given by $(n-1)\Delta$ where $n=4$ for a tetralayer. The remote
hopping term between the interlayer sublattices $A_{l}~{}\&~{}B_{l+1}$ is
$t_{3}=0.293$ eV, and between $A_{l}(B_{l})~{}\&~{}A_{l+1}(B_{l+1})$ is given
by $t_{4}=0.144$ eV, see Fig. 2. The hopping energy between the adjacent
interlayer vertical sublattices ($B_{l}~{}\&~{}A_{l+1}$) is $t_{1}=0.3561$ eV,
and for the sublattices ($A_{l}~{}\&~{}B_{l+2}$) the hopping energy is
$t_{2}=-0.0083$ eV. The diagonal site potentials $u_{A_{l}}$ ($u_{B_{l}}$) at
each sublattice of 4LG are

$\displaystyle u_{A1}$ $\displaystyle=u_{B4}=0~{}\textrm{eV},$ (5)
$\displaystyle u_{B1}$ $\displaystyle=u_{A4}=0.0122~{}\textrm{eV},$
$\displaystyle u_{A2/A3}$ $\displaystyle=u_{B2/B3}=-0.0164~{}\textrm{eV}.$

In Fig. A1 we confirm the close agreement of the DFT-LDA bands with the full-
bands tight binding model from 3LG up to 8LG especially in the low energy
range of $\pm 0.05~{}$eV and near the Dirac point
$|{\bm{k}}|<0.05\left({4\pi}/{3a_{G}}\right)$.

Figure 1: (color online) The atomic structure of four-layer graphene (4LG) on
boron nitride moire superlattices with two possible hexagonal boron nitride
(hBN) orientations is shown in (a) $\xi=1$ (4LG/BN) and (b) $\xi=-1$ (4LG/NB).
The moire superlattice formation due to lattice mismatch is illustrated for
the graphene layers (1LG, 2LG, 3LG, and 4LG) where the local commensurate
stacking are illustrated with AA, AB, BA labels. We show the local atomic
structures illustrating the rhombohedral stacking sublattices Al and Bl, where
$l=1,2,3,4$ is the layer index, together with the intralayer ($t_{0}$) and
interlayer ($t_{1},~{}t_{2},~{}t_{3},~{}t_{4}$) nearest neighbor hopping
terms. (c) We illustrate the moire Brillouin zone (MBZ) without (upper) and
with (lower) a twist angle next to the Brillouin zones (BZ) of graphene (red
hexagon) and hBN (blue hexagon) layers. The shaded hexagons indicate the MBZ
represented at the two valleys $\nu=\pm 1$ of graphene. We define the Dirac
points of the pristine 4LG as the center of the MBZ ($\tilde{\Gamma}$), and
one of the corners of MBZ ($\tilde{K}$) corresponds to the Dirac points of hBN
for each valley $\nu=\pm 1$.  Figure 2: (color online) We present the band
structure plots of 4LG/BN for two hBN orientations ($\xi=\pm 1$) at twist
angles of $\theta=0^{\circ}$ and $\theta=0.54^{\circ}$, with $\Delta=0$ eV and
for valley $\nu=1$ ($\nu=-1$) shown using solid (dashed) lines. To facilitate
analysis, we use a compact representation of the local density of states
(LDOS), $D({\bm{r}},E)$, along the local commensurate stacking (AA, AB, and
BA). We observe that the Van Hove singularities (vHS) are either localized to
AB (AA) stacking for $\xi=1$ ($\xi=-1$), or spread within the moiré unit cell.
Delocalization occurs at the isolated valence bands, where the normalized
density of states (DOS), $\tilde{D}(E)=D(E)/{\rm max}(D(E))$, clearly shows
their band isolation. Additionally, we calculate the Berry curvatures
($\Omega$) of the valence (red) and conduction (blue) bands for valley
$\nu=+1$. The unequal positive and negative weights of Berry curvatures near
the MBZ corners result in non-zero valley Chern numbers, $C_{C1}^{+}=1$ and
$C_{V1}^{+}=4$.

We capture the effect of G/BN moire superlattice by adding the effective
intralayer moire potential $H^{\xi}_{M}$ acting at the bottom layer of 4LG

$\displaystyle H^{\xi}_{M}({\bm{r}})=$ $\displaystyle
V^{\xi}_{AA}({\bm{r}})\left(\frac{\mathbb{1}+\xi\sigma_{z}}{2}\right)+V^{\xi}_{BB}({\bm{r}})\left(\frac{\mathbb{1}-\xi{\sigma}_{z}}{2}\right)$
(6)
$\displaystyle+V_{BA}^{\xi}({\bm{r}})\cdot{\sigma}^{\xi}_{xy}\,\delta_{\nu,1}+V_{AB}^{\xi}({\bm{r}})\cdot{\sigma}^{\xi}_{xy}\,\delta_{\nu,-1}.$

Here, ${\bm{\sigma}}^{\xi}_{xy}=\left(\sigma_{x},\xi\sigma_{y}\right)$ and
$\sigma_{z}$ are the Pauli matrices in the sublattice basis. The $0^{\circ}$
and $60^{\circ}$ orientations of the hBN substrate give rise to different
moire interlayer potentials. We use the $\xi=1$ label for the $0^{\circ}$
orientation of hBN sheet that for $AA$ stacking the B and N atoms are right
below $A_{1}$ and $B_{1}$ sites respectively, while the AA stacking for
$\xi=-1$ corresponding to $60^{\circ}$ orientation and the N and B atoms are
below $A_{1}$ and $B_{1}$ sites, see Fig. 1(a) $\&$ (b).

The modifications of the onsite potential and inter sublattice interaction in
the bottom layer graphene of 4LG contacting hBN to capture the effective
intralayer moire potential is given by the following equations:

$V^{\xi}_{AA(BB)}({\bm{r}})=2\,C_{AA(BB)}\mathrm{Re}\left[e^{i\phi_{AA(BB)}}\,f^{\xi}({\bm{r}})\right],$
(7) $\displaystyle V^{\xi}_{BA}({\bm{r}})$
$\displaystyle=\left(V^{\xi}_{AB}({\bm{r}})\right)^{*}=2\,C_{AB}\left(\mathrm{cos}(\phi)\hat{z}\times\mathbb{1}-\mathrm{sin}(\phi)\mathbb{1}\right)$
(8)
$\displaystyle\times\left(\frac{1}{|{\tilde{\bm{G}}_{1}}|}\nabla\mathrm{Re}\left[e^{i(-\phi_{AB}+\pi/6)}\,f^{\xi}({\bm{r}})\right]\right),$

where the coefficients adopted from Ref. Jung et al. (2014) are
$C_{AA}=-14.88$ meV, $\phi_{AA}=50.19^{\circ}$, $C_{BB}=12.09$ meV,
$\phi_{BB}=-46.64^{\circ}$, $C_{AB}=11.34~{}$ meV, $\phi_{AB}=19.60^{\circ}$,
and the function as
$f^{\xi}({\bm{r}})=\sum_{m=1}^{6}{\left(\frac{1+(-1)^{m+1}}{2}\right)e^{-i\xi{\tilde{\bm{G}}_{m}}\cdot{\bm{r}}}}$.
The moire reciprocal lattice vectors are defined as
${\tilde{\bm{G}}_{m}}=\hat{R}_{\frac{\pi}{6}(m-1)}{\tilde{\bm{G}}_{1}}$ where
${\tilde{\bm{G}}_{1}}=(\frac{a_{\mathrm{G}}}{a_{\mathrm{hBN}}}\hat{R}_{\theta}-\mathbb{1})^{\dagger}[0,\frac{4\pi}{\sqrt{3}a_{\mathrm{G}}}]\approx(\epsilon\mathbb{1}-\theta\hat{z}\times\mathbb{1})[0,\frac{4\pi}{\sqrt{3}a_{{\bm{G}}}}]$
for the lattice mismatch $\epsilon\approx-1.7\%$ and twist angle $\theta$
between the graphene and hBN Jung et al. (2015). In Fig. 1(c), we illustrated
the moire Brillouin zone (MBZ) with the original Brillouin zone of the 4LG and
hBN layers to show the effect of the lattice mismatch and twist angle. The
reduced size of MBZ leads to a repulsion gap at the MBZ zone corners in the
4LG/BN and isolates the bands near the Fermi level.

Figure 3: (color online) The low energy bandwidth variation for the parameter
set of interlayer potential difference ($\Delta$) and twist angle ($\theta$)
is presented for the two orientations of hBN, (a) $\xi=1$ and (c) $\xi=-1$. In
(b) and (d), we present the corresponding primary ($\delta_{p}$) and secondary
gaps for valence ($\delta_{s}~{}(V1)$) and conduction ($\delta_{s}~{}(C1)$)
bands, to estimate the isolation of valence and conduction bands at CNP.

##  III Results and Discussion

In this section we present the valley Chern numbers, and local density of
states that are associated with the large $U_{\rm eff}/W$ values where strong
correlations can be expected in the parameter space of twist angles and
interlayer potential differences.

The low energy band structure of 4LG/BN have extremely narrow bands near the
charge neutrality point (CNP) thanks to the $\sim p^{N}$ behavior of the band
dispersion near the Dirac point of the minimal model rhombohedral multilayer
graphene as the layer number $N$ increases Min and MacDonald (2008). More
realistic band structures include distortions due to trigonal warping and
electron-hole symmetry breaking terms included in the remote hopping terms of
the Hamiltonian as illustrated in the LDA bands from 3LG to 8LG in Fig. A1 in
the appendix. Addition of a moire pattern introduces a band gap both at the
charge neutrality and near the MBZ corners that isolates the low energy
valence and conduction bands from the adjacent bands even in the absence of an
external electric field, see Fig. 2. The band structures are shown for 4LG/BN
for two different orientations of hBN, namely 4LG/BN ($\xi=1$) and 4LG/NB
($\xi=-1$), and we considered the aligned $\theta=0^{\circ}$ and a small
finite twist angle of $\theta=0.54^{\circ}$. The low energy bands have
bandwidths between $W\approx 12\sim 20~{}$meV where the bandwidth of the first
valence ($V1$) or conduction ($C1$) band is calculated as
$W_{V1(C1)}=E^{V1(C1)}_{\max}-E^{V1(C1)}_{\min}$. We note that the low energy
valence and conduction bandwidths of 4LG/BN and all the way up to 8LG/BN are
smaller than those of 3LG/BN Chittari et al. (2019) (See Fig. A2(c) $\&$ (d)).

The real-space distribution of the wave functions associated with the nearly
flat bands in 4LG/BN gives an idea about the kind of broken symmetries that
the system can host. We calculated the density of states (DOS) $D(E)$ through

$\displaystyle D(E)$ $\displaystyle=\int_{\rm
MBZ}\frac{d{\bm{k}}}{(2\pi)^{2}}\,\sum_{n^{\prime}}|\psi_{n^{\prime}}({\bm{k}})|^{2}\,\delta(E-E_{n^{\prime}}({\bm{k}}))$
(9)

and the local density of states (LDOS) $D({\bm{r}},E)$ through

$\displaystyle D({\bm{r}},E)$ $\displaystyle=\int_{\rm
MBZ}\frac{d{\bm{k}}}{(2\pi)^{2}}\,\sum_{n^{\prime}}|\psi_{n^{\prime}}({\bm{k}}:{\bm{r}})|^{2}\,\delta(E-E_{n^{\prime}}({\bm{k}}))$
(10)

using information from the real-space wave functions
$\psi_{n^{\prime}}({\bm{k}}:{\bm{r}})=\sum_{\tilde{G}}\,\psi_{n^{\prime}}({{\bm{k}}+\tilde{\bm{G}}})\,e^{-i({\bm{k}+\tilde{\bm{G}}})\cdot{\bm{r}}}$
where the band indices $n^{\prime}=Vn(Cn)$ refer to the $n^{\rm th}$ valence
(conduction) bands counting from the charge neutrality point. In Fig. 2, we
show the compact representation of the LDOS, $D({\bm{r}},E)$ at each local
commensurate stacking (AA, AB and BA), along with the normalized density of
states defined as $\tilde{D}(E)=D(E)/{\rm max}(D(E))$. The density of states
$\tilde{D}(E)$ plots show the van Hove singularities (vHS) of DOS peaks for
the isolated low energy valence and conduction bands. Additionally, the LDOS
plots show that the vHS of the conduction band is strongly localized at AB(AA)
stacking for $\xi=1$ ($\xi=-1$) at both $\theta=0^{\circ},0.54^{\circ}$.
However, the valence band states at the vHS distribute the carrier densities
almost equally at all local commensurate stacking sites. These differences in
the localization behavior between the conduction and valence band states
should in turn lead to very different spin or charge textures of the
associated ordered phases upon inclusion of Coulomb interactions.

Furthermore, the low energy isolated bands show a non-trivial topological
nature. We calculated the valley resolved Chern number ($C^{\nu}_{V1(C1)}$) of
the first valence (conduction) band using the following equations:

$\displaystyle C_{n^{\prime}}$ $\displaystyle=\int_{\rm
MBZ}d^{2}{{\bm{k}}}\quad\Omega_{n^{\prime}}({\bm{k}})/(2\pi),$ (11)

where the Berry curvature $\Omega_{n^{\prime}}({\bm{k}})$ D. et al. (2010) is
defined as

$\displaystyle\Omega_{n^{\prime}}(\bm{{k}})$ $\displaystyle=-2\sum_{i\neq
n^{\prime}}\mathrm{Im}\left[\frac{\matrixelement{u_{n^{\prime}}}{\frac{\partial
H}{\partial{k}_{x}}}{u_{i}}\matrixelement{u_{i}}{\frac{\partial
H}{\partial{k}_{y}}}{u_{n^{\prime}}}}{(E_{i}-E_{n^{\prime}})^{2}}\right],$
(12)

for the $n^{\prime}=Vn(Cn)$ valence (conduction) bands from the charge
neutrality point where $\ket{u_{n^{\prime}}}$ are the moire Bloch states of
index $n^{\prime}$, and $E_{n^{\prime}}$ are the band energies. The calculated
Berry curvatures for the low energy bands are represented in Fig. 2 for
$\xi=\pm 1$ and $\theta=0^{\circ},0.54^{\circ}$. The Berry curvature plots
have hot spots near the MBZ corners, and the unequal weights of Berry
curvatures give rise to non-zero valley Chern numbers $C^{\nu=\pm
1}_{V1}=\pm~{}4$ for the low energy valence bands of 4LG/BN at twist angles
$\theta=0^{\circ},0.54^{\circ}$ with $\xi=\pm 1$, and we get $C^{\nu=\pm
1}_{C1}=\pm~{}1$ for the low energy conduction bands for $\xi=1$ but zero for
$\xi=-1$. We will discuss later how the valley Chern number phase diagram can
change as a function of interlayer potential difference.

###  III.1 Low energy bands

In the following, we discuss the bandwidth of the low energy bands near the
Fermi level of 4LG aligned to hBN in the parameter space of twist angles
($\theta$) and interlayer potential differences ($\Delta$). In our previous
reports Javvaji et al. (2020); Chebrolu et al. (2019); Park et al. (2020);
Chittari et al. (2019), we systematically demonstrated that the interlayer
potential difference introduced through a perpendicularly applied electric
field is an effective control knob to tune the bandwidth and band isolation of
the low energy bands. The interlayer potential difference ($\Delta$) results
from unequal intralayer potentials $V_{ii}\mathbb{1}$ in each layer of 4LG/BN
in Eq. (3). In Fig. 3(a) and (c), we show the variation of bandwidth as a
function of $\Delta$ and $\theta$. The on-site Coulomb repulsion energy for
graphene on hBN is estimated to be $U\approx 25~{}$meV for the moire length of
$\sim 14.5~{}$nm. From the bandwidth phase diagram, it is evident that a large
set of the parameter space of $\Delta$ and $\theta$ shows a bandwidth smaller
than $\sim$25 meV. Further, we also investigate the isolation of the low
energy flat bands from the adjacent remote bands by calculating the gaps at
CNP ($\delta_{p}$) and secondary gaps ($\delta_{s}$) as a function of $\Delta$
and $\theta$ (see Fig. 3). Here, the primary gap at CNP is defined as
$\delta_{p}=E^{C1}_{\min}-E^{V1}_{\max}$ and the secondary gap (a gap between
the valence (conduction) band from its higher energy bands) is defined as
$\delta_{s}(n^{\prime})=E^{n^{\prime}}_{\min}-E^{n^{\prime}}_{\max}$ with
$n^{\prime}=V1(C1)$. A positive value denotes the magnitude of the gap, and a
negative value indicates the degree of band overlap. The low energy bands
isolate when both secondary and primary gaps are simultaneously positive. The
primary gap at CNP is always open practically for all values of $\theta$ and
$\Delta$ except for a specific value of $\Delta$, see Fig. 3.

A few important observations from the bandwidth and bandgap phase diagrams are
summarized below. (i) Very narrow valence (conduction) bands ($W<$ 5 meV) are
possible for positive(negative) values of $\Delta\approx 0.012\sim
0.013~{}{\rm eV}$ ($\Delta\approx-0.025\sim-0.038~{}{\rm eV}$) at zero-twist
($\theta=0^{\circ}$) for both $\xi=\pm 1$. (ii) The small twist angles
$\theta=0\sim 0.8^{\circ}$ retain the smallest bandwidths ($W<$ 5 meV) for the
valence and conduction bands. (iii) For $\xi=1$, the conduction bands are
isolated for the specific conditions of $\theta\approx 0.5\sim 0.6^{\circ}$
and $\Delta=0.007\sim 0.027$ eV while the valence bands are isolated for the
considered range of $\theta=0\sim 0.8^{\circ}$ for $\Delta=-0.05$ to
$0.01~{}{\rm eV}$. (iv) For $\xi=-1$, the conduction bands are isolated for
the range of $\Delta=$0 to 0.04 eV and $\theta=0^{\circ}\sim 0.5^{\circ}$
while the valence bands are isolated when $\theta\approx 0.5\sim 0.6^{\circ}$
and $\Delta=-0.004$ to $-0.02$ eV. From the bandwidth and band gap phase
diagrams we can observe that the narrowest bandwidths of $W=$ 5$\sim$10 meV
for the isolated valence (conduction) bands is associated with
$\Delta=-0.009~{}(0.009)$ eV for $\xi=1$ and $\Delta=-0.011~{}(0.015)$ eV for
$\xi=-1$ for the twist angle range $\theta=0^{\circ}\sim 0.5^{\circ}$. See
Fig. A3 where we summarize band structures for these selected system
parameters. These observations help us draw a few important conclusions on the
possible Coulomb interation driven phases of 4LG/BN that we will discuss
further in the following subsection.

###  III.2 Nearly flat bands ($U_{\rm eff}/W\gg 1$ )

Figure 4: (color online) We present a phase diagram showing the ratio between
the bare on-site Coulomb repulsion energy ($U$) and low energy band bandwidths
($W$) as a function of an interlayer potential difference ($\Delta$) and twist
angles ($\theta$) for two orientations of hBN: (a) $\xi=1$ and (b) $\xi=-1$.
The colored region represents the narrow bands where $U/W>1$, indicating
strong electronic correlations. We find that almost all the system parameters
of twist angle $\theta=0\sim 0.8^{\circ}$ and interlayer potential difference
$|\Delta|<0.1~{}$eV exhibit narrow bands satisfying $U/W>1$ for both $\xi=\pm
1$ hBN orientation, suggesting that strong electronic correlations are present
in a wide range of experimental conditions.

Here we compare the on-site Coulomb repulsion energy ($U$) and bandwidth ($W$)
of valence and conduction bands in search for the narrow bandwidth regime with
strong effective Coulomb interactions where $U_{\rm eff}/W\gg 1$. The greater
the $U$ compared to the bandwidth $W$ ($U/W>1$) indicates higher chances of
finding Coulomb interaction induced correlated phases. We define the on-site
Coulomb repulsion energy $U$ using the moire superlattice constant
${\ell_{m}}\approx a_{G}/\sqrt{\epsilon^{2}+\theta^{2}}$ and a relative
permittivity $\epsilon_{r}=4$

$\displaystyle U$
$\displaystyle=\frac{e^{2}}{4\pi\epsilon_{r}\epsilon_{0}{\ell_{m}}},$ (13)

where $e$ is the electron charge, and $\epsilon_{0}$ is the permittivity of
vacuum. In Fig. 4, we summarized the ratio ($U/W$) of the low energy bandwidth
$W$ versus the on-site Coulomb repulsion energy $U$ as a function of $\theta$
and $\Delta$, where the colored area manifests the region for the possible
narrow bands ($U/W>1$). Compared to the estimated Coulomb repulsion
energy($\approx 25$ meV), the bandwidths of valence and conduction bands in
4LG/BN are small enough to satisfy the $U/W>1$ condition for most of the
considered parameter space.

Figure 5: (color online) We present a phase diagram illustrating the variation
of the ratio between effective on-site Coulomb repulsion energy ($U_{\rm
eff}$) and bandwidths ($W$) as a function of interlayer potential difference
($\Delta$) and twist angles ($\theta$) for two orientations of hBN: (a)
$\xi=1$ and (b) $\xi=-1$. On the right side of the diagram, we show $U_{\rm
eff}/W$ results at two specific twist angles: $\theta=0^{\circ}$ (cyan) and
$\theta=0.54^{\circ}$ (magenta). For $\theta=0^{\circ}$, $U_{\rm eff}/W>1$ is
only achievable for the valence band with negative $\Delta$ for $\xi=1$, and
for the conduction band with positive $\Delta$ for $\xi=-1$. However, at
$\theta=0.54^{\circ}$, both valence and conduction bands exhibit $U_{\rm
eff}/W\gg 1$, with negative $\Delta$ required for the valence band and
positive $\Delta$ for the conduction band. Our results reveal that $U_{\rm
eff}/W>1$ for both conduction and valence bands can be achieved at
$\Delta=0.009~{}$eV for $\xi=1$ and $\theta=0.54^{\circ}$.

We further consider the screening effects by adding an exponential decay term
$e^{-{\ell_{m}}/\lambda_{D}}$ Chebrolu et al. (2019) adequately when there is
band overlap. Here, the screened effective on-site Coulomb repulsion energy
($U_{\rm eff}$) is redefined as:

$\displaystyle U_{\rm{eff}}$
$\displaystyle=U~{}\textrm{exp}(-{\ell_{m}}/\lambda_{D}),$ (14)

where $\lambda_{D}=2\epsilon_{0}/e^{2}D(\delta_{p},\delta_{s})$ is the Debye
length, which includes the 2D density of states
$D(\delta_{p},\delta_{s})=4[|\delta_{p}|\,u(-\delta_{p})\,+\,|\delta_{s}|\,u(-\delta_{s})]/(W^{2}A_{M})$.
$A_{M}=\sqrt{3}{\ell_{m}}^{2}/2$ is the moire unit cell area, $W$ is the
bandwidth, and $u(x)$ is the Heaviside step function such that
$u(-\delta_{p(s)})$) enhances the screening in the presence of band overlap
($\delta_{p(s)}<0$). We show that the screening effect due to the band overlap
reduces the narrow bandwidth regime ($U_{\rm eff}/W>1$), see Fig. 5 where we
summarized $U_{\rm eff}/W$ for the parameter space of interlayer potential
difference $\Delta$ and twist angle $\theta$, and specifically for
$\theta=0^{\circ}$ (cyan) and $~{}0.54^{\circ}$ (magenta). For $\xi=1$ (Fig. 5
(a)), the valence band show $U_{\rm eff}/W\gg 1$ for a large span of the
parameter spaces ($\Delta\approx-50\sim 10$ meV and $\theta=0\sim
0.8^{\circ}$) while for the conduction band $U_{\rm eff}/W>1$ is only possible
for the small island at $\theta\approx 0.5\sim 0.6^{\circ}$ and $\Delta\approx
20\sim 25$ meV. In another substrate configuration, $\xi=-1$ (Fig. 5 (b)), the
conduction band has greater chances of achieving the flat bands regime at
small positive $\Delta=0\sim 40$ meV compared to the valence band, which has a
limited area around $\Delta\approx 20$ meV and $\theta=0.4\sim 0.6^{\circ}$.

Moreover, as shown in the phase diagram, the 4LG/BN flat bands have different
particle-hole asymmetric behavior for different $\theta$. We found that at
$\theta=0^{\circ}$ the valence (conduction) band only satisfies $U_{\rm
eff}/W\gg 1$ for $\xi=1$($\xi=-1$) substrate orientation. On the contrary, at
$\theta=0.54^{\circ}$, both valence and conduction flat bands are possible,
which requires negative $\Delta$ for valence bands and positive $\Delta$ for
the conduction bands. Also, the maximum values of $U_{\rm eff}/W$ are observed
when $\theta=0.54^{\circ}$ giving $U_{\rm eff}/W=3.98$ for the valence band
with $\xi=1$ orientation and $U_{\rm eff}/W=4.7$ for the conduction band with
$\xi=-1$. Hence, we suggest that $\theta=0.54^{\circ}$ for the 4LG/BN device
has greater chances of achieving Coulomb-interaction driven ordered states in
a wider range of interlayer potential differences.

Figure 6: (color online) The valley Chern number of the low-energy valence and
conduction bands is plotted as a function of interlayer potential difference
($\Delta$) and twist angle ($\theta$) for two substrate orientations: (a)
$\xi=1$ and (b) $\xi=-1$. The dashed lines indicate two specific twist angles,
$\theta=0^{\circ}$ (navy) and $\theta=0.54^{\circ}$ (magenta), and the phase
transition between different Chern numbers is shown on the right-hand side of
the phase diagram. The red (valence) and blue (conduction) dashed lines
enclose regions where the effective on-site Coulomb repulsion energy $U_{\rm
eff}$ is greater than the band widths $W$, i.e., $U_{\rm eff}/W>1$.

###  III.3 Valley Chern numbers

When the moire bands are isolated through gaps, they can acquire finite valley
resolved Chern numbers Chittari et al. (2019); Zhang et al. (2019). The lowest
energy bands of 4LG/BN in the absence of an external electric field have
finite valley Chern numbers for both valence and conduction bands when
$\xi=1$, and for valence bands only for $\xi=-1$, or in other words the
valence band has a finite valley Chern number $C_{V1}=4$ for both $\xi=\pm 1$
substrate orientations, while the conduction band has the valley Chern number
$C_{C1}=1$ only for $\xi=1$.

Here we analyze the phase diagram of the valley Chern numbers in the parameter
space of interlayer potential differences $\Delta$ and twist angles $\theta$
for the two different types of hBN substrate orientations labeled through
$\xi=\pm 1$, specifically for $\theta=0^{\circ}$ (navy dashed) and
$~{}0.54^{\circ}$ (magenta dashed) that we summarize in Fig. 6. Upon
application of a perpendicular electric field, the valence band valley Chern
number $C_{V1}=4$ for $\xi=1$ (Fig. 6 (a)) remains within a large set of
parameter spaces of $\Delta\approx 0\sim 35$ meV and $\theta=0\sim
0.8^{\circ}$ and for the larger electric field switches to $C_{V1}=\pm 1,\pm
2$ at zero twist ($\theta=0^{\circ}$) and $C_{V1}=\pm 1,2$ at
$\theta=0.54^{\circ}$. The conduction band, however, shows a nonzero valley
Chern number $C_{C1}=4$ only for the small angles $\theta\approx 0\sim
0.2^{\circ}$ with negative $\Delta\approx-20\sim-40$ meV and $C_{C1}=1$ for
positive $\Delta\approx 0\sim 10$ meV. In another substrate orientation,
$\xi=-1$, (Fig. 6 (b)), the valence band shows $C_{V1}=4$ with $\Delta=0\sim
35$ meV and moves to $C_{V1}=\pm 1,\pm 2$ with a larger positive interlayer
potential difference, while the nonzero conduction band valley Chern number
$C_{C1}=4$ only appears for the negative $\Delta=-40\sim 0$ meV. Hence, a
single 4LG/BN device can have various topological states ($C_{V1}=\pm 1,\pm
2,4$, $C_{C1}=1,4$), easily controlled through small changes in interlayer
potential differences of the order of $|\Delta|<100$ meV.

Interestingly, we also found that rhombohedral stacked $n$LG/BN ($n=3\sim 8$
layer graphene) for both substrate orientations ($\xi=\pm 1$) have a finite
valence band valley Chern number, $C_{V1}=n$ when $\Delta=0$ eV, see Fig. A2.
Similar to 4LG/BN, the conduction band of $\xi=-1$ have zero valley Chern
number, $C_{C1}=0$ while multilayer 3LG/BN$\sim$8LG/BN systems for $\xi=1$
show various valley Chern numbers, $C_{C1}=2,1,1,-2,$ 0, and 0 at $\Delta=0$
eV. Our calculations suggest that rhombohedral $n$LG/BN multilayers will be an
excellent platform to explore a variety of topological phases with flat enough
bands to find strong correlation effects.

Figure 7: (color online) We present the normalized local density of states
(LDOS) of the isolated valence ($V1$) and conduction ($C1$) bands of 4LG/BN at
their van Hove singularities (vHS), which are defined as
$\tilde{D}({\bm{r}},E)=D({\bm{r}},E)/{\rm max}(D({\bm{r}},E))$. The LDOS is
plotted for two twist angles, $\theta=0^{\circ}$ (left) and
$\theta=0.54^{\circ}$ (right), and two orientations of the hBN substrate,
(a)-(b) $\xi=1$ and (c)-(d) $\xi=-1$. We show the normalized LDOS for selected
interlayer potential differences $\Delta=0.0,0.009,-0.009~{}\rm eV$ for
$\xi=1$, and $\Delta=0.0,0.015,-0.011~{}\rm eV$ for $\xi=-1$. The energy value
at the vHS ($E_{vHS}$) is listed at the bottom of each LDOS panel, and we
label the LDOS maxima as AA, AB, or BA according to the local stacking regions
defined at the top right corner of panels (a) and (c). Our results show that
the LDOS distributions at the vHS are influenced by system parameters such as
the substrate orientation $\xi$, the interlayer potential difference $\Delta$,
and the twist angle $\theta$, and exhibit distinct features for the valence
and conduction bands.

###  III.4 Local density of states

In this subsection, we discuss the local density of states (LDOS) of 4LG/BN at
the van Hove singularities (vHS) of the low energy valence and conduction
bands that are most prone to form ordered phases. Knowledge of their
localization centers is helpful for anticipating and understanding the nature
of the ground states that we can expect when Coulomb interactions are
accounted for. We show that the LDOS profiles associated to the valence or
conduction bands can be tuned by varying the system parameters such as
electric fields, the twist angles, and the orientation of the hBN substrate.
Our results are summarized in Fig. 7, where we show the unit renormalized LDOS
function $\tilde{D}({\bm{r}},E)$ for twist angles
$\theta=0^{\circ},\,0.54^{\circ}$ for both orientations of the hBN substrate
$\xi=\pm 1$. We explicitly label with AA, AB, and BA the local stacking
configurations where the LDOS profiles show maxima values. For zero electric
fields and $\xi=1$ substrate orientation we notice that the LDOS accumulates
at the AB local stacking sites, while for $\xi=-1$ the localization center
switches to AA stacking sites for both valence and conduction bands. This
behavior is observed for both $\theta=0^{\circ},\,0.54^{\circ}$ twist angles
considered and can be modified by applying an interlayer potential difference
$\Delta$.

In the following, we discuss in more detail how an interlayer potential
difference can alter the localization properties. The electron accumulation at
the topmost layer favored by $\Delta>0$ should predominantly concentrate at
the low energy $B_{4}$ site, and likewise for $\Delta<0$ will concentrate
mostly at $A_{1}$, and in opposite senses for holes. We will show that
$\Delta$ can modify the degree of localization of the states around particular
local stacking configurations. Let us consider the $\xi=1$ orientation for the
hBN substrate (Fig. 7(a)) and focus on the $\theta=0^{\circ}$ case for sake of
definiteness. The valence band LDOS spread over all stacking configurations
and are delocalized, mildly concentrating at the AB regions when $\Delta=0$
eV. The initially widespread LDOS profile concentrates at the AB sites in the
presence of interlayer potential differences either positive or negative, see
the plots for $\Delta=0.009$, $-0.009$ eV. We thus expect that a finite
$\Delta$ will enhance the chances of triggering a Coulomb driven transition.
For the conduction band, the LDOS concentrate initially at both AB to BA
regions when $\Delta\neq 0$, where AB is slightly favored. A positive
$\Delta=$0.009 eV shifts the carrier densities from AB to predominantly
concentrate at the BA regions and also increases the population of the AA
regions, while for negative $\Delta=-0.009$ eV the LDOS profile spreads almost
uniformly for all possible stacking regions. Hence, vertical fields allow
transitions from a delocalized to AB concentrated LDOS profiles in the valence
bands, or from AB/BA centered profiles to delocalized phases for conduction
bands depending on applied electric fields.

Figure 8: (color online) We calculate the $k$-point-projected LDOS
$|\psi_{V1(C1)}({\bm{k}}:{\bm{r}})|^{2}$ at twist angle $\theta=0^{\circ}$ and
high-symmetry points $\tilde{K}$, $\tilde{K}^{\prime}$, $\tilde{\Gamma}$, and
$\overline{M}$, for two different hBN substrate orientations: (a) $\xi=1$ and
(b) $\xi=-1$. We normalize the values by the maximum value, and define
$\overline{M}$ as $(\tilde{M}+\tilde{M}^{\prime}+\tilde{M}^{\prime\prime})/3$.
We consider selected interlayer potential differences
$\Delta=0.0,0.009,-0.009~{}\rm eV$ for $\xi=1$, and
$\Delta=0.0,0.015,-0.011~{}\rm eV$ for $\xi=-1$, and plot the LDOS for a given
$k$ point that contributes significantly to the LDOS, as labeled in the upper
left corner of each figure. Our results reveal that the LDOS projected onto
$\tilde{K}$, $\tilde{K}^{\prime}$, and $\overline{M}$ tends to localize in one
of the local stackings, AA, AB, or BA. However, the $\tilde{\Gamma}$-projected
LDOS, especially associated with the $B_{4}$ band, shows delocalization for
zero $\Delta$ in the valence bands and negative $\Delta$ in the conduction
bands, regardless of the hBN substrate orientation $\xi=\pm 1$.

In the other substrate orientation, $\xi=-1$, (Fig. 7(b)) with $\Delta=0$, the
LDOS slightly(strongly) concentrate at AA stacking for valence(conduction)
band. Unlike the $\xi=1$ orientation, a positive $\Delta=0.015$ eV makes the
valence band LDOS further spread over all stacking configurations while the
conduction band LDOS from AA stacking concentrates at AB stacking strongly. A
negative $\Delta=-0.011$ eV shifts the conduction band states towards the top-
most layer, especially at $B_{4}$, delocalizing the conduction band LDOS,
while the valence band LDOS concentrating at $A_{1}$ in the bottom layer is
strongly localized at AA.

A twist angle tilts the moire unit cell and affects the localization of the
valence and conduction band LDOS in 4LG/BN. For $\xi=1$ at
$\theta=0.54^{\circ}$ (Fig. 7(c)), similar to the zero twist angle system, the
valence band LDOS spreads over the moire unit cell and the conduction band
LDOS resides at AB and BA slightly favoring AB stacking for
$\Delta=0,~{}0.009~{}$eV, while the valence band LDOS resides on AB stacking
sites and conduction band LDOS spreads within the moire unit cell for
$\Delta=-0.009~{}$eV. A positive $\Delta=0.009$ eV shifts the localization
center of the conduction band LDOS from BA to the AB. However, for $\xi=-1$ at
$\theta=0.54^{\circ}$ (Fig. 7(d)), the carrier densities of valence and
conduction bands move to the AA stacking site regardless of the electric
fields ($\Delta=0,0.015,-0.011$ eV) where the conduction band LDOS at the
negative $\Delta=-0.011$ eV spreads almost uniformly over the moire unit cell
like in the other hBN substrate orientation and twist angles.

In the remaining part of this subsection, we further discuss the
delocalization in the isolated nearly flat bands of 4LG/BN, that has some
analogies to double bilayer graphene Zhang et al. (2021). We found that the
selective sublattice and $k$-point-projected LDOS at the vHS is helfpul for
the study of the delocalization in 4LG/BN. From the sublattice-projected LDOS
(Fig. A4), we found that most of the low-energy states are associated with the
$A_{1}$ site in the bottom-most graphene layer and the $B_{4}$ site in the
top-most layer, and the states from $B_{4}$ are most important for the flat
bands’ delocalization where the conduction or valence band LDOS almost equally
spreads over all stacking configurations. Also, from the layer-dependent
sublattice-projected LDOS calculation of 3LG/BN, 4LG/BN, and 5LG/BN (See Fig.
A5), we verified that the delocalization of the states are associated with the
$B_{n}$ site located far away from the hBN substrate but the degree of
localization shows a non-monotonic behavior as $n$ increases in $n$LG/BN for
both $\xi=\pm 1$. When we analyzed the $k$-point-projected LDOS we found that
the states from the $\Gamma$ are the most delocalized. In Fig. 8, we
illustrate the $k$-point-projected LDOS at the specific $\bm{k}$ points,
$\tilde{\Gamma}$ or $\tilde{M}$, where the vHS states dominantly come from.
Here, $\bar{M}$ is defined as
$\bar{M}=(\tilde{M}+\tilde{M}^{\prime}+\tilde{M}^{\prime\prime})/3$ and we
calculate $|\psi({\bm{k}}:{\bm{r}})|^{2}$ in Eq. (10) as the $k$-point-
projected LDOS by extracting the wave function at the select $k$-points,
${\bm{k}}=\tilde{K},\tilde{K}^{\prime},\tilde{\Gamma},$or
$\tilde{M}^{(\prime,\prime\prime)}$. The $\tilde{\Gamma}$-projected LDOS,
especially associated to the $B_{4}$ site, show delocalization for
zero(negative) $\Delta$ in the valence(conduction) bands regardless of the hBN
substrate orientation $\xi=\pm 1$, while the LDOS contribution from the other
symmetric $k$-points $\tilde{K}$, $\tilde{K}^{\prime}$, and
$\tilde{M}^{(\prime,\prime\prime)}$ are localized at either AA, AB, or BA.

##  IV Summary and Conclusions

In summary, we investigated the parameters giving rise to isolated nearly flat
bands in rhombohedral four-layer aligned to boron nitride (4LG/BN) and
$n$LG/BN multilayer systems using full-bands continuum models. We found that
4LG/BN ($\xi=1$) and 4LG/NB ($\xi=-1$) alignments with zero degree twists have
naturally narrow low energy bandwidths (12-20 meV), and further explored the
parameter space of interlayer potential difference ($\Delta$) and twist angles
($\theta$). The band widths of the 4LG/BN are generally smaller than the
estimated on-site Coulomb repulsion energy $U\approx 25~{}$meV for a large
parameter space up to $\Delta=\pm$ 0.1 eV and $\theta=0\sim 0.8^{\circ}$,
where $U/W>1$ indicates regions of possible onset of Coulomb ordered phases in
4LG/BN. Very narrow band widths ($<5$ meV) are possible with appropriate
$\Delta$ values for a range of twist angles $\theta=0\sim 0.5^{\circ}$ for
both $\xi=\pm 1$.

We also found that isolated flat bands are achievable for specific ranges of
$\theta$ and $\Delta$ where both primary ($\delta_{p}$) and secondary
($\delta_{s}$) gaps are open. In the narrow band regime ($U_{\rm eff}/W\gg
1$), we estimated the screened on-site Coulomb repulsion energy $U_{\rm
eff}=U~{}e^{-{\ell_{m}}/\lambda_{D}}$, taking into account the screening
effect due to band overlap. As a result, we observe strong particle-hole
asymmetric behavior in the 4LG/BN narrow bands. Narrow valence bands for the
$\xi=1$ substrate orientation can be found for a large range of parameters
$\Delta\approx-50\sim 10$ meV and $\theta=0\sim 0.8^{\circ}$, while the
conduction bands flatten only near $\Delta\approx 7\sim 27$ meV and
$\theta\approx 0.5^{\circ}$. In contrast, the $\xi=-1$ alignment has a larger
parameter space for conduction narrow band regime than the valence narrow band
regime. Conduction narrow bands can exist for $\Delta=0\sim 40$ meV and
$\theta=0\sim 0.8^{\circ}$, while valence narrow bands are only possible near
$\Delta\approx 20\sim 25$ meV and $\theta\approx 0.5\sim 0.6^{\circ}$. We
found that the maximum value of $U_{\rm eff}/W$ occurs when
$\theta=0.54^{\circ}$, in the valence band with $\xi=1$ alignment taking
$U_{\rm eff}/W=3.98$, and in the conduction band with $\xi=-1$ alignment
giving rise to $U_{\rm eff}/W=4.7$. Therefore, we suggest that the twist angle
of $\theta\approx 0.5^{\circ}$ for the 4LG/BN device has the chance of
exhibiting Coulomb-induced ordered states for a wider range of perpendicular
electric fields.

Furthermore, we showed that the isolated valence and conduction flat bands
have well-defined valley Chern numbers of $C_{V1}^{\nu=\pm 1}=\pm 4$ for both
$\xi=\pm 1$ and $C_{C1}^{\nu=\pm 1}=\pm 1(0)$ for $\xi=1$($\xi=-1$) even
without any external electric fields. The valence band valley Chern number
remains $C_{V1}^{\nu=\pm 1}=\pm 4$ for both $\xi=\pm 1$ orientations for
$\Delta=0\sim 0.035$ eV and similarly for the conduction bands we have
$C_{C1}^{\nu=\pm 1}=\pm 1$ for $\xi=1$ for $\Delta=0\sim 0.010$ eV. By
increasing $|\Delta|<100$ meV, the valence and conduction bands can have
various topological states ($C_{V1}=\pm 1,\pm 2,4$, $C_{C1}=1,4$). Also, we
found in $n$LG/BN (3$\sim$8 layer graphene) that the valence band valley Chern
number of $n$LG/BN is the same as the number of graphene layers ($C_{C1}=n$)
for both $\xi=\pm 1$ alignments when $\Delta=0$ eV and $\theta=0^{\circ}$.

We have then investigated the local density of states associated with the van
Hove singularities (vHS) of the 4LG/BN flat bands. In the absence of an
interlayer potential difference and twist angle, the valence band states at
the vHS are spread over all local stackings, slightly favoring the AB stacking
for $\xi=1$ and AA stacking for $\xi=-1$. Meanwhile, the conduction band
states for $\xi=-1$ are strongly localized at AA stacking, while for the
$\xi=1$ substrate orientation, they reside mainly at the AB and BA stackings,
with the AB stacking being slightly favored. The position and degree of
localization vary with the substrate orientation, interlayer potential
difference, and twist angle. A negative $\Delta=-0.009$ eV for $\xi=1$ and
$\Delta=-0.011$ eV for $\xi=-1$ enhances the localization strength of the
valence band vHS states without moving the localization center, while a
positive $\Delta=0.009$ eV for $\xi=1$ and $\Delta=0.015$ eV for $\xi=-1$ move
the localization center from AB (AA) to BA (AB) for $\xi=1$ ($\xi=-1$)
substrate orientation.

We can modify the localization of the 4LG/BN nearly flat band states with an
electric field. Through sublattice-projected LDOS calculations, we confirmed
that the states at the vHS mostly come from the low energy non-dimer
sublattice sites, specifically $A_{1}$ from the bottom-most graphene layer
contacting hBN and $B_{4}$ from the top-most layer. The delocalized states are
related with the $B_{4}$ sites in 4LG/BN. We examined the LDOS data for other
multilayers such as 3LG/BN and 5LG/BN and found that the states populating the
farthest $B_{n}$ site away from the hBN substrate tend to be delocalized at
all local stackings, either for the conduction or valence band depending on
the sign of $\Delta$. However, there is seemingly no direct relation between
layer number and localization, considering that the $B_{3}$ in 3LG/BN or
$B_{5}$ in 5LG/BN both show stronger localization compared to 4LG/BN. Further,
we examine the $k$-point projected LDOS in 4LG/BN, which shows that the
delocalized $B_{4}$ state mainly originates from the $\tilde{\Gamma}$ point in
the MBZ.

In summary, the 4LG/BN has narrow bands in a large parameter space of twist
angles and interlayer potential differences, has various associated valley
Chern bands, and the localization of the wave functions can be easily
controlled by applying relatively small interlayer potential differences of
$|\Delta|<0.1$ eV. The low energy states in 4LG/BN are particularly
delocalized compared to other layer numbers that will impact the ground states
when we explicitly account for Coulomb interactions. Given the high tunability
of the system and narrow bandwidths as small as 10 meV we believe that the
4LG/BN and larger $n$ multilayer $n$LG/BN are excellent systems to study flat
band physics provided that the hurdles for preparing rhombohedral graphene
multilayers can be overcome.

##  V Acknowledgments.

This work was supported by the National Research Foundation of Korea (NRF)
with grant number No. NRF-2021R1A3A13045898 for Y. P., No.
NRF-2020R1A2C3009142 for Y. K., No. NRF-2020R1A5A1016518 for J. J., and the
SERB with grant no. SRG/2022/001102 for B. L. C.. We acknowledge computational
support from KISTI Grant No. KSC-2022-CRE-0514 and by the resources of Urban
Big data and AI Institute (UBAI) at UOS. J.J. also acknowledges support by the
Korean Ministry of Land, Infrastructure and Transport(MOLIT) from the
Innovative Talent Education Program for Smart Cities.

## References

  * Hass et al. (2008) J. Hass, F. Varchon, J. E. Millán-Otoya, M. Sprinkle, N. Sharma, W. A. de Heer, C. Berger, P. N. First, L. Magaud, and E. H. Conrad, Phys. Rev. Lett. 100, 125504 (2008).
  * Miller et al. (2009) D. L. Miller, K. D. Kubista, G. M. Rutter, M. Ruan, W. A. de Heer, P. N. First, and J. A. Stroscio, Science 324, 924 (2009).
  * Miller et al. (2010) D. L. Miller, K. D. Kubista, G. M. Rutter, M. Ruan, W. A. de Heer, P. N. First, and J. A. Stroscio, Phys. Rev. B 81, 125427 (2010).
  * Sadowski et al. (2006) M. L. Sadowski, G. Martinez, M. Potemski, C. Berger, and W. A. de Heer, Phys. Rev. Lett. 97, 266405 (2006).
  * De Heer et al. (2010) W. A. De Heer, C. Berger, X. Wu, M. Sprinkle, Y. Hu, M. Ruan, J. A. Stroscio, P. N. First, R. Haddon, B. Piot, et al., Journal of Physics D: Applied Physics 43, 374007 (2010).
  * Brihuega et al. (2012) I. Brihuega, P. Mallet, H. González-Herrero, G. Trambly de Laissardière, M. M. Ugeda, L. Magaud, J. M. Gómez-Rodríguez, F. Ynduráin, and J.-Y. Veuillen, Phys. Rev. Lett. 109, 196802 (2012).
  * Ohta et al. (2012) T. Ohta, J. T. Robinson, P. J. Feibelman, A. Bostwick, E. Rotenberg, and T. E. Beechem, Phys. Rev. Lett. 109, 186807 (2012).
  * Lopes dos Santos et al. (2012) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 86, 155449 (2012).
  * Lopes dos Santos et al. (2007) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
  * Shallcross et al. (2010) S. Shallcross, S. Sharma, E. Kandelaki, and O. A. Pankratov, Phys. Rev. B 81, 165105 (2010).
  * Shallcross et al. (2008) S. Shallcross, S. Sharma, and O. A. Pankratov, Phys. Rev. Lett. 101, 056803 (2008).
  * Landgraf et al. (2013) W. Landgraf, S. Shallcross, K. Türschmann, D. Weckbecker, and O. Pankratov, Phys. Rev. B 87, 075433 (2013).
  * Shallcross et al. (2013) S. Shallcross, S. Sharma, and O. Pankratov, Phys. Rev. B 87, 245403 (2013).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Moon and Koshino (2013) P. Moon and M. Koshino, Phys. Rev. B 87, 205404 (2013).
  * Moon and Koshino (2012) P. Moon and M. Koshino, Phys. Rev. B 85, 195458 (2012).
  * Jung et al. (2014) J. Jung, A. Raoux, Z. Qiao, and A. H. MacDonald, Phys. Rev. B 89, 205414 (2014).
  * San-Jose et al. (2012) P. San-Jose, J. González, and F. Guinea, Phys. Rev. Lett. 108, 216802 (2012).
  * San-Jose and Prada (2013) P. San-Jose and E. Prada, Phys. Rev. B 88, 121408 (2013).
  * Stauber et al. (2013) T. Stauber, P. San-Jose, and L. Brey, New Journal of Physics 15, 113050 (2013).
  * Bistritzer and MacDonald (2010) R. Bistritzer and A. H. MacDonald, Phys. Rev. B 81, 245412 (2010).
  * Wang et al. (2012) Z. Wang, F. Liu, and M. Chou, Nano letters 12, 3833 (2012).
  * Schmidt et al. (2014) H. Schmidt, J. C. Rode, D. Smirnov, and R. J. Haug, Nature communications 5 (2014).
  * Carr et al. (2018) S. Carr, S. Fang, P. Jarillo-Herrero, and E. Kaxiras, Phys. Rev. B 98, 085144 (2018).
  * Koshino et al. (2018) M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys. Rev. X 8, 031087 (2018).
  * Kang and Vafek (2018) J. Kang and O. Vafek, Phys. Rev. X 8, 031088 (2018).
  * Kang and Vafek (2019) J. Kang and O. Vafek, Phys. Rev. Lett. 122, 246401 (2019).
  * Tarnopolsky et al. (2019) G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Phys. Rev. Lett. 122, 106405 (2019).
  * Po et al. (2019) H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Phys. Rev. B 99, 195455 (2019).
  * Cao et al. (2018a) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, et al., Nature 556, 80 (2018a).
  * Kim et al. (2017) K. Kim, A. DaSilva, S. Huang, B. Fallahazad, S. Larentis, T. Taniguchi, K. Watanabe, B. J. LeRoy, A. H. MacDonald, and E. Tutuc, Proceedings of the National Academy of Sciences 114, 3364 (2017).
  * Sharpe et al. (2019) A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * Cao et al. (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Yankowitz et al. (2019) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Cao et al. (2020) Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020).
  * Lopez-Bezanilla and Lado (2020) A. Lopez-Bezanilla and J. L. Lado, Phys. Rev. Research 2, 033357 (2020).
  * Shin et al. (2021) J. Shin, Y. Park, B. L. Chittari, J.-H. Sun, and J. Jung, Phys. Rev. B 103, 075423 (2021).
  * Park et al. (2021) J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 590, 249 (2021).
  * Hao et al. (2021) Z. Hao, A. M. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Science 371, 1133 (2021).
  * Ramires and Lado (2021) A. Ramires and J. L. Lado, Phys. Rev. Lett. 127, 026401 (2021).
  * Qin and MacDonald (2021) W. Qin and A. H. MacDonald, Phys. Rev. Lett. 127, 097001 (2021).
  * Phong et al. (2021) V. o. T. Phong, P. A. Pantaleón, T. Cea, and F. Guinea, Phys. Rev. B 104, L121116 (2021).
  * Kim et al. (2022) H. Kim, Y. Choi, C. Lewandowski, A. Thomson, Y. Zhang, R. Polski, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Nature 606, 494 (2022).
  * Fischer et al. (2022) A. Fischer, Z. A. H. Goodwin, A. A. Mostofi, J. Lischner, D. M. Kennes, and L. Klebl, npj Quantum Materials 7, 5 (2022).
  * Shen et al. (2023) C. Shen, P. J. Ledwith, K. Watanabe, T. Taniguchi, E. Khalaf, A. Vishwanath, and D. K. Efetov, Nature Materials 22, 316 (2023).
  * Yankowitz (2023) M. Yankowitz, Nature Materials 22, 286 (2023).
  * Park et al. (2020) Y. Park, B. L. Chittari, and J. Jung, Phys. Rev. B 102, 035411 (2020).
  * Rademaker et al. (2020) L. Rademaker, I. V. Protopopov, and D. A. Abanin, Phys. Rev. Res. 2, 033150 (2020).
  * Ma et al. (2021) Z. Ma, S. Li, Y.-W. Zheng, M.-M. Xiao, H. Jiang, J.-H. Gao, and X. Xie, Science Bulletin 66, 18 (2021).
  * Xu et al. (2021) S. Xu, M. M. Al Ezzi, N. Balakrishnan, A. Garcia-Ruiz, B. Tsim, C. Mullan, J. Barrier, N. Xin, B. A. Piot, T. Taniguchi, et al., Nature Physics 17, 619 (2021).
  * Chen et al. (2021) S. Chen, M. He, Y.-H. Zhang, V. Hsieh, Z. Fei, K. Watanabe, T. Taniguchi, D. H. Cobden, X. Xu, C. R. Dean, et al., Nature Physics 17, 374 (2021).
  * He et al. (2021) M. He, Y.-H. Zhang, Y. Li, Z. Fei, K. Watanabe, T. Taniguchi, X. Xu, and M. Yankowitz, Nature Communications 12, 4727 (2021).
  * Li et al. (2022) S.-y. Li, Z. Wang, Y. Xue, Y. Wang, S. Zhang, J. Liu, Z. Zhu, K. Watanabe, T. Taniguchi, H.-j. Gao, et al., Nature Communications 13, 4225 (2022).
  * Tong et al. (2022) L.-H. Tong, Q. Tong, L.-Z. Yang, Y.-Y. Zhou, Q. Wu, Y. Tian, L. Zhang, L. Zhang, Z. Qin, and L.-J. Yin, Phys. Rev. Lett. 128, 126401 (2022).
  * Chebrolu et al. (2019) N. R. Chebrolu, B. L. Chittari, and J. Jung, Phys. Rev. B 99, 235417 (2019).
  * Liu et al. (2020) X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, Y. Ronen, H. Yoo, D. Haei Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, et al., Nature 583, 221 (2020).
  * Lee et al. (2019) J. Y. Lee, E. Khalaf, S. Liu, X. Liu, Z. Hao, P. Kim, and A. Vishwanath, Nature Communications 10, 5333 (2019).
  * Shen et al. (2020) C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, K. Watanabe, et al., Nature Physics 16, 520 (2020).
  * Choi and Choi (2019) Y. W. Choi and H. J. Choi, Phys. Rev. B 100, 201402 (2019).
  * Koshino (2019) M. Koshino, Phys. Rev. B 99, 235406 (2019).
  * Li et al. (2020) X. Li, F. Wu, and S. Das Sarma, Phys. Rev. B 101, 245436 (2020).
  * Sinha et al. (2020) S. Sinha, P. C. Adak, R. S. Surya Kanthi, B. L. Chittari, L. D. V. Sangani, K. Watanabe, T. Taniguchi, J. Jung, and M. M. Deshmukh, Nature Communications 11, 5548 (2020).
  * Song et al. (2013) J. C. W. Song, A. V. Shytov, and L. S. Levitov, Phys. Rev. Lett. 111, 266801 (2013).
  * Javvaji et al. (2020) S. Javvaji, J.-H. Sun, and J. Jung, Phys. Rev. B 101, 125411 (2020).
  * Naik and Jain (2018) M. H. Naik and M. Jain, Phys. Rev. Lett. 121, 266401 (2018).
  * Naik et al. (2020) M. H. Naik, S. Kundu, I. Maity, and M. Jain, Phys. Rev. B 102, 075413 (2020).
  * Wu et al. (2019) F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, Phys. Rev. Lett. 122, 086402 (2019).
  * Zhan et al. (2020) Z. Zhan, Y. Zhang, P. Lv, H. Zhong, G. Yu, F. Guinea, J. A. Silva-Guillén, and S. Yuan, Phys. Rev. B 102, 241106 (2020).
  * Wang et al. (2020) L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, et al., Nature Materials 19, 861 (2020).
  * Zhang et al. (2020) Z. Zhang, Y. Wang, K. Watanabe, T. Taniguchi, K. Ueno, E. Tutuc, and B. J. LeRoy, Nature Physics 16, 1093 (2020).
  * Shabani et al. (2021) S. Shabani, D. Halbertal, W. Wu, M. Chen, S. Liu, J. Hone, W. Yao, D. N. Basov, X. Zhu, and A. N. Pasupathy, Nature Physics 17, 720 (2021).
  * Rademaker (2022) L. Rademaker, Phys. Rev. B 105, 195428 (2022).
  * Wong et al. (2015) D. Wong, Y. Wang, J. Jung, S. Pezzini, A. M. DaSilva, H.-Z. Tsai, H. S. Jung, R. Khajeh, Y. Kim, J. Lee, et al., Phys. Rev. B 92, 155409 (2015).
  * Kerelsky et al. (2019) A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, et al., Nature 572, 95 (2019).
  * Choi et al. (2019) Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, et al., Nature Physics 15, 1174 (2019).
  * Moon and Koshino (2014) P. Moon and M. Koshino, Phys. Rev. B 90, 155406 (2014).
  * Yankowitz et al. (2012a) M. Yankowitz, J. Xue, D. Cormode, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, P. Jacquod, and B. J. LeRoy, Nature Physics 8, 382 (2012a).
  * Yankowitz et al. (2012b) M. Yankowitz, J. Xue, D. Cormode, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, P. Jacquod, and B. J. LeRoy, Nature Physics 8, 382 (2012b).
  * Dean et al. (2013) C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, et al., Nature 497, 598 (2013).
  * Ponomarenko et al. (2013) L. A. Ponomarenko, R. V. Gorbachev, G. L. Yu, D. C. Elias, R. Jalil, A. A. Patel, A. Mishchenko, A. S. Mayorov, C. R. Woods, J. R. Wallbank, et al., Nature 497, 594 (2013).
  * Amet et al. (2013) F. Amet, J. R. Williams, K. Watanabe, T. Taniguchi, and D. Goldhaber-Gordon, Phys. Rev. Lett. 110, 216601 (2013).
  * Hunt et al. (2013) B. Hunt, J. Sanchez-Yamagishi, A. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, et al., Science 340, 1427 (2013).
  * Woods et al. (2014) C. Woods, L. Britnell, A. Eckmann, R. Ma, J. Lu, H. Guo, X. Lin, G. Yu, Y. Cao, R. Gorbachev, et al., Nature physics 10, 451 (2014).
  * Jung et al. (2015) J. Jung, A. M. DaSilva, A. H. MacDonald, and S. Adam, Nature communications 6, 6308 (2015).
  * San-Jose et al. (2014) P. San-Jose, A. Gutiérrez-Rubio, M. Sturla, and F. Guinea, Phys. Rev. B 90, 075428 (2014).
  * Leconte et al. (2017) N. Leconte, J. Jung, S. Lebègue, and T. Gould, Phys. Rev. B 96, 195431 (2017).
  * Kim et al. (2018) H. Kim, N. Leconte, B. L. Chittari, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Jung, and S. Jung, Nano Letters 18, 7732 (2018).
  * Leconte and Jung (2020) N. Leconte and J. Jung, 2D Materials 7, 031005 (2020).
  * Chen et al. (2019a) G. Chen, L. Jiang, S. Wu, B. Lyu, H. Li, B. L. Chittari, K. Watanabe, T. Taniguchi, Z. Shi, J. Jung, et al., Nature Physics 15, 237 (2019a).
  * Chen et al. (2019b) G. Chen, A. L. Sharpe, P. Gallagher, I. T. Rosen, E. J. Fox, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, et al., Nature 572, 215 (2019b).
  * Chittari et al. (2019) B. L. Chittari, G. Chen, Y. Zhang, F. Wang, and J. Jung, Phys. Rev. Lett. 122, 016401 (2019).
  * Zhang et al. (2019) Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T. Senthil, Phys. Rev. B 99, 075127 (2019).
  * Chen et al. (2020) G. Chen, A. L. Sharpe, E. J. Fox, Y.-H. Zhang, S. Wang, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, et al., Nature 579, 56 (2020).
  * González et al. (2021) D. A. G. González, B. L. Chittari, Y. Park, J.-H. Sun, and J. Jung, Phys. Rev. B 103, 165112 (2021).
  * Yang et al. (2022) J. Yang, G. Chen, T. Han, Q. Zhang, Y.-H. Zhang, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, et al., Science 375, 1295 (2022).
  * Calderón et al. (2022) M. J. Calderón, A. Camjayi, and E. Bascones, Phys. Rev. B 106, L081123 (2022).
  * Min and MacDonald (2008) H. Min and A. H. MacDonald, Progress of Theoretical Physics Supplement 176, 227 (2008).
  * D. et al. (2010) X. D., C. M.-C., and N. Q., Rev. Mod. Phys. 82, 1959 (2010).
  * Zhang et al. (2021) C. Zhang, T. Zhu, S. Kahn, S. Li, B. Yang, C. Herbig, X. Wu, H. Li, K. Watanabe, T. Taniguchi, et al., Nature Communications 12, 2516 (2021).
  * Jung and MacDonald (2013) J. Jung and A. H. MacDonald, Phys. Rev. B 87, 195450 (2013).

## Appendix

## A1. F1G0 tight binding model for multilayer graphene

The effective tight-binding F1G0 model Jung and MacDonald (2013) utilized in
this study is based on grouping the hopping amplitudes by the hopping
distances, where the effective Hamiltonian only includes the structure factor
functions of $g_{0}({\bm{k}})$ for the zero distance intra-sublattice hopping
and $f_{1}({\bm{k}})$ for the first nearest neighbor inter-sublattice hopping.
To construct the ABC stacked rhombohedral multilayer graphene Hamiltonian, we
adopted the trilayer graphene (3LG) hopping parameters since they provide a
highly accurate tight-binding fit to the first-principles calculations
bandstructure in the low energy range near the charge neutrality point (CNP).
The resulting band structures of the F1G0 tight-binding model for multilayer
graphene, as shown in Fig. A1, are in good agreement with those obtained from
density functional theory (DFT) calculations, within the energy range of
$|E|<0.05$eV and the $k$ vector range of
$|{\bm{k}}|<0.05~{}\left(\frac{4\pi}{3a_{G}}\right)$ measured from the Dirac
points ${\bm{K}}_{\nu=\pm 1}$. We also extended the effective tight-binding
model Hamiltonian to 4LG, 5LG, 6LG, 7LG, and 8LG by modifying the adopted 3LG
hopping parameters. To obtain the DFT band structures, we used $Quantum\
ESPRESSO$ with local density approximation (LDA) parametrization (C.pz-
rrkjus.UPF), where the graphene lattice constant and the interlayer separation
are set to $a=2.46\AA$ and $c=3.35\AA$. A 60 Ry energy cutoff and a $30\times
30\times 1$ Monkhorst-Pack $k$ point grid were used.

Figure A1: (color online) We compare the F1G0 tight-binding model bands
(TB(F1G0)) to those obtained from first-principles calculations within the
local density approximation (DFT(LDA)) for 3LG-8LG along the $k$ path from M
to K to $\Gamma$. The bands from both methods are in good agreement in the
energy range of $-0.05\sim 0.05~{}$eV and the $k$-space range of
$|{\bm{k}}|<0.05\left(\frac{4\pi}{3a}\right)$ near the Dirac points
${\bm{K}}_{\nu=1}$.

## A2. Bandwidth in multilayer graphene boron-nitride moire superlattices
($n$LG/BN)

Our investigation revealed that the low-energy valence and conduction bands in
4LG/BN exhibit significantly flatter bandwidths compared to those in 3LG/BN
Chittari et al. (2019). To gain deeper insights into this trend, we examined
the bandwidths for $n$LG/BN where $n$ varies from 3 to 8. For this purpose, we
added the effective moiré potential $H^{\xi}_{M}({\bm{k}}:{\bm{r}})$ described
in Eq.(6) to the F1G0 Hamiltonian of $n$LG, as detailed in Appendix A1. Our
analysis is summarized in Fig.A2, where we show (a)-(b) the band structures of
3LG/BN, 4LG/BN, and 5LG/BN and (c)-(d) the bandwidths of the low-energy
valence and conduction bands as a function of the number of graphene layers
for both orientations of the BN substrate ($\xi=1$). We observed a nonzero
valley Chern number $C^{+}_{V1}=+n$ for the valence bands in $n$LG/BN for both
$\xi=\pm 1$, matching with the number of graphene layers. Meanwhile, the
conduction bands for $\xi=1$ have a nonzero valley Chern number
$C^{+}_{C1}=+2$ for 3LG/BN and $C^{+}_{C1}=+1$ for 4LG/BN and 5LG/BN, while
$C^{+}_{C1}=0$ for $\xi=-1$. Our analysis also revealed that the bandwidths of
the valence and conduction bands decrease from $W\approx 0.03\sim 0.04~{}$eV
(3LG/BN) to $W\approx 0.005\sim 0.010$ eV (5LG/BN) and remain unchanged for
further increases in the number of graphene layers up to eight (8LG/BN), as
shown in Fig.A2 (c)-(d). Among the $n$LG/BN moiré superlattices, the 5LG/BN
with a narrow bandwidth of $W\approx 0.005\sim 0.010$ eV is the most
attractive candidate, but the 4LG/BN with a bandwidth of $W\approx 0.01$ eV is
still one of the most attractive candidates. This is because, as the number of
layers increases, it becomes more challenging to prepare the samples.

Figure A2: (color online) We present the band structure of $n$LG/BN, where
$n=3\sim 5$ layers of graphene on (a) $\xi=1$ and (b) $\xi=-1$ hBN substrates.
The band structures for two valleys are shown with solid ($\nu=1$) and dashed
($\nu=-1$) lines. We observe that the valley Chern number for the $n$LG/BN
valence band matches the number of layers, $C^{\nu=\pm 1}_{V1}=\pm n$, while
those for the conduction bands are $C_{C1}=2$ for the 3LG/BN and $C_{C1}=1$
for the 4LG/BN and 5LG/BN. To investigate the trend in bandwidths of low
energy valence and conduction bands, we plot the bandwidths as a function of
the number of graphene layers for $n$LG/BN, where $n=3\sim 8$. In (c) and (d),
we show the bandwidths for $\xi=1$ and $\xi=-1$, respectively. We observe that
the bandwidths of the valence and conduction bands decrease from $W\approx
40~{}$meV (3LG/BN) to a minimum value of $W\approx 5~{}$meV (5LG/BN) for
$\xi=\pm 1$. The bandwidths remain constant when further increasing the number
of graphene layers for $n$LG/BN, where $n=6\sim 8$.  Figure A3: (color online)
We investigate the impact of interlayer potential difference $\Delta$ on the
band structure of 4LG/BN for the two hBN substrate orientations ($\xi=\pm 1$).
We consider two twist angles (left: $\theta=0^{\circ}$, right:
$\theta=0.54^{\circ}$) and three values of $\Delta$ for each orientation:
(a)-(c) $\Delta=0.0$ eV, $\Delta=0.009$ eV, and $\Delta=-0.009$ eV for
$\xi=1$, and (d)-(f) $\Delta=0.0$ eV, $\Delta=0.015$ eV, and $\Delta=-0.011$
eV for $\xi=-1$. We plot the band structure for both valleys, denoted by
solid/dashed lines for $\nu=1$ and $\nu=-1$. Our results show that the
bandwidth of the valence band decreases with increasing twist angle and
interlayer potential difference, while the valley Chern numbers for each low
energy valence (conduction) band are unaffected by the twist angle at a given
interlayer potential difference.  Figure A4: (color online) We present the
sublattice-projected local density of states (LDOS) of 4LG/BN for various
interlayer potential differences ($\Delta$) at the van Hove singularity (vHS)
of low-energy valence and conduction flat bands. Specifically, we consider (a)
$\Delta=0.0,0.009,-0.009~{}$eV for $\xi=1$, and (b)
$\Delta=0.0,0.015,-0.011~{}$eV for $\xi=-1$, at $\theta=0^{\circ}$. We also
show the corresponding results for $\theta=0.54^{\circ}$ in panels (c) and
(d). Our analysis reveals that the dominant states in the LDOS stem from
$B_{4}$ and $A_{1}$ sublattice sites, as denoted by the black-filled and red-
open triangle symbols, respectively. Furthermore, the dominant $A_{1}$ states
are primarily localized at AB (AA) stacking sites for $\xi=1$ ($\xi=-1$),
while the $B_{4}$ states spread within the moire unit cell.

## A3. Effects of a twist angle

Our study demonstrates that a small twist angle $\theta=~{}0.54^{\circ}$
between graphene and boron nitride has a significant impact on the 4LG/BN band
structure. Specifically, we observe that this twist angle leads to the opening
of a secondary band gap for both conduction and valence bands, compared to the
$\theta=0^{\circ}$ case. These changes are facilitated by the presence of
finite interlayer potential differences. To further illustrate the effect of
the twist angle and interlayer potential differences, we present the band
structures of 4LG/BN in Fig. A3 for twist angles of $\theta=0^{\circ}$ and
$\theta=0.54^{\circ}$, with select interlayer potential differences of
$\Delta=0,0.009,-0.009$ eV for $\xi=1$ and $\Delta=0,0.015,-0.011$ eV for
$\xi=-1$. Notably, we find that the conduction (valence) band isolation for
$\xi=1$ ($\xi=-1$) is present at $\theta=0.54^{\circ}$ in the presence of a
small positive (negative) $\Delta$, which is difficult to achieve at
$\theta=0^{\circ}$. This underscores the critical role of the small finite
twist angle in achieving band isolation remaining the band flatness. In
addition, we discover that the twist angle $\theta=0.54^{\circ}$ reduces the
asymmetry between the two minivalleys at $\tilde{K}$ and $\tilde{K}^{\prime}$
without modifying their topological phases. This adds another intriguing
feature to the impact of the twist angle on the electronic properties of
4LG/BN.

Figure A5: (color online) We present the normalized local density of states
(LDOS) at the van Hove singularity (vHS) of the conduction (C1) and valence
(V1) bands for 3LG/BN, 4LG/BN, and 5LG/BN for two different hBN orientations:
(a) $\xi=1$ and (b) $\xi=-1$. Interestingly, we find no clear trend in the
localization of states with increasing number of graphene layers.
Additionally, we compare the sublattice-projected LDOS of these systems
calculated at the vHS of flat bands for the (c) $\xi=1$ and (d) $\xi=-1$ hBN
orientations. Our analysis reveals that the dominant states originate from the
low energy $B_{n}$ and $A_{1}$ sublattice sites. We find that the $B_{n}$
states far from the hBN susbstrate exhibit delocalization in the case of
$n$LG/BN.

## A4. Sublattice projected local density of states

To understand the delocalization observed in 4LG/BN, we analyzed the
sublattice-projected LDOS of the low-energy flatbands. Our results, shown in
Fig. A4, indicate that the localized states originate from non-dimer low
energy sites, such as $A_{1}$ and $B_{4}$. In particular, $A_{1}$ states are
well localized at AB(AA) stacking sites for $\xi=+1$($\xi=-1$), while $B_{4}$
states exhibit delocalization. The population of the sublattice-projected LDOS
depends on the interlayer potential difference ($\Delta$), with more $B_{4}$
states occupying the valence bands for zero and positive $\Delta$, and the
conduction bands for negative $\Delta$. To confirm that the delocalization of
$B_{4}$ states is related to the distance from the hBN substrate, we
calculated the sublattice-projected LDOS for 3LG/BN, 4LG/BN, and 5LG/BN, as
shown in Fig. A5. Our analysis reveals that valence band states are spread
over the moire unit cell for all systems, with the delocalized states mostly
originating from $B_{3}$, $B_{4}$, and $B_{5}$ for 3LG/BN, 4LG/BN, and 5LG/BN,
respectively. However, the extent of localization is not linear with the
distance from the hBN substrate, as the delocalized states of 4LG/BN are much
more spread over the moire unit cell than those of 5LG/BN. This result
suggests that the delocalization of $B_{4}$ states is influenced not only by
the distance from the hBN substrate but also by other factors, such as the
interlayer potential difference.