[supp-]supporting

# Moiré flat bands and antiferroelectric domains in lattice relaxed twisted
bilayer hexagonal boron nitride under perpendicular electric fields

Fengping Li  Department of Physics, University of Seoul, Seoul 02504, Korea
Dongkyu Lee  Department of Physics, University of Seoul, Seoul 02504, Korea
Department of Smart Cities, University of Seoul, Seoul 02504, Korea  Nicolas
Leconte  Department of Physics, University of Seoul, Seoul 02504, Korea
Srivani Javvaji  Department of Physics, University of Seoul, Seoul 02504,
Korea  Jeil Jung<EMAIL_ADDRESS>Department of Physics, University of
Seoul, Seoul 02504, Korea  Department of Smart Cities, University of Seoul,
Seoul 02504, Korea

###### Abstract

Local interlayer charge polarization of twisted bilayer hexagonal boron
nitride (t2BN) is calculated and parametrized as a function of twist angle and
perpendicular electric fields through tight-binding calculations on lattice
relaxed geometries Lattice relaxations tend to increase the bandwidth of the
nearly flat bands, where widths smaller than $\sim 1$ meV are expected for
$\theta\leq 1.08^{\circ}$ for parallel BN/BN alignment, and for
$\theta<1.5^{\circ}$ for the antiparallel BN/NB alignment. Local interlayer
charge polarization maxima of $\sim$2.6 pC/m corresponding to interlayer
electron density differences of $\sim 1.3\times 10^{12}\rm cm^{-2}$ are
expected at the AB and BA stacking sites of BN/BN aligned t2BN in the long
moire period limit for $\theta\ll 1^{\circ}$, and evolves non-monotonically
with a maximum of $\sim$3.5 pC/m at $\theta=1.6^{\circ}$ before reaching
$\sim$2 pC/m for $\theta=6^{\circ}$. The electrostatic potential maxima due to
the t2BN moiré patterns are overall enhanced by $\sim 20\%$ with respect to
the rigid system assuming potential modulation depths of up to $\sim 300$ mV
near its surface. In BN/BN aligned bilayers the relative areas of the AB or BA
local stacking regions can be expanded or reduced through a vertical electric
field depending on its sign.

##  I  Introduction

Twisted bilayer and multilayer graphene typically have their magic angles
between $1^{\circ}$–$2^{\circ}$ [1, 2, 3, 4] where the low energy bandwidths
become nearly flat. When the Fermi level is tuned to lie within those bands
they often give rise to correlated phases [5, 6, 7, 8, 9, 10, 11, 12],
superconducting phases [13, 14, 15], and phase transition cascades as a
function of carrier density [11, 12]. Twisted gapped Dirac models [16] that
approximate some types of twisted semiconductor transition-metal
dichalcogenides (TMDs) [17, 18, 19, 20, 21, 22, 23, 24, 25, 26] and twisted
bilayer hexagonal boron nitride (t2BN) [27] are expected to give rise to
narrower low energy bands prone to correlation effects in larger twist angles
than the ones from the graphene family. In particular, t2BN has drawn
attention thanks to its interfacial local layer ferroelectricity whose charge
polarization can be switched by perpendicular electric fields [28, 29, 30, 31,
32, 33, 34] by modifying the stacking geometry. Moreover, the moiré potential
resulting from interlayer charge polarization offers an effective means to
engineer the properties of adjacent 2D material layers [35]. However, most
tight-binding (TB) electronic structures currently rely on highly approximate
models. For example, their TB hopping terms do not often distinguish the
atomic species and assume a simplified distance dependence [36], or utilize
interlayer hopping terms essentially similar to those of graphite [37].
Moreover, the geometry for the atomic structure is often assumed rigid [27,
24].

In this work, we report the atomic and electronic structure of t2BN using
models informed by local density approximation (LDA) density functional theory
(DFT) to systematically improve the theoretical framework to describe twisted
bilayer and multilayer hexagonal boron nitride (hBN) systems. Specifically, we
use exact exchange random phase approximation (EXX-RPA)-fitted force fields in
LAMMPS [38, 39] and calculate the electronic structure distinguishing the
intralayer TB Hamiltonian modeled with distant neighbor hopping terms, and the
interlayer coupling modeled through two-center (TC) approximations [40] that
we call the hybrid TC (HTC) model. More specifically, the TC model introduced
for twisted bilayer graphene (tBG) [41, 42] is reparametrized to account for
the species-specific interlayer interactions between B and N atoms. This HTC
model reproduces more reliably the DFT band gap, the bandwidth, and layer-
resolved local charge distributions for different twist angles.

We will show that while lattice relaxations in t2BN tend to widen the band
widths of the low energy nearly flat bands, they reduce below $\sim 1$ meV for
twist angles smaller than $1.08^{\circ}$ and $1.5^{\circ}$ for the parallel
BN/BN and anti-parallel alignments BN/NB respectively. These two types of
alignments have different lattice relaxation behaviors and responses to a
perpendicular electric field that alters the relative areas of different local
stacking configurations. The external electric field has been simulated by
adding layer-dependent onsite energy terms in the HTC model, inferred from the
forces acting on the atoms during the LAMMPS calculation, offering an
effective means to adjust the electronic properties of the system due to the
changes in the moiré stacking domain sizes. Moreover, significant interlayer
charge polarization domains in parallel BN/BN t2BN have been parametrized
through a third-order harmonic expansion, which allows modeling of the
electrostatic potential induced in adjacent 2D material layers.

The paper is organized as follows. Section II introduces the methodology for
the TB electronic structure parametrization and the atomic structure energy
minimization calculations. Section III presents the lattice reconstruction
effects and electronic structure results including band structures in the
absence and presence of an electric field as well as charge transfer and
induced moiré potential estimates, and in section IV we summarize our main
findings.

##  II Methodology: systems and methods

###  II.1 Systems

In the following, we describe the t2BN systems considered in our study. Fig.
1(a) illustrates the moiré pattern formed by twisting the two layers by a
value of $\theta=3.89^{\circ}$, starting from parallel zero degree stacked
BN/BN and choosing the rotation center at the first atomic site of AA stacking
111We choose the rotation axis to be at the first atom in the small unit cell
that defines the local stacking. in which boron and nitrogen atoms in one
layer lie on top of the boron and notrogen atoms in the other layer.
Similarly, the BN/NB moiré pattern can be constructed starting from anti-
parallel BN/NB with the AA′ local stacking axis as the rotation center, where
the boron and nitrogen atoms in one layer lie on top of the nitrogen and boron
atoms in the other layer (BN/NB). The high symmetry local stacking domain
regions appearing in the moiré pattern of BN/BN and BN/NB systems are shown in
Fig. 1(b). A gradual variation of stacking from AA ($\rm AA^{\prime}$) to AB
($\rm AB^{\prime}$) and BA($\rm BA^{\prime}$) domain regions can be
recognized, reminiscent of the well-known moiré pattern strains in tBG moiré
systems where the unstable high energy local stacking regions shrink and the
low energy regions expand. In this work, we use the approach outlined in Ref.
[42] to find commensurate cells for rotation angles from $21.78^{\circ}$ down
to $1.08^{\circ}$. Fig. 1(c) displays the extended Brillouin zone for the
moiré system resulting from a twist angle $\theta$ where ${\bm{b}}_{1}$ and
${\bm{b}}_{2}$ are the reciprocal lattice vectors. The electronic band
structures in our calculations are shown for variations of the
$\mathbf{\Gamma}$-$\mathbf{M}$-$\mathbf{K}$ path of the hexagonal MBZ region
in reciprocal space while the local density of states is plotted at each of
these high-symmetry points.

Figure 1: Schematic representation of t2BN and local stacking geometries. (a)
Real space moiré pattern for BN/BN with $\theta=3.89^{\circ}$. The gray
rhombus represents one unit cell of the moiré pattern. $L$ and $R$ correspond
to the moiré length and the radius of the AA moiré stacking region. (b)
Illustration of these high-symmetry stacking regions for BN/BN with AA-
stacking rotation center (top row) and BN/NB with AA′-stacking rotation center
(bottom row). In both cases, the rotation center is located exactly where the
sub-lattices of the bottom and top layers coincide. Pink and blue balls denote
the B and N atoms, respectively, and the top (bottom) layer is denoted by
smaller (bigger) balls. (c) Schematic illustration of the moiré Brillouin zone
(MBZ) in t2BN. The blue and pink hexagons denote the Brillouin zone
corresponding to the top and bottom hBN layers, respectively. The smaller gray
hexagons represent the MBZ formed by connecting the high symmetry points of
the BZs of the constituent layers (e.g. $\mathbf{K_{1}}$ and
$\mathbf{K_{2}}$). $\theta$ here matches the real space twist angle.

###  II.2 Methods

The atomistic calculations of the relaxed lattice structures are combined with
electronic structure calculations within the TB approximation. The atomic
positions are relaxed using LAMMPS [38, 39], using the ExTeP (Extended Tersoff
Potential) potentials for intralayer boron and nitrogen interactions [44],
while the interlayer interactions are described by the Dihedral Registry-
Dependent Interlayer Potential (DRIP) [45], reparameterized using EXX-RPA data
[46], that can be used for both graphene and hBN interactions [42]. The first-
principles calculations based on DFT performed with the software package
Quantum ESPRESSO (QE) are interfaced with Wannier90 to construct the Wannier
functions [47, 48, 49]. We use the norm-conserving Perdew-Zunger LDA
pseudopotential [50], a $42\times 42\times 1$ k-points sampling density and
$70$ Ry plane wave basis cutoff for the microscopic calculations of the
bilayer with four atoms per unit cell. A $30$ Å perpendicular direction box
size for each bilayer leaves enough vacuum space that suppresses the vertical
interlayer coupling between the repeated cells [51]. The electronic structures
and local densities for t2BN are obtained using a reparametrized hybrid TC
model [41, 52, 42], hereafter indicated as the HTC model, where the intralayer
hopping terms are characterized by considering up to six nearest neighbors for
intra- and intersublattices [53, 54, 55, 56], while we use the TC
approximation for the interlayer hopping terms.

For the HTC model parametrization, considering that the low energy conduction
and valence bands (CB and VB) mostly stem from $p_{z}$ orbitals it is
sufficient to use ten localized Wannier orbitals including four $p_{z}$
orbitals at the atomic sites and six $\sigma$ orbitals at the bonds. The
calculated hopping amplitudes for intra-sublattice (B to B atom or N to N
atom) and inter-sublattice (B to N atom) are listed in Table 1 and are closely
similar to the ones in the model presented in Ref. [53].

Table 1:  Six nearest neighbor intralayer hopping amplitudes in the top or bottom layer in the HTC model as defined in the text, where we take LDA values for the on-site energies of boron and nitrogen are $1.7666$ eV and $-2.1843$ eV, respectively. The indices in the first column refer to the $n$-th nearest neighbors as defined in Ref. [54]. | $t_{\text{BB}}$ | $t_{\text{NN}}$ | $t_{\text{BN}}$
---|---|---|---
$0$ | $1.7666$ | $-2.1843$ | $0.0000$
$1$ | $0.0053$ | $0.1923$ | $-2.7205$
$2$ | $0.0223$ | $0.0193$ | $-0.2102$
$3$ | $-0.0483$ | $-0.0373$ | $0.0797$
$4$ | $-0.0029$ | $-0.0027$ | $0.0047$
$5$ | $-0.0033$ | $0.0000$ | $-0.0081$
$6$ | $0.0002$ | $-0.0009$ | $0.0155$

The diagonal intralayer Hamiltonian elements, corresponding to row 0 in Table
1, are simulated taking into account the interlayer displacement-dependent
Hamiltonian matrix element developed in Refs. [54, 55] whose sliding-dependent
maps are summarized in Fig. S.1 of the Supplemental Material and whose
expressions are reminded here as

$H_{\alpha\beta}({\bm{d}})=C_{0\alpha\beta}+2C_{1\alpha\beta}\operatorname{Re}\left[f({\bm{d}})\exp(i\varphi_{\alpha\beta})\right],\\\
$ (1)

with $\alpha\beta$ = BtBt, NtNt, BbBb or NbNb and

$\displaystyle
f({\bm{d}})=\exp(-iG_{1}d_{y})+2\exp\left(i\frac{G_{1}d_{y}}{2}\right)\cos\left(\frac{\sqrt{3}}{2}G_{1}d_{x}\right)$

where ${d}_{x}$ and ${d}_{y}$ are the vector-components of the local stacking
registry vector [51, 57]

${\bm{d}}({\bm{r}})=(\mathcal{R}(\theta)-1){\bm{r}}+{\bm{u}}({\bm{r}})$ (2)

where $G_{1}=4\pi/\sqrt{3}a$, taking $a$ as the reference lattice constant,
where $\mathcal{R}$ is the rotation operator acting on an atom at site
${\bm{r}}$ and where ${\bm{u}}({\bm{r}})$ captures the displacement of an atom
due to lattice reconstruction. The three parameters controlling this harmonic
approximation (HA) are summarized in Table 2.

Table 2:  The numerical values of $C_{0}$, $C_{1}$ and $\varphi$ in Eq. (1), where ${\rm B_{t}}$, ${\rm B_{b}}$, ${\rm N_{t}}$ and ${\rm N_{b}}$ denote respectively the top layer boron, bottom layer boron, top layer nitrogen, and bottom layer nitrogen atoms. | ${\rm B_{t}}$ | ${\rm B_{b}}$ | ${\rm N_{t}}$ | ${\rm N_{b}}$
---|---|---|---|---
$C_{0}$ | $1.7157$ | $1.7157$ | $-2.2560$ | $-2.2560$
$C_{1}$ | $0.0099$ | $0.0150$ | $0.0150$ | $0.0150$
$\varphi$ | $0.5437$ | $-0.5437$ | $-0.5437$ | $-0.6525$

For relaxed systems we account for variable intralayer bond length-dependence
of the hopping terms through [58, 42, 53]

$t_{\alpha\beta}({r}_{ij})=t_{\alpha\beta}({r}_{0,ij})\exp\left[-2.45\left(\frac{{r}_{ij}-{r}_{0,ij}}{{r}_{0,ij}}\right)\right]$
(3)

where $t_{\alpha\beta}({r}_{0,ij})$ are the hopping terms of the rigid
lattice’s in-plane intralayer interatomic distances ${r}_{0,ij}$ between $i$
and $j$ atoms and ${r}_{ij}$ is the relaxed distance. The $\alpha$ and $\beta$
refer to the sublattices in a single layer hBN. We used a lattice constant of
$2.4795$ Å for the reference rigid-structure DFT calculations to feed our
model but the electronic structure corresponding to lattice relaxed structures
which are close to the experimental lattice constant of $a=2.504$ Å [59] can
be modeled through Eq. (3).

After establishing the intralayer hopping terms, the interlayer hopping terms
are determined using the Slater-Koster approximation [60], which is expressed
as

$t_{\text{inter}}({r}_{ij})={n^{2}_{z}}V_{pp\sigma}({r}_{ij})+\left(1-{n^{2}_{z}}\right)V_{pp\pi}({r}_{ij})$
(4)

where ${n_{z}}=\ z_{ij}/{r}_{ij}$ with $z_{ij}$ the z-coordinate of ${r_{ij}}$
and the coefficients used by the original two center approximation model [41]
for graphite are given by

$V_{\mathrm{pp}\pi}({r}_{ij})=\gamma_{0}\exp\left[q_{\pi}\left(1-\frac{{r}_{ij}}{a_{\rm
BN}}\right)\right],$ (5)

$V_{\mathrm{pp}\sigma}({r}_{ij})=\gamma_{1}\exp\left[q_{\sigma}\left(1-\frac{{r}_{ij}}{c_{\perp}}\right)\right].$
(6)

We adopt this form for the hBN using $a_{\text{BN}}=a/\sqrt{3}=1.43~{}$Å for
the intralayer nearest neighbor distance of hBN, $c_{\perp}=3.261~{}$Å for the
vertical interlayer distance, and $\gamma_{0}=-2.7$ eV for the first nearest
coupling term. The second nearest neighbor hopping term $\gamma_{0}^{\prime}$
is set to $0.1\gamma_{0}$, which is used to fix the ratio of $q_{\pi}/a$ such
that

$\frac{q_{\sigma}}{c_{\perp}}=\frac{q_{\pi}}{a_{BN}}=\frac{\ln\left(\gamma_{0}^{\prime}/\gamma_{0}\right)}{a_{\text{BN}}-a}$
(7)

following the conventions and values from Ref. [41] for tBG.

Figure 2:  Benchmark comparison of TB and DFT Hamiltonian and associated
bands, as well as sliding-dependent energy landscapes, for the 4-atom
commensurate cells. (a) In-plane distance-dependence of the nearest-neighbor
interlayer hopping parameters used in the HTC model defined in the text.
$t_{\rm B_{t}\rm B_{b}}$($t_{\rm N_{t}\rm N_{b}}$) denote the interlayer
hopping terms between top B (N) atom and bottom B (N) atoms. The hopping terms
for the interaction between top B or N atom with bottom N or B atoms are
denoted as $t_{B_{t}N_{b}}$. (b) Comparison of the low energy bands from DFT
and using the HTC model on BN/BN ($\theta=21.78^{\circ}$). The red (blue) real
line and orange (cyan) lines represent the relaxed (rigid) configuration from
DFT and HTC, respectively. The horizontal dashed line at $0$ eV is the Fermi
level. (c) Illustration of the local registry vector
${\bm{d}}=(d_{x}$,$d_{y})$ showing the sliding between the bottom and top
layer atoms in the same sublattice where the lattice is represented by the
rectangle unit cell $a_{\text{BN}}$($1\times\sqrt{3}$) where $a$ is the
lattice parameter. (d) The relative total energy landscape as a function of
$d_{x}$ and $d_{y}$ sliding values for the same rectangle unit cell relative
to AA and AA′ stackings. (e) and (f) Real and imaginary parts of $\rm B_{t}\rm
N_{b}$ interlayer matrix elements as a function of sliding based on AA
stacking. DFT, HTC and Diff represent the tunneling data calculated by DFT
Wannier representation, HTC model, and the difference between the DFT and the
HTC model data, respectively.

To improve on this generic model proposed originally for tBG, we reparametrize
$\gamma_{1}$ in Eq. (6) to match the specific interactions that exist between
B and N atoms. This is done by matching the DFT interlayer tunneling at the
$\mathbf{K}$ point [51, 42] defined as

$\displaystyle
H_{ss^{\prime}}(\bm{K}:\bm{d})=\sum_{j_{s^{\prime}}}t_{i_{s}j_{s^{\prime}}}^{\mathrm{inter}}\exp\left[\mathrm{i}\bm{K}\cdot\left(\bm{d}+\bm{r}_{i_{s}j_{s^{\prime}}}\right)\right]$
(8)

where $s$ and $s^{\prime}$ refer to the two different bottom-layer sublattices
$\alpha$ and $\beta$ and top-layer sublattices $\alpha^{\prime}$ and
$\beta^{\prime}$. Together with Eq. (4) this leads to the species-specific
terms [53] given by $\gamma_{1}=t_{\text{BB}^{\prime}}=0.831$ eV,
$\gamma_{1}=t_{\text{NN}^{\prime}}=0.6602$ eV, or
$\gamma_{1}=t_{\text{BN}^{\prime}}=t_{\text{NB}^{\prime}}=0.3989$ eV, see Fig.
2(a). Fixing the bottom layer sublattice by choosing one arbitrary $i_{s}$
site, we sum over all possible $j_{s^{\prime}}$ sites of sublattice
$s^{\prime}$ in the top layer until the sum is converged. As $\mathbf{K}$ is a
2D vector, only the in-plane components of the 3D distance vector
$\bm{r}_{i_{s}j_{s^{\prime}}}$ contribute to the phase factor. We can safely
ignore the contributions proportional to $\gamma_{0}$ in this species-specific
fitting procedure because the contribution from Eq. (5) is negligible in Eq.
(4) for distances of the order of the interlayer distance $c_{\perp}$ [41].
Using BN/BN 4-atom commensurate cells for zero twist angle, we illustrate the
DFT-tunneling maps, the fitted TB-tunneling maps, and their difference, for
the real and imaginary parts of $t_{\rm{B_{t}N_{b}}}$ as a function of
${\bm{d}}$ in Fig. 2(e) and 2(f). These displacement vectors
${\bm{d}}=(d_{x},d_{y})$ are illustrated in Fig. 2(c). The fitting of real and
imaginary parts for the other terms $t_{\rm{B_{t}B_{b}}}$,
$t_{\rm{N_{t}N_{b}}}$ and $t_{\rm{B_{t}N_{b}}}$ is illustrated in Fig. S.2 of
the Supplemental Material. We note that for all these terms, both for their
real and imaginary contributions, the maximum difference between DFT and TB
modeling is less than a few percent and we can safely assume that our TB model
offers an accurate parametrization of the main contribution to the interlayer
tunneling.

We further benchmark the accuracy of our models by comparing the low-energy
bands obtained from DFT and our HTC model for large twist angles. As shown in
Fig. 2(b) for a $21.78^{\circ}$ BN/BN system, we find excellent agreement
between the DFT and HTC $\pi$-bands especially at the $\mathbf{M}$ and
$\mathbf{K}$ points both for relaxed and rigid configurations. The small
differences near the $\mathbf{\Gamma}$ point are due to the $\sigma$-bands,
which are not included in our HTC model. To confirm the accuracy of our HTC
model for different stacking configurations, we compare the low-energy bands
obtained from DFT and HTC for six high-symmetry stackings. The results, shown
in Fig.S.3 of the Supplemental Material, reveal only minor differences at the
lowest CB and the highest VB at the $\mathbf{K}$ point. Based on these
benchmark calculations using commensurate cells, we can safely use the HTC
model for the rest of the manuscript to calculate the electronic band
structures and the local density of states in real space for the moiré
systems.

To investigate the impact of an external electric field in the moiré cells of
t2BN we introduce an additional correction proportional to the force
introduced by an electric field ($F_{z}=qE_{z}$) on a given B or N atom. Using
the electron charge $e$, the Bader charge [61] values of $q=0.82275|e|$ are
used for B and $-q$ for N atoms, respectively [42]. The $z$ atomic positions
relaxation under an external electric field using DFT (QE) and classical
energy minimization simulations (LAMMPS) show both an approximately linear
relationship and are in reasonable agreement with each other (see Fig. S.4(a)
of the Supplemental Material). Fig. S.4(b) in turn displays a schematic of the
atom position moving with the perpendicular electric field $E_{z}$ for AA and
AA′ stacking. respectively. The site-resolved electric potential can be
modeled through

$V_{i}=z_{i}\frac{E_{z}}{\epsilon_{r}}$ (9)

where $z_{i}$ is the z-position of atom $i$ within each relaxed structure. The
origin point of the z-axis is defined such that $\braket{z}=0$. We take
relative dielectric constant of BN/BN to be $\epsilon_{r}=2.7$ because it
gives the best fit of the band gap as a function of perpendicular electric
field in bilayer hBN. This value is comparable with the static perpendicular
dielectric constant value of $\epsilon_{r}=2.6$ reported in Ref. [62] and
somewhat smaller than $\epsilon_{r}=3.4$ in Ref. [63] for bilayer hBN.

##  III Results

In the following, we analyze the atomic and electronic structure of t2BN in
the absence and presence of a perpendicular electric field, together with the
local ferroelectric charge distribution.

###  III.1 Atomic structure

The atomic structure of twisted BN/BN systems undergoes relaxations depending
on the energetics of local stacking total energies at zero twist angle, in a
manner similar to tBG. In panel Fig. 2(d), we show the energy landscape for
zero twist angle BN/BN and BN/NB systems. The results reveal that for BN/BN
systems, the AB and BA stackings have the lowest energy, whereas for BN/NB
systems, the $\rm AA^{\prime}$ stacking is the most stable one due to the
attractive interaction between the top B(N) atoms and the bottom N(B) atoms
[64, 65].

Figure 3:  Various representations of lattice reconstruction effects,
schematics of charge transfer, and sliding-dependent gap estimates for the
BN/BN and BN/NB systems, for a twist angle equal to $1.08^{\circ}$. (a) and
(d) Stacking area maps of the different high-symmetry stacking regions
determined by the in-plane displacement (also, see Fig.S.5 of the Supplemental
Material) for the rigid configurations of BN/BN and BN/NB, respectively. (b)
and (e) Similar maps for the relaxed configurations of BN/BN and BN/NB. (c)
and (f) The in-plane displacement distance $d^{xy}_{ij}$ for the relaxed BN/BN
and BN/NB systems, where $i$ and $j$ indices refer to the closest pair of
atoms in different layers belonging to the same sub-lattice, where $\rm
B_{t}\rm B_{b}$ denote the top layer B atoms and bottom layer B atoms. (g) and
(h) The out-of-plane disp lacement distance $d^{z}_{ij}$ for relaxed BN/BN and
BN/NB configurations, respectively. (i) Schematic representation of the side
views of the six high symmetry stackings under an external electric field.
$\rm\bm{E}_{z}$ denotes the external electric field, $p$ denotes the dipole
polarization between B and N atoms. The vertical motion directions of B and N
atoms upon the application of $\rm\bm{E}_{z}$ are indicated by red and blue
arrows. (j) The band gap as a function of stacking relative to AA (top) and
AA′ (bottom) stackings where the sliding path corresponds to the y direction
in Fig. 1 (d). $\Delta$ denotes the variation in interlayer distance where
$\Delta=0$ corresponds to the equilibrium interlayer distance. Figure 4:  Low-
energy bands, band gaps, lattice-reconstructed stacking areas depending twist
angle for the BN/BN and BN/NB systems. (a) and (b) Low energy CB and VB for
various twist angles ranging from $3.14^{\circ}$ to $1.08^{\circ}$ in BN/BN.
(c) and (d) Similar plots as (a) and (b) but for BN/NB. The red-to-yellow
scale to denote the angle is used for the relaxed systems and the blue-to-
green scale is used for the rigid systems. (e) and (f) Band gap as a function
of twist angle for rigid (Rig) and relaxed (Rel) configurations in BN/BN and
BN/NB, respectively. (g) and (h) Twist angle dependence of the VB and CB for
rigid (Rig) and relaxed (Rel) moiré pattern in BN/BN and BN/NB. (i) The
stacking area percentage with twist angles of relaxed and rigid configurations
for BN/BN. (j) Similar plots for BN/NB. Figure 5:  Stacking regions and
electronic properties under external electric fields for the BN/BN and BN/NB
systems. (a) and (b) Percentage variation of stacking areas versus twist angle
for BN/BN (AA, AB, and BA stacking regions) and BN/NB (AA′, AB′, and BA′
stacking regions) systems under a perpendicular electric field ranging from
$-$0.5 V/nm to 0.5 V/nm. (c) Energy per atom dependence on sliding stacking
relative AA stacking under external electric fields ranging from $-0.5$ V/nm
to 0.5 V/nm. (d) Band gap of the BN/BN twisted system under various electric
fields and different twist angles. (e) and (f) Electronic polarization
distributions in rigid and relaxed $\theta=1.08^{\circ}$ BN/BN moiré systems.
(g) The maximum and electronic polarization simulated by the TB model as a
function of twist angle in the BN/BN system. Red and blue dots correspond to
the relaxed and rigid configurations, respectively, while the corresponding HA
is given with dashed lines. The individual first and third-order contributions
to this HA for the relaxed system are presented as continuous dotted lines.

We then assign a high symmetry local stacking to each atom in the moiré system
following the approach illustrated in Fig.S.5 of the Supplemental Material.
This allows us to transform continuous in-plane displacement ${d}^{xy}_{ij}$
maps, as seen in Fig. 3(c) and (f), to discrete maps in Fig. 3(a), (b), (d),
and (e), which divide the relaxed domain regions into AA($\rm AA^{\prime}$),
AB($\rm AB^{\prime}$)/BA($\rm BA^{\prime}$) stacking regions. We use these
discrete maps to have a simple and intuitive representation of the lattice
reconstruction.

Continuous colormaps of the out-of-plane displacement ${d}^{z}_{ij}$, for both
BN/BN and BN/NB systems, are shown in Figs. 3(g) and (h). We notice a
qualitative difference in the relaxation profile depending on the overall
alignment sign. For BN/BN, the size of the AA-stacking region shrinks in favor
of the other two local stackings, while the $\rm AA^{\prime}$ stacking in
BN/NB has a larger domain region size after relaxation. The size of these
different domain regions can be understood from the total energy landscape as
a function of stacking shown in Fig. 2(d).

Specifically, a more stable stacking gives rise to a larger domain size, while
the unstable stacking shrinks to a smaller one. In addition to the lattice
reconstructions observed without an electric field, we also find that a
vertical external electric field significantly impacts the system. Assuming
that B and N atoms move in opposite directions due to their positive and
negative effective charges, respectively, it makes sense to classify the
interlayer interactions into four different contributions $\rm N_{t}\rm
B_{b}$, $\rm B_{t}\rm N_{b}$, $\rm N_{b}\rm B_{t}$, and $\rm B_{b}\rm N_{t}$,
where $t$ and $b$ indicate the top and bottom layers, respectively. This
separation is illustrated in Fig. S.4 of the Supplemental Material. By
observing the evolution of the electric field-dependent interlayer distance
$d^{z}_{ij}$ in Fig. S.6 of the Supplemental Material for both BN/BN and BN/NB
systems, we can track the sizes of the six different local stacking regions of
the BN/BN and BN/NB systems that respond differently depending on the sign of
the electric fields, as illustrated in Fig. 3(i). In this schematic, we
illustrate the local interlayer charge polarization resulting from the overall
positive charge accumulation on B atoms and negative charge accumulation on N
atoms. The diagram demonstrates how the charge distribution within the
interlayer region is influenced by the contrasting charges for B and N atoms.

For BN/BN the layer charge polarization is achieved clearly for AB and BA
stacking regions where the interlayer dipole pointing along the $z$-direction
is denoted by $p$ in Fig. 3(i). The opposite sign of the dipoles for different
local stacking AB and BA domains implies asymmetric out-of-plane response
under opposite electric fields. In contrast, for BN/NB no local polarization
effect exists in the $\rm AA^{\prime}$, $\rm AB^{\prime}$, and $\rm
BA^{\prime}$ stacking configurations, which explains why the sign of the
electric field does not alter the local out-of-plane deformation nor the local
interlayer separation. This is shown more clearly in Fig. S.6(c) and (d) of
the Supplemental Material where we show the role of the electric field in
modifying the stacking regions for BN/BN and BN/NB configurations, examining a
range of electric fields spanning from $-0.5$ V/nm to $0.5$ V/nm. Due to the
reversed polarity of AB and BA stackings in BN/BN, as shown in Fig 3(i), their
respective sizes evolve in opposite ways for positive and negative electric
fields. Conversely, for BN/NB, the different local stackings have almost no
interlayer polarization, and result in nearly the same domain sizes and
stacking registry maps regardless of the applied electric field. In the next
sections, we relate these changes in the local stacking registry with the
behavior of the low-energy bands, including the changes in the band gap size
under an electric field.

###  III.2 Electronic structure of twisted bilayer BN

Nearly flat low energy bands without magic angles are expected near charge
neutrality in sufficiently small twist angle t2BN thanks to the presence of a
band gap that reduces the group velocity of the electrons near the band edges
[66, 67, 36, 27, 16]. Here, we revisit the electronic structure of this system
with systematically improved atomic and electronic structure models to explore
a wider range of twist angles, layer alignments, and perpendicular electric
fields. Specifically, as our HTC model distinguishes the hopping terms between
different B or N atomic species, we study the BN/BN and BN/NB alignments in
their relaxed and rigid configurations, in the presence and absence of an
electric field. We refer to Fig. S.6(a) and (b) of the Supplemental Material
for real-space local stacking distributions in different commensurate twist
angles for BN/BN and BN/NB. In Fig. 3(j) we show the local band gap for zero
twist systems for different interlayer sliding and distances. The vertical
distribution of the positively and negatively charged boron and nitrogen sites
can create a finite local interlayer dipole that leads to interlayer charge
differences that we explore later. As expected, from an atomic structure point
of view, the local stacking distributions undergo changes upon lattice
relaxations that tend to expand the local stacking areas that are
energetically more stable.

In particular, we investigate the impact of twist angle and lattice
reconstruction on the behavior of the low energy VB and CB. We present their
evolution as a function of twist angle in Fig. 4(a)-4(d) where we compare the
results of rigid and relaxed geometries. For rigid geometries, we use a fixed
interlayer separation of $c_{\perp}$ = 3.261 Å as in Eq. (6). We observe that
the band gaps decrease for twisted systems after relaxation as shown in Figs.
4(e)-4(f). As expected, the changes are more substantial for small twist
angles when the moiré pattern strain profiles become more important. In the
BN/BN case, we observe that the VBs shift to a lower energy region and the CBs
shift to a higher energy region after relaxation. In contrast, in the BN/NB
case, the CBs and VBs shift upwards in energy upon relaxation compared to the
rigid case.

To further understand the electronic properties in relation to the twist
angle, we present in Figs. 4(g) and 4(h) the variations of the bandwidths as a
function of the twist angle for both rigid and relaxed superlattices in our
two moiré patterns. The VB and CB that form the band gaps reach their maximum
flatness below $1.08^{\circ}$ and $1.5^{\circ}$ degree for BN/BN and BN/NB,
respectively. The reduced bandwidths of the VB and CB can be correlated with
an increase in the size of the domain regions that contribute to those bands.
Indeed, to check the wave function distribution on VBs and CBs, the electron
density distribution at high symmetry points of $\mathbf{M}$,
$\mathbf{\Gamma}$ and $\mathbf{K}$ are plotted in Fig. S.7 of the Supplemental
Material. We show that the states at those high symmetry points distribute
both at B and N atoms, near the AA stacking sites for BN/BN both for VB and CB
states, while for BN/NB, the valence band states accumulate mainly at the N
atoms near BA′ stacking and CB states mainly accumulate at the B atoms near
AB′ stacking [67]. Thus, this distribution of electron density, in conjunction
with the stacking area percentage allows us to understand the bandwidth
variation. In Fig. 4(i), the angle-dependence of stacking areas is presented
for the BN/BN configuration. This diagram illustrates the decrease in the AA
stacking region, which can be associated with the increased valence and
conduction bandwidths shown in Fig. 4(g). Conversely, in the BN/NB
configuration shown in Fig. 4(j), the reduction in the BA′ region is
associated with the increased valence bandwidth while the expansion of the AB′
region to the decreased conduction bandwidth.

In the following, we consider the effect of applying a perpendicular electric
field as a means to further control the structural and associated electronic
properties.

###  III.3 Electric field

The electric field dependence of the atomic structure for the BN/BN and BN/NB
systems is shown in Fig. 5(a) and (b) through the stacking area evolution in a
BN/BN system. The local stacking regions gradually shrink or expand with
varying twist angles. Using the conventions of stacking area definition
outlined in Fig. S.5 of the Supplemental Material, the AA stacking reduces
down to about 10% of the total area, while AB/BA stackings almost make up 45%
in BN/BN system without external perpendicular electric field. When an
external perpendicular electric field is introduced, the AA stacking region
remains mostly unaffected. On the contrary, the AB and BA regions show an
inverse trend for their relative ratios with respect to the total area.
Specifically, the BA stacking area expands to 60%, while the AB stacking area
contracts to 25% when subjected to a $0.5$ V/nm electric field, for small
twist angles of around $\sim 1^{\circ}$. The expansion and shrinking of the AB
and BA stacking area ratios depend on the direction of the applied electric
field. These changes in local domain sizes stem from the impact of the
electric field on the site potentials within specific stacking regions, as
illustrated in Fig. S.4(b) of the Supplemental Material.

Notably, B and N atoms move in opposite directions in response to an electric
field due to their opposite charge. This behavior can be understood by
examining the energy dependence of the different stacking geometries when a
vertical external electric field is applied, as shown in Fig. 5(c). The AB and
BA stacking configurations demonstrate opposite stability behaviors under the
influence of positive and negative electric fields. Contrary to the BN/BN
alignment we find a symmetric response to the electric field direction for the
three representative local stacking configurations. Likewise, the energy-
dependent sliding behavior of stacking can elucidate this relative
insensitivity to an electric field. Indeed, the AA′, AB′, and BA′
configurations demonstrate consistent stability under the influence of both
positive and negative electric fields, see Fig. S.4(c) of the Supplemental
Material.

Our results show that the electric field can modify the overall interface
interlayer ferroelectricity in t2BN by changing the size of the different
stacking domains that have opposite layer charge polarization. This is further
illustrated in Fig. S.8 of the Supplemental Material. The panels (d) in Fig. 5
show the corresponding electric field evolution of the band-gap size relaxed
BN/BN structures at different twist angles. The gap sizes decrease upon
applying an electric field.

We note that the dependence of electronic structure on perpendicular electric
fields has been studied for a variety of 2D heterostructures combining
graphene and hBN [68, 69, 70, 71], and for hBN bilayers [30, 72, 31] where
perpendicular electric fields play a role in modifying the local stacking.
Experimentally, it was shown that in-plane sliding can switch the interlayer
polarization in bilayer transition metal dichalcogenides and the specific
stacking domains have associated layer-resolved charge polarization
corresponding to those local stacking domains [72]. Our current calculations
on t2BN show qualitatively similar trends. The observed behavior stems from
the repulsive potential difference between B-B or N-N, which linearly
correlates with the change in the amplitude of the bandwidths. As a result,
the electric field mediated change of local stacking in hBN bilayers [30, 72,
31] can be an interesting pathway for engineering the electronic and optical
properties in layer-charge polarized moiré systems.

Figure 6: Moiré potential observed in an adjacent material arising from its
proximity with the polarized BN/BN moiré system. (a) Schematic representation
of the induced moiré potential in the adjacent material using a blue/red
colormap. The gray map depicts the electron density difference between the top
and bottom layers in the BN/BN moiré system. (b) Moiré period dependence of
the maximally induced potential for different vertical distances in relaxed
and rigid BN/BN system with the gray gradient dashed lines representing the
single HA (SHA) where S refers to using a constant electronic polarization of
$2.25$ pC/m. We represent $\Delta V=V_{\rm max}-V_{\rm min}$ the the maximum
potential modulation depths for different moire periods.

###  III.4 Interlayer polarization induced moiré potential

Following our discussion of the angle-dependence of the flat bands in our two
moiré systems, we now shift our focus to the interfacial charge polarization,
which has been observed in recent experiments [29, 30, 31]. Interlayer
electron density redistribution is expected in BN/BN moiré systems [36, 31,
30, 29] but not in the BN/NB system. Indeed, for BN/NB alignment the
interlayer charge polarization becomes zero at AA′, AB′ and BA′ stacking
region, which corresponds to zero intrinsic interlayer charge polarization as
illustrated in Fig. 3(i). Hence, in the following electron density analysis,
we focus only on the BN/BN system that has a significant interface domain
polarization. Inverse electronic polarization can be observed between two
stacking regions corresponding to the states located at either the AB or BA-
stacking regions. In order to obtain the electron density difference and
electronic polarization in the small twist angle regime, we fit the local
interlayer electron density difference $n(\bm{r})$ in BN/BN as a function of
twist angle $\theta$ for each layer through a third-order HA for relaxed
BN/BN, see Eqs. (18)-(20) in Ref. [57], using here the real-space coordinate
$\bm{r}=(x,y)$-dependence which reads, after simplification, as

$\displaystyle\Delta n({\bm{r}})$ $\displaystyle=$
$\displaystyle\sum_{j=1,3}2c_{j}\left[\sin(G_{j}y)\right.$ (10)
$\displaystyle-$
$\displaystyle\left.4\cos(\sqrt{3}G_{j}x/2)\sin(G_{j}y/2)\right]$

where $c_{1}$ and $c_{3}$ are the coefficients of the first-order and third-
order HA. Considering that only $c_{1}$ terms play a significant role in the
rigid BN/BN system, the first-order HA is used for the rigid BN/BN simulation.
In contrast, both the first-order and third-order HA are considered in the
relaxed BN/BN simulation due to their concurrent sizeable contributions (see
Supplemental Material in Fig.S.9). The carrier densities are given in
electrons per unit $\rm cm^{2}$. We use $G_{1}=4\pi/\sqrt{3}a_{m}$, and
$G_{3}=2G_{1}$, where $a_{m}$ is the real space periodicity of the moiré
pattern. The total interlayer electron density difference for top layer
coordinate $\bm{r}$, can be obtained by adding the contributions from B and N
atoms calculated separately using

$\Delta n(\bm{r})=\Delta n_{\rm B}({\bm{r}})+\Delta n_{\rm N}({\bm{r}})$ (11)

where the species-dependent electron density differences are calculated from
the electron density at each sublattice in the top and bottom layers (see Eq.
10) as

$\Delta n_{\rm B}({\bm{r}})=n_{\rm B_{t}}({\bm{r}})-n_{\rm B_{b}}({\bm{r}})$
(12)

and

$\Delta n_{\rm N}({\bm{r}})=n_{\rm N_{t}}({\bm{r}})-n_{\rm N_{b}}({\bm{r}}).$
(13)

This separation is more convenient for our calculations considering the
strongly distinct electronegativities of the B and N sites.

The $\theta$-dependence of the parameters $c_{1}(\theta)$ for interlayer
electron density of the rigid system and $c_{1}(\theta)$ and $c_{3}(\theta)$
for relaxed BN/BN systems are illustrated in Fig. S.10 of the Supplemental
Material and can be fitted using the exponential function

$c_{j}(\theta)=\alpha e^{\beta\theta}+\gamma e^{\delta\theta^{2}}\\\ $ (14)

whose individual fitting parameters $\alpha$, $\beta$, $\gamma$ and $\delta$
are listed in Table 3. We provide parameterization sets both for rigid and
relaxed geometries. Moreover, charge density at each sublattice in the top and
bottom layers for rigid and relaxed BN/BN also provided in Fig. S.11 and the
fitting coefficients are listed in Table I of the Supplemental Material.

Table 3:  The $\alpha$, $\beta$, $\gamma$, and $\delta$ fitting parameters in units of $cm^{-2}$ for $c_{\Delta_{\rm B/\rm N},1}$, $c_{\Delta_{\rm B/\rm N},3}$ in Eq. (14) to find $\Delta n_{\rm B}$ and $\Delta n_{\rm N}$ areal density differences between top and bottom B and N atoms using Eq. (10), for both rigid (Rig) and relaxed (Rel) BN/BN configurations. The $1$ and $3$ subscripts refer to the first and third-order coefficients. It is worth noting that the first-order HA proportional to $c_{1}$ cannot accurately fit small twist angles interlayer electron distributions in the relaxed model and we must include the third order $c_{3}$ term.  $\mathrm{Rig}$ | $\alpha$ | $\beta$ | $\gamma$ | $\delta$ |
---|---|---|---|---|---
${c}_{\rm\Delta B,1}$ | $0.00$ | - | $5.81\times 10^{11}$ | $-0.020$ |
${c}_{\rm\Delta N,1}$ | $0.00$ | - | $1.11\times 10^{12}$ | $-0.006$ |
$\rm Rel$ | $\alpha$ | $\beta$ | $\gamma$ | $\delta$ |
${c}_{\rm\Delta B,1}$ | $1.24\times 10^{12}$ | $-0.13$ | $-2.46\times 10^{11}$ | $-0.41$ |
${c}_{\rm\Delta N,1}$ | $1.93\times 10^{12}$ | $-0.11$ | $6.46\times 10^{11}$ | $-0.58$ |
${c}_{\rm\Delta B,3}$ | $1.73\times 10^{11}$ | $-0.31$ | $3.64\times 10^{12}$ | $0.49$ |
${c}_{\rm\Delta N,3}$ | $-3.57\times 10^{11}$ | $-0.03$ | $1.42\times 10^{12}$ | $-0.64$ |

The electronic polarization

$P(\textbf{r})=-e\Delta n(\textbf{r})c_{\perp}(\textbf{r})$ (15)

where $e=1.602\times 10^{-19}$ C, can be obtained from the local interlayer
electron density difference $\Delta n$ between top and bottom layers and the
interlayer distance $c_{\perp}$. Fig. 5(e) and (f) show TB-simulated
electronic polarization occurring at the interface for rigid and relaxed
$\theta=1.08^{\circ}$ BN/BN. The local stacking registry
${\bm{d}}=(d_{x},d_{y})$ dependent electronic polarization and interlayer
electron density difference $\Delta n({\bm{d}})$ for 4-atoms bilayer
commensurate cells is shown in Fig. S.12 for calculations based on TB, DFT,
and the first harmonic (HA) fitting for zero twist angle. We find a nearly
perfect agreement between TB and DFT, and the Harmonic approximation (HA) fit
to first-order shows a difference in the density smaller than 1% near density
maxima points. We find a good HA fitting of the electronic polarization for
finite twisted BN/BN systems in Fig. 5(g) where we show the maximum electronic
polarization for both relaxed and rigid systems as a function of twist angle,
where the lattice relaxation reveals a more pronounced variation in the
electronic polarization with twist, which is likely due to the larger AB and
BA stacking regions in the relaxed model. For relaxed and rigid BN/BN, the
electronic polarization can be converged to $2.6$ pC/m and $2.23$ pC/m for
large moiré patterns, respectively, where the rigid electronic polarization is
in good agreement with theoretical and experimental values ranging from $2.00$
to $2.25$ pC/m in Refs. [30, 73, 74] for AB-stacked bilayer h-BN. The relaxed
BN/BN reaches a maximal interlayer polarization of $3.5$ pC/m at about
$2.13^{\circ}$ ($\sim 7.5$ nm) due to the increasing importance of the third
order terms at small angles that accounts for lattice reconstruction effects
(purple dotted line in Fig. 5(g)).

This interlayer electronic polarization in the BN/BN moiré configuration leads
to induced moiré potentials that can influence the properties of adjacent
materials [35]. We thus calculate such induced electronic static potentials
using [73]

$V(\mathbf{r},z)\approx\sum_{j=1,3}-{\rm sgn}(z)\frac{e\Delta
n_{j}(\textbf{r})c_{\perp}(\textbf{r})}{2\varepsilon_{0}}e^{-G_{j}|z|}\\\ $
(16)

where $\varepsilon_{0}$ stands for the vacuum permittivity, and
$n_{j}(\mathbf{r})$ is the electron density difference at $\bm{r}$ for the
$j^{\text{th}}$ order of the harmonic approximation. $z$ represents the
vertical distance measured from the middle of the t2BN to the adjacent
material, as illustrated in Fig. 6(a). Fig. 6(b) shows the moiré period
dependence of the maximum potential difference value of $\Delta V=V_{\rm
max}-V_{\rm min}$ for vertical interlayer distances ranging from $1$ nm to
$20$ nm, in which $n_{j}$ corresponds to the electron density difference at
the AB stacking in each BN/BN moiré pattern. The polarization is determined
using a first-order HA ($j=1$) simulated electronic polarization based on the
electron density distribution obtained through Eq. (10) (indicated by the blue
to green gradient lines), yielding $240$ mV at a $1$ nm vertical distance. In
contrast, relaxed BN/BN, based on the first and third-order HA ($j=1,3$)
(depicted by red gradient lines), exhibits an increased induced moiré
potential, reaching $296$ mV. This value closely aligns with experimental
findings [31, 30, 29]. We further illustrate with the SHA curve in Fig. 6(b)
(represented by gray gradient dashed lines) that using a constant polarization
[73] of $2.25$ pC/m instead of the twist-angle-dependent values from Fig. 5(g)
is a reasonable approximation when considering the rigid configuration.

We note that the induced potential does not sizeably change anymore for moiré
patterns with periods larger than $600$ nm, i.e. twist angles smaller than
$\sim 0.02^{\circ}$. Furthermore, the simulated relaxed moiré potential
demonstrates near vertical distance independence for large moiré lengths, with
variations smaller than $10$ mV. This closely matches the distance dependence
of moiré potential for large moiré pattern estimated in experiments. Based on
the observations in this Section, we speculate that the moiré potential
arising from t2BN may significantly enhance the versatility of moiré
engineering, thus potentially expanding the functionalities of electronic and
photonic layered materials.

##  IV Conclusions

In this work, we have quantitatively obtained the local interlayer charge
polarization in t2BN that gives rise to moiré antiferroelectric domains when
local stacking changes in real space. This calculation is carried out through
a tight-binding model, informed by DFT calculations where we model separately
the intralayer contributions with up to six nearest neighbor hopping terms and
the interlayer tunneling is captured through a systematically improved
distance-dependent two-center approximation that distinguishes the different
atomic species. We have carried out our TB calculations on real-space lattice
rigid and relaxed geometries and have made quantitative estimates on the
effects that lattice relaxations have on the electronic structure of the
system, generally enhancing the gaps and widening the bandwidths. We have
further studied the atomic and electronic structures in the presence of an
external electric field that can be used to expand or shrink certain local
stacking regions. Finally, the interface electron density profiles of the
moiré patterns have been quantified through parametrization employing third-
order harmonic expansions, thus enhancing their adaptability for future
research endeavors. Based on these electron density profiles, we explore the
moiré potential induced in the adjacent material contacted by t2BN in the
BN/BN configurations, indicating potential avenues for manipulating the
properties of neighboring functional layers by leveraging the surface
potential of a t2BN substrate.

###### Acknowledgements.

This work was supported by the National Research Foundation of Korea (NRF)
through grant numbers NRF2020R1A5A1016518 (F.L.) and RS-2023-00249414 (N.L.).
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. D.L.
and J.J. also acknowledge support by the Korean Ministry of Land,
Infrastructure and Transport (MOLIT) from the Innovative Talent Education
Program for Smart Cities.

## References

  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Park  _et al._ [2022] J. M. Park, Y. Cao, L.-Q. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature Materials 21, 877 (2022).
  * Ledwith  _et al._ [2021] P. J. Ledwith, E. Khalaf, Z. Zhu, S. Carr, E. Kaxiras, and A. Vishwanath, arXiv:2111.11060 (2021).
  * Leconte  _et al._ [2022a] N. Leconte, Y. Park, J. An, A. Samudrala, and J. Jung, 2D Materials 9, 044002 (2022a).
  * 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).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Physical review letters 124, 097601 (2020).
  * 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).
  * Xie  _et al._ [2019] Y. Xie, B. Lian, B. Jäck, X. Liu, C.-L. Chiu, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 572, 101 (2019).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * Serlin  _et al._ [2020] M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Science 367, 900 (2020).
  * Wong  _et al._ [2020] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 582, 198 (2020).
  * Zondiner  _et al._ [2020] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern, _et al._ , Nature 582, 203 (2020).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, _et al._ , Nature 574, 653 (2019).
  * 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).
  * Javvaji  _et al._ [2020] S. Javvaji, J.-H. Sun, and J. Jung, Physical Review B 101, 125411 (2020).
  * Yang  _et al._ [2018] Q. Yang, M. Wu, and J. Li, The journal of physical chemistry letters 9, 7160 (2018).
  * Naik and Jain [2018] M. H. Naik and M. Jain, Physical review letters 121, 266401 (2018).
  * Pan  _et al._ [2020] H. Pan, F. Wu, and S. D. Sarma, Physical Review Research 2, 033087 (2020).
  * Shimazaki  _et al._ [2020] Y. Shimazaki, I. Schwartz, K. Watanabe, T. Taniguchi, M. Kroner, and A. Imamoğlu, Nature 580, 472 (2020).
  * An  _et al._ [2020] L. An, X. Cai, D. Pei, M. Huang, Z. Wu, Z. Zhou, J. Lin, Z. Ying, Z. Ye, X. Feng, _et al._ , Nanoscale horizons 5, 1309 (2020).
  * Regan  _et al._ [2020] E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, _et al._ , Nature 579, 359 (2020).
  * Zhang  _et al._ [2020] Y. Zhang, N. F. Yuan, and L. Fu, Physical Review B 102, 201115 (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).
  * Tang  _et al._ [2020] Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, _et al._ , Nature 579, 353 (2020).
  * Li  _et al._ [2021] H. Li, S. Li, M. H. Naik, J. Xie, X. Li, J. Wang, E. Regan, D. Wang, W. Zhao, S. Zhao, _et al._ , Nature materials 20, 945 (2021).
  * Xian  _et al._ [2019] L. Xian, D. M. Kennes, N. Tancogne-Dejean, M. Altarelli, and A. Rubio, Nano letters 19, 4934 (2019).
  * Yao  _et al._ [2021] K. Yao, N. R. Finney, J. Zhang, S. L. Moore, L. Xian, N. Tancogne-Dejean, F. Liu, J. Ardelean, X. Xu, D. Halbertal, _et al._ , Science Advances 7, eabe8691 (2021).
  * Vizner Stern  _et al._ [2021] M. Vizner Stern, Y. Waschitz, W. Cao, I. Nevo, K. Watanabe, T. Taniguchi, E. Sela, M. Urbakh, O. Hod, and M. Ben Shalom, Science 372, 1462 (2021).
  * Yasuda  _et al._ [2021] K. Yasuda, X. Wang, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Science 372, 1458 (2021).
  * Woods  _et al._ [2021] C. Woods, P. Ares, H. Nevison-Andrews, M. Holwill, R. Fabregas, F. Guinea, A. Geim, K. Novoselov, N. Walet, and L. Fumagalli, Nature communications 12, 347 (2021).
  * Rode  _et al._ [2017] J. C. Rode, D. Smirnov, C. Belke, H. Schmidt, and R. J. Haug, Annalen der Physik 529, 1700025 (2017).
  * Moore  _et al._ [2021] S. Moore, C. Ciccarino, D. Halbertal, L. McGilly, N. Finney, K. Yao, Y. Shao, G. Ni, A. Sternbach, E. Telford, _et al._ , Nature Communications 12, 5741 (2021).
  * Ni  _et al._ [2019] G. Ni, H. Wang, B.-Y. Jiang, L. Chen, Y. Du, Z. Sun, M. Goldflam, A. Frenzel, X. Xie, M. Fogler, _et al._ , Nature communications 10, 4360 (2019).
  * Kim  _et al._ [2024] D. S. Kim, R. C. Dominguez, R. Mayorga-Luna, D. Ye, J. Embley, T. Tan, Y. Ni, Z. Liu, M. Ford, F. Y. Gao, _et al._ , Nature materials 23, 65 (2024).
  * Walet and Guinea [2021] N. R. Walet and F. Guinea, Physical Review B 103, 125427 (2021).
  * Anđelković  _et al._ [2020] M. Anđelković, S. P. Milovanović, L. Covaci, and F. M. Peeters, Nano Letters 20, 979 (2020).
  * Plimpton [1995] S. Plimpton, Journal of computational physics 117, 1 (1995).
  * Plimpton and Thompson [2012] S. J. Plimpton and A. P. Thompson, MRS bulletin 37, 513 (2012).
  * Javvaji  _et al._ [2024] S. Javvaji, F. Li, and J. Jung, arXiv:2404.13041 [physics.hist-ph] (2024).
  * De Laissardiere  _et al._ [2012] G. T. De Laissardiere, D. Mayou, and L. Magaud, Physical Review B 86, 125413 (2012).
  * Leconte  _et al._ [2022b] N. Leconte, S. Javvaji, J. An, A. Samudrala, and J. Jung, Physical Review B 106, 115410 (2022b).
  * Note [1] We choose the rotation axis to be at the first atom in the small unit cell that defines the local stacking.
  * Los  _et al._ [2017] J. Los, J. Kroes, K. Albe, R. Gordillo, M. Katsnelson, and A. Fasolino, Physical Review B 96, 184108 (2017).
  * Wen  _et al._ [2018] M. Wen, S. Carr, S. Fang, E. Kaxiras, and E. B. Tadmor, Physical Review B 98, 235404 (2018).
  * Leconte  _et al._ [2017] N. Leconte, J. Jung, S. Lebègue, and T. Gould, Physical Review B 96, 195431 (2017).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, _et al._ , Journal of physics: Condensed matter 21, 395502 (2009).
  * Mostofi  _et al._ [2008] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer physics communications 178, 685 (2008).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Reviews of Modern Physics 84, 1419 (2012).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Physical Review B 23, 5048 (1981).
  * Jung  _et al._ [2014] J. Jung, A. Raoux, Z. Qiao, and A. H. MacDonald, Physical Review B 89, 205414 (2014).
  * Trambly de Laissardière  _et al._ [2010] G. Trambly de Laissardière, D. Mayou, and L. Magaud, Nano letters 10, 804 (2010).
  * Javvaji  _et al._ [2023] S. Javvaji, F. Li, and J. Jung, under preparation (2023).
  * Jung and MacDonald [2013] J. Jung and A. H. MacDonald, Physical Review B 87, 195450 (2013).
  * Jung and MacDonald [2014] J. Jung and A. H. MacDonald, Physical Review B 89, 035405 (2014).
  * Carr  _et al._ [2019] S. Carr, S. Fang, H. C. Po, A. Vishwanath, and E. Kaxiras, Physical Review Research 1, 033072 (2019).
  * Jung  _et al._ [2015] J. Jung, A. M. DaSilva, A. H. MacDonald, and S. Adam, Nature communications 6, 6308 (2015).
  * Pereira  _et al._ [2009] V. M. Pereira, A. C. Neto, and N. Peres, Physical Review B 80, 045401 (2009).
  * Lynch and Drickamer [1966] R. Lynch and H. Drickamer, The Journal of Chemical Physics 44, 181 (1966).
  * Slater and Koster [1954] J. C. Slater and G. F. Koster, Physical Review 94, 1498 (1954).
  * Henkelman  _et al._ [2006] G. Henkelman, A. Arnaldsson, and H. Jónsson, Computational Materials Science 36, 354 (2006).
  * Laturia and Vandenberghe [2017] A. Laturia and W. G. Vandenberghe, in  _2017 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD)_ (IEEE, 2017) pp. 337–340.
  * Laturia  _et al._ [2018] A. Laturia, M. L. Van de Put, and W. G. Vandenberghe, npj 2D Materials and Applications 2, 6 (2018).
  * Constantinescu  _et al._ [2013] G. Constantinescu, A. Kuc, and T. Heine, Physical review letters 111, 036104 (2013).
  * Marom  _et al._ [2010] N. Marom, J. Bernstein, J. Garel, A. Tkatchenko, E. Joselevich, L. Kronik, and O. Hod, Physical review letters 105, 046801 (2010).
  * Ochoa and Asenjo-Garcia [2020] H. Ochoa and A. Asenjo-Garcia, Physical Review Letters 125, 037402 (2020).
  * Zhao  _et al._ [2020] X.-J. Zhao, Y. Yang, D.-B. Zhang, and S.-H. Wei, Physical Review Letters 124, 086401 (2020).
  * Chen  _et al._ [2019] 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 (2019).
  * Long  _et al._ [2023] M. Long, Z. Zhan, P. A. Pantaleón, J. A. Silva-Guillén, F. Guinea, and S. Yuan, Physical Review B 107, 115140 (2023).
  * Chittari  _et al._ [2019] B. L. Chittari, G. Chen, Y. Zhang, F. Wang, and J. Jung, Physical review letters 122, 016401 (2019).
  * Zhu  _et al._ [2022] Z. Zhu, S. Carr, Q. Ma, and E. Kaxiras, Physical Review B 106, 205134 (2022).
  * Wang  _et al._ [2022] X. Wang, K. Yasuda, Y. Zhang, S. Liu, K. Watanabe, T. Taniguchi, J. Hone, L. Fu, and P. Jarillo-Herrero, Nature nanotechnology 17, 367 (2022).
  * Zhao  _et al._ [2021a] P. Zhao, C. Xiao, and W. Yao, npj 2D Materials and Applications 5, 38 (2021a).
  * Li and Wu [2017] L. Li and M. Wu, ACS nano 11, 6382 (2017).
  * Gonzalez and Stauber [2019] J. Gonzalez and T. Stauber, Physical review letters 122, 026801 (2019).
  * Xu and Balents [2018] C. Xu and L. Balents, Physical review letters 121, 087001 (2018).
  * Dos Santos  _et al._ [2007] J. L. Dos Santos, N. Peres, and A. C. Neto, Physical review letters 99, 256802 (2007).
  * Morell  _et al._ [2010] E. S. Morell, J. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Physical Review B 82, 121407 (2010).
  * Moon and Koshino [2012] P. Moon and M. Koshino, Physical Review B 85, 195458 (2012).
  * Wu  _et al._ [2018] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Physical review letters 121, 026402 (2018).
  * Kundu  _et al._ [2021] S. Kundu, I. Maity, R. Bajaj, H. Krishnamurthy, and M. Jain, arXiv:2108.01112 [physics.hist-ph] (2021).
  * Zheng  _et al._ [2020] Z. Zheng, Q. Ma, Z. Bi, S. de La Barrera, M.-H. Liu, N. Mao, Y. Zhang, N. Kiper, K. Watanabe, T. Taniguchi, _et al._ , Nature 588, 71 (2020).
  * Zhao  _et al._ [2021b] X.-J. Zhao, Y. Yang, D.-B. Zhang, and S.-H. Wei, Physical Review Materials 5, 014007 (2021b).
  * de la Barrera  _et al._ [2021] S. C. de la Barrera, Q. Cao, Y. Gao, Y. Gao, V. S. Bheemarasetty, J. Yan, D. G. Mandrus, W. Zhu, D. Xiao, and B. M. Hunt, Nature communications 12, 1 (2021).
  * Chiodini  _et al._ [2022] S. Chiodini, J. Kerfoot, G. Venturi, S. Mignuzzi, E. M. Alexeev, B. Teixeira Rosa, S. Tongay, T. Taniguchi, K. Watanabe, A. C. Ferrari, _et al._ , ACS nano 16, 7589 (2022).
  * Alden  _et al._ [2013] J. S. Alden, A. W. Tsen, P. Y. Huang, R. Hovden, L. Brown, J. Park, D. A. Muller, and P. L. McEuen, Proceedings of the National Academy of Sciences 110, 11256 (2013).
  * Van Wijk  _et al._ [2015] M. Van Wijk, A. Schuring, M. Katsnelson, and A. Fasolino, 2D Materials 2, 034010 (2015).
  * Kazmierczak  _et al._ [2021] N. P. Kazmierczak, M. Van Winkle, C. Ophus, K. C. Bustillo, S. Carr, H. G. Brown, J. Ciston, T. Taniguchi, K. Watanabe, and D. K. Bediako, Nature materials 20, 956 (2021).
  * Strand  _et al._ [2019] J. Strand, L. Larcher, and A. L. Shluger, Journal of Physics: Condensed Matter 32, 055706 (2019).
  * King-Smith and Vanderbilt [1993] R. King-Smith and D. Vanderbilt, Physical Review B 47, 1651 (1993).
  * Chanana and Waghmare [2019] A. Chanana and U. V. Waghmare, Physical review letters 123, 037601 (2019).