#  Quantum spin Hall effect from multi-scale band inversion in twisted bilayer
Bi2(Te1-xSex)3

Ikuma Tateishi<EMAIL_ADDRESS>RIKEN Center for Emergent Matter
Science, Wako, Saitama 351-0198, Japan  Motoaki Hirayama
<EMAIL_ADDRESS>RIKEN Center for Emergent Matter Science, Wako,
Saitama 351-0198, Japan  Department of Applied Physics, University of Tokyo,
Tokyo 113-8656, Japan  JST, PRESTO, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan

###### Abstract

Moiré materials have become one of the most active fields in material science
in recent years due to their high tunability, and their unique properties
emerge from the Moiré-scale structure modulation. Here, we propose twisted
bilayer Bi2(Te1-xSex)3 as a new Moiré material where the Moiré-scale
modulation induces a topological phase transition. We show, in twisted bilayer
Bi2(Te1-xSex)3, a topological insulator domain and a normal insulator domain
coexist in the Moiré lattice structure, and edge states on the domain boundary
make nearly flat bands that dominate the material properties. The edge states
further contribute to a Moiré-scale band inversion, resulting in Moiré-scale
topological states. There are corresponding Moiré-scale edge states and they
are so to speak “edge state from edge state”, which is a unique feature of
twisted bilayer Bi2(Te1-xSex)3. Our result not only proposes novel quantum
phases in twisted bilayer Bi2Te3-family, but also suggests the twisting of
stacking sensitive topological materials paves an avenue in the search for
novel quantum materials and devices.

##  I Introduction

The twisted van der Waals heterostructure materials, or Moiré materials, have
been studied very intensively in recent years as a platform for exploring
novel quantum phases [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20]. In those materials, Moiré superlattices are formed by the
lattice misalignment with a small twist angle, and the Moiré superlattices
produce flat electric bands and various strongly correlated phases. In
particular, the experimental reports on the magic-angle twisted bilayer
graphene have stimulated this field [1, 2]. They reported that the bilayer
graphene stacked with a twist angle $1.08^{\circ}$ (magic angle), which has
been known to have flat bands near the Fermi level, shows correlated
insulating phases and superconducting phases when the filling factor is tuned.
Because the behavior resembles the phase diagram of the high-temperature
cuprate superconductors, twisted bilayer graphene has attracted great
attention. The unique feature of the Moiré material is its high tunability.
The twist angle is a tunable parameter specific to Moiré materials.
Furthermore, because the system is two-dimensional (2D) and has a large Moiré
unit cell, the filling factor can be tuned easily and significantly. This high
tunability allows us to find various quantum phases in a single Moiré
material. Inspired by the twisted bilayer graphene, Moiré systems of some
other layered materials have been studied. In the twisted bilayer transition
metal dichalcogenides (TMD), nearly flat bands have been theoretically
predicted on the valence band edge, and an experimental signature of a
correlated insulator phase has been reported [21, 22, 23, 24, 25, 26, 27, 28,
29, 30, 31, 32]. Also for other materials, such as hexagonal boron nitride
(hBN), the existence of nearly flat bands has been suggested [33, 34].
Although the tunability of Moiré materials is remarkable, research of Moiré
materials has concentrated on the layered materials above. Hence, the physics
that describes their low-energy electronic states can be qualitatively
categorized into two groups, semimetallic one (graphene) and insulating ones
(TMD, hBN).

In this paper, we theoretically propose twisted bilayer Bi2(Te1-xSex)3 (Fig.1)
as a new Moiré material that hosts novel low-energy electronic states
described by a topological phase transition and corresponding topological edge
states. Three-dimensional bulk Bi2(Te1-xSex)3 is one of the van der Waals
heterostructure materials, and is well known as a typical strong topological
insulator [35, 36, 37, 38]. For the thin film Bi2(Te1-xSex)3 case, it has been
suggested that the topological invariant strongly depends on the number of
stacked layers [39]. Therefore, topological phase transitions are expected to
occur when the stacking order or interlayer distance is changed. Generally in
a Moiré material, the local stacking order and interlayer distance are
modulated by the lattice misalignment [20]. Combining the stacking modulation
and the stacking sensitive topological insulator, we propose a Moiré material
with mixed topological insulator domain and normal insulator domain.

##  II Effective model for Moiré materials

In this paper, to avoid confusion, an atomic-scale lattice structure in an
untwisted system are explicitly referred to as “atomic lattice” and a Moiré
superlattice structure in a twisted system is referred as “Moiré lattice”
(Fig.1). The words “atomic” and “Moiré” are used in the same way for other
terms, such as atomic(Moiré) Brillouin zone (BZ).

The Moiré electronic states are calculated with an effective model within a
small-angle approximation [40, 41].

Figure 1: Lattice structure of twisted bilayer Bi2(Te1-xSex)3. (a) Schematic
picture of the Moiré pattern in twisted bilayer Bi2(Te1-xSex)3. There are
three in-plane atomic sites, A-, B-, and C-site, in the untwisted
Bi2(Te1-xSex)3. The atomic lattices in the upper and lower layers are drawn
with red and blue markers, respectively. The solid rhombus represents the
Moiré unit cell. The three sampling points used in the effective model are
denoted with circles. (b) Schematic picture of the Moiré BZ. The red and blue
dashed hexagons are twisted atomic BZ of the upper and lower layers,
respectively. $\bm{G}_{0}\sim\bm{G}_{6}$ are the Moiré reciprocal lattice
vectors that we used in the effective model. (c) Top and horizontal views of
the atomic lattice structures of untwisted bilayer Bi2(Te1-xSex)3 for AA-
(left), AB- (middle), and AC-stacking (right), drawn by VESTA [42]. Each of
the upper and lower layers has five atoms in the atomic unit cell.

The Hamiltonian of the effective model is given as

$\begin{split}H&=\int
d\bm{k}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{\bm{G}_{l}}t_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}~{}c_{\alpha\sigma,\bm{k}}^{\dagger}c_{\beta^{\prime}\sigma^{\prime},\bm{k}-\bm{G}_{l}}^{\phantom{\dagger}},\end{split}$
(1)

where $\alpha\sigma$ and $\beta^{\prime}\sigma^{\prime}$ are orbital-spin
indices, and $c^{\dagger}$ ($c$) is a creation (annihilation) operator. The
sum $\sum_{\bm{G}_{l}}$ is taken over seven Moiré reciprocal lattice vectors
$\bm{G}_{0}\sim\bm{G}_{6}$ (Fig. 1(b)). The
$t_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}$
are determined to satisfy

$\begin{split}\sum_{\bm{G}_{l}}e^{i\bm{G}_{l}\cdot\bm{r}_{j}}~{}t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}-\bm{G}_{m},\bm{G}_{l}}=&~{}h_{\bm{r}_{j}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{m}/2),\end{split}$
(2)

where $h_{\bm{r}_{j}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})$ is
the matrix elements calculated from hopping parameters around position
$\bm{r}_{j}$ (See Appendix A, B, and C for more detail of the model
derivation). The local atomic lattice structure around $\bm{r}_{j}$ is
approximated by an untwisted lattice with a particular stacking order, and
thus $h_{\bm{r}_{j}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})$ and
the electronic states on it are also estimated by calculation in the untwisted
lattice. We take three three sampling points of $\bm{r}_{j}$, and interpolate
the intermediate region by the discrete Fourier transform. We call the three
stacking orders at the three sampling points AA-, AB-, and AC-stacking (Fig.
1(a),(c)). The AB-stacking is the most stable one and thus it is realized in
the 3D bulk Bi2(Te1-xSex)3. With the three sampling points approximation, the
explicit definition of
$t_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}$
are given as

$\begin{split}t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k},\bm{G}_{0}}=&\frac{1}{3}\left[h_{\bm{r}_{AA}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})+h_{\bm{r}_{AB}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})+h_{\bm{r}_{AC}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})\right],\\\
t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}=&\frac{1}{9}\left[h_{\bm{r}_{AA}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\right.\\\
&~{}~{}~{}~{}~{}+e^{-i\frac{2\pi}{3}}h_{\bm{r}_{AB}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\\\
&~{}~{}~{}~{}~{}\left.+e^{i\frac{2\pi}{3}}h_{\bm{r}_{AC}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\right]~{}~{}(l=1,3,5),\\\
t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}=&\frac{1}{9}\left[h_{\bm{r}_{AA}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\right.\\\
&~{}~{}~{}~{}~{}+e^{i\frac{2\pi}{3}}h_{\bm{r}_{AB}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\\\
&~{}~{}~{}~{}~{}\left.+e^{-i\frac{2\pi}{3}}h_{\bm{r}_{AC}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\right]~{}~{}(l=2,4,6).\end{split}$
(3)

$h_{\bm{r}_{AA}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})$,
$h_{\bm{r}_{AB}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})$, and
$h_{\bm{r}_{AC}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})$ are
matrix elements of Hamiltonians calculated in the AA-, AB-, and AC-stacking
untwisted bilayer Bi2(Te1-xSex)3, respectively.

Figure 2: Electronic band structure and Wilson loop spectra of the untwisted bilayer Bi2Te3 and Bi2Se3 for each stacking orders. (a), (b), and (c) are for AA-, AB-, and AC-stacking bilayer Bi2Te3, and (d), (e), and (f) are for those of bilayer Bi2Se3, respectively. The magenta and yellow dots in the band structure figures represent the projected weight on the $p_{z}$ orbitals of Bi and Te/Se atoms, respectively. Only AB-stacking bilayer Bi2Te3 (b) is a topological insulator, and the others are normal insulators. Bi2Te3 |  AA |  AB |  AC
---|---|---|---
in-plane lat. const. $a$ (Å) | 4.40 | 4.41 | 4.39
interlayer distance $D$ (Å) | 3.18 | 2.73 | 3.98
Bi2Se3 | AA | AB | AC
in-plane lat. const. $a$ (Å) | 4.15 | 4.16 | 4.15
interlayer distance $D$ (Å) | 3.05 | 2.75 | 3.78
Table 1: Optimized in-plane lattice constant $a$ and interlayer distance $D$
for bilayer Bi2Te3 and Bi2Se3 with three stacking orders. The AA-, AB-, and
AC-stacking are defined in Fig.1.

$h_{\bm{r}_{AA}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})$,
$h_{\bm{r}_{AB}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})$, and
$h_{\bm{r}_{AC}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})$ are
obtained by the first-principles calculation. All first-principles
calculations are implemented in the Vienna ab initio simulation package (VASP)
[43]. We use the projector augmented wave (PAW) potential sets recommended by
VASP and set the kinetic-energy cutoff to 500 eV. For each of the three
stacking orders, we optimize the lattice structure and obtain electronic band
structures. The lattice optimizations are performed in the strongly
constrained and appropriately normed (SCAN) meta-generalized gradient
approximation for short- and intermediate-range interactions with the long-
range vdW interaction (rVV10) with the spin-orbit interaction [44]. We
calculate the electronic band structure in the B3LYP with the VWN3-correlation
[45]. We construct Wannier functions for the Bi and Te/Se $p$ orbitals with
WANNIER90 package [46]. We use $8\times 8\times 1$ $k$-mesh for the lattice
optimization and $9\times 9\times 1$ $k$-mesh for the electronic calculation
and the Wannier function. The obtained band structures reproduce well the
results of angle-resolved photoemission spectroscopy (ARPES) [47, 48] and the
GW approximation [49]. The Fermi level is set in the averaged Hamiltonian of
three stacking orders
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k},\bm{G}_{0}}$ by
the filling factor, i. e. the middle of the 36th and 37th band in the $\Gamma$
point of the 60 bands (2 layers $\times$ 5 atoms $\times$ $p$ orbitals
$\times$ spin). The obtained Wannier functions and matrix elements are also
used in the calculations of the Wilson loop spectra with the WannierTools
package [50]. We also calculate the Sb2Te3 under the same condition (see
Appendix D).

Figure 3: Moiré band structure of twisted bilayer Bi2(Te1-xSex)3. (a) Moiré
band dispersion of Bi2(Te0.85Se0.15)3 with a twist angle
$\theta=0.50^{\circ}$. The top (VB1), second (VB2), third (VB3), and fourth
(VB4) valence band pairs are shown as orange, red, green, and purple lines,
respectively. The symbols of the high-symmetry points are defined in Fig.1(b).
(b)-(e) Real space plot of the wave functions (upper layer, lower Bi, $p_{z}$
orbital, spin up component, as shown in (f)) of VB1-4 in the $\Gamma_{M}$
point. The brightness and color indicate the absolute value (normalized to
have the maximum value of 1) and phase of the wave function as shown in the
right of (c). The Moiré unit cell is shown as a white dashed rhombus. The
angular momentum $l$ is also shown in each figure. (f) Bi atom that is focused
in the plots.  Figure 4: Moiré wave function and local electronic state. The
rhombus is a wave function plot of VB2 in the Moiré unit cell (same as
Fig.3(c)). The lower plot shows the real space dependence of the absolute
value of the wave function of VB2 (cyan line, right vertical axis) and the
bandgap calculated in the interpolated untwisted Hamiltonian (violet line,
left vertical axis). The horizontal axis is the real space position along the
AA-AB-AC-AA-stacking line, which is the longer diagonal of the Moiré unit cell
as shown with a white arrow in the top figure. The negative bandgap means that
the bands are inverted. The two gapless points around the AB-stacking area are
the domain boundary, and the topological insulator domain is shown as a green-
shaded range. Figure 5: Twist angle dependence of Moiré band dispersion of
Bi2(Te0.85Se0.15)3. (a), (b), and (c) are the cases with a twist angle
$\theta=1.00^{\circ}$, $1.50^{\circ}$, and $2.00^{\circ}$, respectively. VB1-4
are shown with the same color as Fig.3(a). (d) Schematic picture to explain
the twist angle dependence of the flatness and gap of the edge-state-
originated nearly flat bands.

##  III Bilayer Bi2(Te1-xSex)3

First, we show the result of calculations on the untwisted atomic lattices of
bilayer Bi2Te3 and Bi2Se3 for the three stacking orders. The atomic positions
are shown in Fig. 1(c). All these untwisted bilayer lattices belong to the
layer group #72 (or the space group #166 with infinitely long c-axis). The
optimized in-plane lattice constant $a$ and interlayer distance $D$, which is
defined as the vertical distance between the two Te/Se atoms in the twist face
(See Fig.1(c), right), are listed in Table 1. We neglect the stacking order
dependence in the in-plane lattice constant in each material and use an
averaged value in the effective model calculations. For Te/Se doping, the in-
plane lattice constant is linearly interpolated. In both materials, the AB-
stacking has the smallest interlayer distance and the AC-stacking has the
largest. The obtained electronic band structures are shown in Figs.2(a)-(f)
for both of Bi2Te3 and Bi2Se3, where the Fermi level is determined by the
filling factor. The magenta and yellow dots represent the projected weight on
the $p_{z}$ orbitals of Bi and (Te,Se) atoms, respectively. In these
materials, the overlap of the $p_{z}$ orbitals contributes to a topological
phase transition. Therefore, the smaller the interlayer distance, the more
likely it is to be a topological insulator. To evaluate their topological
invariants, we make Wannier functions for them and calculate the Wilson loop
spectra as shown in the right panel of each figure of Figs.2(a)-(f). We can
see only AB-stacking Bi2Te3 is a topological insulator and all of the others
are normal insulators. These results indicate that a twisted bilayer Bi2Te3
system has a topological insulator domain around the AB-stacking region, and a
normal insulator domain in the other region. Although Bi2Se3 is a normal
insulator in the three stacking orders, we consider Se doping in Bi2Te3 by the
linear interpolation to tune the system parameters. However, note that the
topological non-triviality of the AB-stacking Bi2Te3 plays an essential role
even in doped cases.

##  IV Moiré band and quantized edge state in twisted bilayer Bi2(Te1-xSex)3

Next, we show the Moiré band dispersion of twisted bilayer Bi2(Te1-xSex)3. Due
to the twisting, the inversion symmetry is broken and twisted bilayer
Bi2(Te1-xSex)3 belongs to the layer group #67 (or the space group #149 with
infinitely long c-axis). The in-plane $C_{2}$ axis exists along the AA-AB-AC-
stacking line (see Fig. 1(a)), which corresponds to the $\Gamma_{M}$-MM line
in the reciprocal space. The symbols of the high symmetry points in the Moiré
BZ are defined as Fig. 1(b). Figure 3(a) shows the band dispersion of the
twisted bilayer Bi2(Te0.85Se0.15)3 with a twist angle $\theta=0.50^{\circ}$.
Here, to obtain a clear domain boundary, the amount of Se ($x$) is determined
so that the gap in the AA- and AB-stacking would be roughly the same (See
Appendix D for the detail of Se-doping dependence of the gap and Moiré band
dispersion). All bands are doubly degenerate at time-reversal invariant
momenta (TRIM), the $\Gamma_{M}$ and MM points. Because of the absence of the
inversion symmetry, the Kramers degeneracy in the untwisted bilayer splits at
a general momentum. The split is easy to see in the conduction bands above
$0.020$ eV in the KM point, while it is too small to see in the valence bands.
It is worth noting that there are nearly flat bands around the Fermi level.
For the nearly flat valence bands (from the top to fourth valence band pairs
VB1, VB2, VB3, and VB4 in the Fig. 3(a)), real space plots of the wave
functions (upper layer, lower Bi, $p_{z}$ orbital, spin up component. See
Fig.3(f)) at the $\Gamma_{M}$ point are shown in Fig.3(b)-(e). The brightness
and color indicate the absolute value (normalized to have the maximum value of
1) and phase of the wave function, respectively (as shown in the right of
Fig.3(c)). The wave function has a ring-shaped density distribution, clearly
indicating that these nearly flat bands originate from the edge state
corresponding to the topological insulator domain around the AB-stacking area.
To compare with the wave functions of the flat bands, we calculate the real
space dependence of the bandgap between the valence top and conduction bottom
bands in the $\Gamma$ point $\Delta E(\Gamma)$ in the interpolated untwisted
Hamiltonian obtained by the discrete Fourier transform in Eq. (2) (See
Appendix D for more detail). In Fig.4, the real space dependence of the band
bap (violet line) and the absolute value of the wave function of VB2 (cyan
line) along the AA-AB-AC-AA-stacking line (the longer diagonal of the Moiré
unit cell shown as a white arrow) are shown. The negative bandgap means that
the bands are inverted. The gapless points are the domain boundary and the
topological insulator domain is shown as a green-shaded range. We can confirm
that the wave function has a large amplitude around the domain boundary. Next,
we focus on the difference between the wave functions of VB1-4 (Fig.3(b)-(e)).
They show ring-shaped density distributions in common, but their phase
structures are different from each other. We can see a lower energy band has a
larger angular momentum $l$, which is calculated as a winding number of the
phase along the ring. Note that because these states are coupled with spin,
the spin-down component has one different $l$ and the total angular momentum
is a half-integer (See Appendix E for more detail). The phase structure
indicates that the nearly flat bands are formed by quantization of topological
edge states due to the finite size effect of the domain boundary. The flatness
of these bands and the gap size between them depend on the twist angle.
Figures 5(a)-(c) show the twist angle dependence of the Moiré band dispersion.
It can be seen that as the twist angle increases, the nearly flat bands get
more dispersive and the energy gaps between them get larger. This tendency is
due to the fact that the length of the Moiré reciprocal lattice vector, i.e.
the width of the Moiré BZ, is proportional to $\theta$. When the original
Dirac cone on the domain boundary is folded in Moiré BZ and quantized, folding
with large Moiré BZ results in the folding of bands with widely different
energies (Fig.5(d)). Further, folding with large Moiré BZ relatively keeps the
original band dispersion, and thus the obtained Moiré bands tend to have
larger bandwidth (See Appendix E for more detail). With these points of view,
the twist angle dependence of the flatness and energy gap of the edge-state-
originated nearly flat bands is roughly explained by the size of Moiré BZ, but
in more detail, the effect of the Rashba splitting and hybridization with
other bands must be taken into account. Finally, we mention the conduction
bands. In Bi2(Te1-xSex)3 case, since the conduction bands are more dispersive
than the valence bands in the band structures of untwisted lattices, the Moiré
conduction bands are also more likely to have dispersion in the higher energy
region. For the same reason, a bandgap between Moiré conduction bands tends to
be larger than that of Moiré valence bands.

Figure 6: Quantum spin Hall effect in twisted bilayer Bi2(Te1-xSex)3. (a)
Moiré band dispersion of Bi2(Te0.92Se0.08)3 with a twist angle
$\theta=1.55^{\circ}$. The top (VB1), second (VB2), third (VB3), fourth (VB4),
eighth (VB8), and ninth (VB9) valence band pairs are shown as orange, red,
green, purple, blue, and yellow lines, respectively. (b)-(g) Wilson loop
spectra calculated on VB1-4, VB8, and VB9. While (b)-(e) VB1-4 are
topologically trivial, (f) VB8 and (g) VB9 are topologically nontrivial. (h)
Magnified view of the VB8 and VB9. In this figure. the numerical error due to
the first-principles calculation is corrected to recover the exact time-
reversal symmetry. (i),(j) Real space plot of the wave functions (upper layer,
lower Bi, $p_{z}$ orbital, spin up component) of VB8 and VB9 in the
$\Gamma_{M}$ point, respectively. VB8 has a ring-shaped density distribution
around the AB-stacking area, but VB9 has no longer ring-shaped one but
localized density distribution around the AA-stacking area. Figure 7:
Topological edge state in twisted bilayer Bi2(Te1-xSex)3. (a) Moiré-scale edge
BZ that is used in the edge calculation. (b) Edge band spectrum around VB8 and
VB9. The symbols of the high-symmetry points are defined in (a). Dirac cones
of the Moiré-scale helical edge state can be seen appearing around the
$\bar{\Gamma}_{M}$ point. (c) Schematic picture of “edge state from edge
state”. The red and blue arrows correspond to spin currents having opposite
spin directions to each other. The helical edge state obtained in Moiré bands
(Moiré-scale helical edge state) is running along the edge of the Moiré-scale
lattice system. The Moiré-scale edge state is made from the ring-shaped edge
state on the domain boundary around the AB-stacking area (atomic-scale helical
edge state).

##  V Moiré-scale topological band and edge state

We move on to a discussion of the emergence of topological bands in the Moiré
bands of twisted bilayer Bi2(Te1-xSex)3. As we discussed above, the formation
of the edge-state-originated nearly flat bands are explained by a simple
theory and they are ordered with angular momentum $l$. Therefore, we could not
find a case to have a band inversion between the edge-state-originated bands.
However, in the region farther from the Fermi level, a bulk-state-originated
band appears and it can hybridize with an edge-state-originated to form
topologically nontrivial bands. For example, the band dispersion of
Bi2(Te0.92Se0.08)3 with a twist angle $\theta=1.55^{\circ}$ is shown in
Fig.6(a). Here, the amount of Se ($x$) is determined so that the gap in the
AA-stacking area would be smaller than that of the AB-stacking area to obtain
the bulk and edge states in different areas. To evaluate the topological
invariants for each band, the Wilson loop spectra are calculated and given in
Fig.6(b)-(g). While the nearly flat band around the Fermi level (VB1-4,
Fig6(b)-(e)) are evaluated as topologically trivial, the bands around $-0.065$
eV (the eighth and ninth valence band pairs VB8 and VB9, Fig. 6(f),(g)) are
evaluated as topologically nontrivial. Figure 6(h) is a magnified view of the
VB8 and VB9. Here, we correct the numerical error due to the first-principles
calculation to recover the exact time-reversal symmetry, because the Moiré
system generally deals with a small energy scale and thus the numerical errors
get more noticeable. It can be seen that the two bands are gapped and thus if
the Fermi level is tuned and placed between them the system becomes a
topological insulator protected by the time-reversal symmetry. The real space
plots of the wave functions in the $\Gamma_{M}$ point of the two bands are
shown in Fig.6(i),(j). The wave function of VB8 (Fig.6(i)) has a ring-shaped
density distribution as well as the nearly flat bands around the Fermi level.
On the other hand, VB9 (Fig.6(j)) does not have a ring-shaped one, but has a
density distribution localized around the AA-stacking area (the corner of the
Moiré unit cell). Considering that the AA-stacking Bi2Te3 has the smallest gap
among all stacking orders (Fig.2), this state is understood as a state
originated from the bulk state of the AA-stacking area. Due to a symmetry
restriction of the hybridization, the simple bulk-state-originated band can
hybridize with a band with a small $l$. That is why the VB8 has $l=1$ (as we
explained, the spin-down component has one different $l$). To obtain edge
states corresponding to these topological bands, we make Wannier functions for
these bands and calculate the edge band spectra with Green’s function method
(technical details are described in Appendix.C.). The Moiré-scale edge BZ that
is used in the edge calculation is shown in Fig.7(a). The obtained edge band
spectrum is shown in Fig.7(b). As in a typical topological insulator, a Dirac-
cone of the helical edge state can be seen appearing around the
$\bar{\Gamma}_{M}$ point (See Appendix E for a symmetry restriction on the
Dirac-cone). These topological bands are obtained in the Moiré band
dispersion, and thus the corresponding helical edge state is running along the
edge of the Moiré-scale lattice system (Fig.7(c)). As we have seen, the Moiré-
scale edge state is partially made from the ring-shaped edge state on the
domain boundary around the AB-stacking area (let us say atomic-scale helical
edge state). Therefore, we can say the Moiré edge state is “edge state from
edge state” (Fig.7(c)), and this is a novel and unique phenomenon that is
observed in twisted bilayer Bi2(Te1-xSex)3. These topological properties also
have twist angle dependence. In the small-angle limit, more edge-state-
originated nearly flat bands appear around the Fermi level due to the band
folding with the small Moiré BZ, as seen in Fig.5. Therefore, as twist angle
is decreased, topological phase transitions occur on VB8 and VB9 and
eventually they become topologically trivial bands. The topological phase
transition goes through a gap-closing at general momenta as described in a
general theory of topological phase transitions in 2D systems [51, 52] (See
Appendix E and F for more detail). From this perspective, when the twist angle
gets smaller, topological bands tend to appear in lower (or higher) energy
bands.

The edge-state-originated nearly flat bands, Moiré-scale edge states, and
other properties obtained in twisted bilayer Bi2(Te1-xSex)3 are also
reproduced in a more simplified model that we propose, twisted Bernevig-
Hughes-Zhang model (See Appendix G). This fact not only allows us to easily
explore the properties of these systems, but also indicates that our essential
strategy can be applied to other topological materials.

##  VI Discussion

In conclusion, we have theoretically studied the electronic structure of
twisted bilayer Bi2(Te1-xSex)3. It is revealed that twisted bilayer
Bi2(Te1-xSex)3 has a topological insulator domain and a normal insulator
domain in the Moiré unit cell due to the stacking modulation of the Moiré
superlattice structure. We have obtained a Moiré-scale band inversion and
corresponding edge states that are made from the atomic-scale edge state on
the domain boundary, thus it can be called “edge state from edge state”.  With
these results, we propose twisted bilayer Bi2(Te1-xSex)3 as a new topological
Moiré material that hosts novel low-energy states and unique multi-scale band
inversion. Not only as a proposal of new material, but this proposal also
provides a new platform to study topological phases both in atomic- and Moiré-
scale. The atomic-scale helical edge states are expected to be used as an
ideal platform to realize the Chalker-Coddington network model [53, 54, 55,
56, 57] without an external field. The Moiré-scale topological states also
propose a realization of Moiré topological spin Hall materials that do not
require external fields. Further, because the wave function of the Moiré-scale
helical edge state has characteristic density distribution and far larger real
space size than that of helical edge states in previous topological insulator
crystals, it allows us to observe the nontriviality of wave function with a
real space observation such as scanning tunneling microscope (STM). While
those phenomena proposed by some previous studies require a strong external
field [58, 59, 60, 61, 62, 63, 64] or a particular symmetry-broken valley
setup [24], the role is replaced by the intrinsic strong spin-orbit coupling
in twisted bilayer Bi2(Te1-xSex)3. Moreover, twisted bilayer Bi2(Te1-xSex)3
can be a new platform to investigate correlated phases. As shown in the Moiré
band dispersion, large Rashba splits are predicted around KM points in our
calculation. It indicates that the ratio between bandwidth and the size of the
Rashba split can be tuned by changing the twist angle.

Although we have studied a specific material twisted bilayer Bi2(Te1-xSex)3
here, the strategy underlying this study — making a Moiré material with a
stacking sensitive topological material — is quite versatile. There is an
abundance of stacking-sensitive topological materials, and it should be
possible to design Moiré materials in which multiple topological phase domains
coexist with this strategy. The combination of a novel quantum system with
several different topological phase domains and the high tunability of the
Moiré materials provides an avenue in the search for new quantum materials and
devices.

## Acknowledgements

We are grateful to Akira Furusaki for fruitful discussion. This work was
supported by JST CREST (Grants No. JPMJCR19T2). M.H. was supported by PRESTO,
JST (JPMJPR21Q6) and JSPS KAKENHI Grants No. 20K14390.

##  Appendix A Moiré pattern

To make a 2D Moiré superlattice structure, there are two fundamental methods,
lattice constant mismatch and twisting identical lattices. General Moiré
superlattice systems are made by combining the two methods. In this paper, we
focus on the cases of twisting identical lattices.

###  A.1 General cases

Generally, given atomic lattice vectors $\bm{a}_{1}$, $\bm{a}_{2}$, and twist
angle $\theta$, the Moiré lattice vectors $\bm{L}_{1}$ and $\bm{L}_{2}$
(Fig.8(a)) satisfy

$\begin{split}2\sin\frac{\theta}{2}\bm{e}_{z}\times\bm{L}_{1}&=-\bm{a}_{2},\\\
2\sin\frac{\theta}{2}\bm{L}_{2}\times\bm{e}_{z}&=\bm{a}_{1},\end{split}$ (4)

within the small angle approximation (Fig.8(b)), where $\bm{e}_{z}$ is a unit
vector along the $z$ axis. By solving these equations, the Moire lattice
vectors are obtained as

$\begin{split}\bm{L}_{1}&=\bm{a}_{2}\times\frac{1}{2\sin\frac{\theta}{2}}\bm{e}_{z}\approx\bm{a}_{2}\times\frac{1}{\theta}\bm{e}_{z},\\\
\bm{L}_{2}&=\frac{1}{2\sin\frac{\theta}{2}}\bm{e}_{z}\times\bm{a}_{1}\approx\frac{1}{\theta}\bm{e}_{z}\times\bm{a}_{1}.\\\
\end{split}$ (5)

It is also proved that $\bm{a}_{l}$ and $\bm{L}_{m}$ satisfies

$\begin{split}\bm{a}_{l}\cdot\bm{L}_{m}=\delta_{lm}\frac{|\bm{a}_{1}||\bm{a}_{2}|}{2\sin\frac{\theta}{2}}\sin\gamma\approx\delta_{lm}\frac{|\bm{a}_{1}||\bm{a}_{2}|}{\theta}\sin\gamma,\end{split}$
(6)

where $\gamma$ is the angle between $\bm{a}_{1}$ and $\bm{a}_{2}$.

This relation indicates that the Moiré lattice (Moiré reciprocal lattice) of
the twisted system is similar to the atomic reciprocal lattice (atomic
lattice) in the original untwisted system.

Figure 8: (a) Atomic lattice vectors $\bm{a}_{1},\bm{a}_{2}$ and Moiré lattice
venctors $\bm{L}_{1},\bm{L}_{2}$. (b) Relation between the atomic lattice
vector, the Moiré lattice vector, and twist angle $\theta$.

###  A.2 Trigonal lattice

For a trigonal lattice system, atomic lattice vectors are given with the
lattice constant $a_{0}$ as

$\begin{split}\bm{a}_{1}&=a_{0}\left(1,0\right),\\\
\bm{a}_{2}&=a_{0}\left(-1/2,\sqrt{3}/2\right).\end{split}$ (7)

Moiré lattice vectors are given as

$\begin{split}\bm{L}_{1}&=L\left(\sqrt{3}/2,1/2\right),\\\
\bm{L}_{2}&=L\left(0,1\right),\end{split}$ (8)

where the Moiré lattice constant $L$ is written as

$L=\frac{a_{0}}{2\sin\frac{\theta}{2}}.$ (9)

##  Appendix B Effective model of Moiré superlattice systems

In this section, we derive a model of a Moiré superlattice system. In the
following, the atomic lattice constant is set to $a_{0}=1$. In subsection B.1,
we describe the general and exact part of the model derivation that is
independent of the detail of the target Moiré system. The description in this
subsection is consistent with previous studies [40, 41]. In subsection B.2, we
introduce the small angle approximation, which is generally used in a Moiré
system. In this subsection, we improve the approximation method of the
previous studies to make them more efficient. In subsection B.3, we describe
the system dependent part of the model derivation. In this subsection, we make
a correction to make the model satisfy symmetry restrictions of the twisted
bilayer Bi2(Te1-xSex)3 system.

###  B.1 General formula of model derivation

The exact Hamiltonian of a Moiré superlattice system is given by the following
with a microscopic picture.

$\begin{split}H=&\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{pq}t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)\\\
&\times
c^{\dagger}_{\alpha\sigma}\left(\bm{r}_{j}+\bm{d}_{\alpha}\left(\bm{r}_{j}\right)\right)c^{\phantom{\dagger}}_{\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}+\bm{d}_{\beta^{\prime}}\left(\bm{r}_{j}\right)+\bm{a}_{pq}\right).\end{split}$
(10)

$\begin{array}[]{ccl}\bm{r}_{j}&:&\mathrm{atomic~{}unit~{}cell~{}position}\\\
\alpha\sigma,\beta^{\prime}\sigma^{\prime}&:&\mathrm{orbital,~{}spin~{}index}\\\
\bm{d}_{\alpha,\beta^{\prime}}\left(\bm{r}_{j}\right)&:&\mathrm{site~{}position~{}of~{}}\alpha,\beta^{\prime}\mathrm{~{}in~{}the~{}unit~{}cell~{}}\bm{r}_{j}\\\
pq&:&\mathrm{hopping~{}lattice~{}index~{}}(\bm{a}_{pq})\\\
c^{\dagger},c&:&\mathrm{creation/annihilation~{}operator}\\\ \end{array}$ (11)

Here, $\bm{r}_{j}$ is a position of an atomic unit cell, and
$\alpha\sigma,\beta^{\prime}\sigma^{\prime}$ are indices to specify orbital
and spin. Now $\bm{r}_{j}$ is defined as the corner of the atomic unit cell,
and thus we define $\bm{d}_{\alpha}(\bm{r}_{j})$ as the site position of the
orbital $\alpha$ measured from the corner of the atomic unit cell
$\bm{r}_{j}$. Due to the effect of twisting, $\bm{d}_{\alpha}(\bm{r}_{j})$
generally depends on the atomic unit cell position $\bm{r}_{j}$.
$c^{\dagger}_{\alpha\sigma}(\bm{r}_{j}+\bm{d}_{\alpha}(\bm{r}_{j}))$ is a
creation operator of a fermion on the orbital $\alpha$ spin $\sigma$ in the
atomic unit cell $\bm{r}_{j}$, and $c$ is an annihilation operator.
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)$
is a hopping parameter which is defined as a hopping between “orbital $\alpha$
spin $\sigma$ in the atomic unit cell $\bm{r}_{j}$” and “orbital
$\beta^{\prime}$ spin $\sigma^{\prime}$ in the atomic unit cell
$\bm{r}_{j}+\bm{a}_{pq}$” (Fig.9). The $\bm{a}_{pq}$ are defined with integers
$p,q$ as

$\bm{a}_{pq}=p\bm{a}^{\prime}_{1}+q\bm{a}^{\prime}_{2},$ (12)

where $\bm{a}^{\prime}_{1}$, $\bm{a}^{\prime}_{2}$ are atomic lattice vectors
defined differently in intralayer and interlayer cases. The definitions of
$\bm{a}^{\prime}_{1}$ and $\bm{a}^{\prime}_{2}$ are

$\begin{array}[]{ccl}\bm{a}^{\prime}_{1}=C_{\theta/2}\bm{a}_{1}&\bm{a}^{\prime}_{2}=C_{\theta/2}\bm{a}_{2}&\mathrm{(intralayer,~{}upper~{}layer)}\\\
\bm{a}^{\prime}_{1}=C_{-\theta/2}\bm{a}_{1}&\bm{a}^{\prime}_{2}=C_{-\theta/2}\bm{a}_{2}&\mathrm{(intralayer,~{}lower~{}layer)}\\\
\bm{a}^{\prime}_{1}=\bm{a}_{1}&\bm{a}^{\prime}_{2}=\bm{a}_{2}&\mathrm{(interlayer)}\\\
\end{array}.$ (13)

We assume that long range hoppings are negligible and $\sum_{pq}$ is a finite
sum. There are two points to note about $\bm{a}_{pq}$ and $\sum_{pq}$. First,
in the interlayer case, $(\bm{a}_{1},\bm{a}_{2})$ does not match either the
atomic lattice period of the upper layer or that of the lower layer. As a
result, it can be possible that there is no $\beta^{\prime}$ site (or there
are two $\beta^{\prime}$ sites) in the atomic unit cell
$\bm{r}_{j}+\bm{a}_{pq}$. However, because the small angle limit is
interpreted as a continuous limit, those inconvenient cases are neglected.
Second, because $\bm{d}_{\alpha}(\bm{r}_{j})$ depends on the position in the
Moiré unit cell $\bm{r}_{j}$, a site can disappear on one end of the summation
range of $\sum_{pq}$ and appear on the other end when $\bm{r}_{j}$ is changed.
This will essentially make a “discontinuity” in the Hamiltonian. However, as
long as the range of $\sum_{pq}$ is large enough, the disappearing and
appearing hoppings have small absolute values and thus the discontinuity is
negligible.

Figure 9: Definition of the hopping
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)$.
$\bm{r}_{j}$ is a position that indicates atomic unit cell index. The
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)$
is a hopping between orbital $\alpha$ spin $\sigma$ in the atomic unit cell
$\bm{r}_{j}$ (red dashed cell) and prbital $\beta^{\prime}$ spin
$\sigma^{\prime}$ in the atomic unit cell $\bm{r}_{j}+\bm{a}_{pq}$ (black
dashed cell).

We perform a Fourier transform of Eq. (10). We assume that the wannier
functions are independent of position in the Moiré unit cell, and thus the
annihilation operator is transformed as

$c_{\alpha\sigma}(\bm{r}_{j})=\int
d\bm{k}~{}e^{i\bm{k}\cdot\bm{r}_{j}}c_{\alpha\sigma,\bm{k}}.$ (14)

The Fourier transform of Eq. (10) is obtained as

$\begin{split}H=&\int
d\bm{k}d\bm{k}^{\prime}\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\\\
&\left[\sum_{pq}t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)e^{-i\bm{k}\cdot\bm{d}_{\alpha}\left(\bm{r}_{j}\right)}e^{i\bm{k}^{\prime}\cdot\left(\bm{d}_{\beta^{\prime}}\left(\bm{r}_{j}\right)+\bm{a}_{pq}\right)}\right]\\\
&\times
e^{i\left(-\bm{k}^{\prime}+\bm{k}\right)\bm{r}_{j}}~{}c_{\alpha\sigma,\bm{k}}^{\dagger}c_{\beta^{\prime}\sigma^{\prime},\bm{k}^{\prime}}^{\phantom{\dagger}}.\end{split}$
(15)

We focus on the bracket part. The $\bm{d}_{\alpha}(\bm{r}_{j})$ and
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)$
are Moiré periodic. Therefore, the bracket part is Moiré periodic and can be
decomposed into components of the Moiré reciprocal lattice vectors
$\bm{G}_{l}$ as

$\begin{split}&\left[\sum_{pq}t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)e^{-i\bm{k}\cdot\bm{d}_{\alpha}\left(\bm{r}_{j}\right)}e^{i\bm{k}^{\prime}\cdot\left(\bm{d}_{\beta^{\prime}}\left(\bm{r}_{j}\right)+\bm{a}_{pq}\right)}\right]\\\
&=\sum_{\bm{G}_{l}}e^{i\bm{G}_{l}\cdot\bm{r}_{j}}~{}t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}},\end{split}$
(16)

where
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$
is the Fourier transform of the bracket part. Note that the
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$
is not the Fourier transform of
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)$
because the bracket part includes a sum on $p,q$ and phase factors.

By using Eq. (16), the Hamiltonian is writetn as

$\begin{split}H&=\int
d\bm{k}d\bm{k}^{\prime}\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{\bm{G}_{l}}~{}t_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}~{}e^{i\left(-\bm{k}+\bm{k}^{\prime}+\bm{G}_{l}\right)\cdot\bm{r}_{j}}\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times
c_{\alpha\sigma,\bm{k}}^{\dagger}c_{\beta^{\prime}\sigma^{\prime},\bm{k}^{\prime}}^{\phantom{\dagger}}\\\
&=\int
d\bm{k}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{\bm{G}_{l}}t_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}~{}c_{\alpha\sigma,\bm{k}}^{\dagger}c_{\beta^{\prime}\sigma^{\prime},\bm{k}-\bm{G}_{l}}^{\phantom{\dagger}}.\end{split}$
(17)

Note that $\bm{k}^{\prime}$ is restricted to $\bm{k}-\bm{G}_{l}$.

So far, the model derivation is exact and general. However, in the Eq. (16),
all
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)$
in the Moiré unit cell are needed to obtain exact value of
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$.
Because it is usually unrealistic, we use a small angle approximations to
obtain
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$,
as described in the following.

###  B.2 Small angle approximation

In a small twist angle case,
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$
can be approximately obtained as following.

First, we obtain
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)$
at finite number of sampling points $\\{\bm{r}_{j}\\}$ in the Moiré unit cell.
In a small twist angle case, the local atomic lattice structures around each
sampling points can be approximated by untwisted atomic lattices with proper
layer stacking orders. We perform first-principles calculation for each
untwisted structure and use the results as approximate value of
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)$
at the sampling points. With the approximated
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)$,
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$
are obtained by the discrete Fourier transform Eq. (16).

Figure 10: Examples of parameter sampling and $\bm{G}_{l}$ adopted in the
discrete Fourier transform in a trigonal lattice system. The left figure shows
the Moiré unit cell (red), sampling points, and sampling lattice vectors (blue
arrows) in the real space. The right figure shows the Moiré BZ (red dashed),
reciprocal vectors of the sampling lattice (blue arrows), and $\bm{G}_{l}$
adopted in the discrete Fourier transform (white and colored dots). For
$\bm{G}_{l}$ shown with the same colored dot,
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$
have the same value. Figure 11: Examples of parameter sampling and
$\bm{G}_{l}$ adopted in the discrete Fourier transform in a square lattice
system. The left figure shows the Moiré unit cell (red), sampling points, and
sampling lattice vectors (blue arrows) in the real space. The right figure
shows the Moiré BZ (red dashed), reciprocal vectors of the sampling lattice
(blue arrows), and $\bm{G}_{l}$ adopted in the discrete Fourier transform
(white and colored dots). For $\bm{G}_{l}$ shown with the same colored dot,
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$
have the same value.

Generally, the discrete Fourier transform guarantees that the parameters in
the sampling points are reproduced, but those in the intermediate area are not
necessarily estimated in a physically reasonable way. To obtain a reasonable
estimation for the whole Moiré unit cell by discrete Fourier transform, we
need to adopt appropriate $\bm{G}_{l}$ and sampling points must be determined
to capture the “feature values”, for example, the maximum and minimum values
of
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)$.
Therefore, the sampling points should be determined from the atomic positions
in the untwisted atomic lattice. Here, we explain general way to decide the
sampling points and how to choose $\bm{G}_{l}$ in detail, and then we show
some examples in Figs.10 and 11. As the first step, we decide untwisted
bilayer structures with interlayer in-plane shift $\\{\Delta\bm{r}_{j}\\}$
that should be sampled. For the convenience of the discrete Fourier transform,
the list $\\{\Delta\bm{r}_{j}\\}$ should be taken as a mesh in the atomic unit
cell. As explained above, the $\\{\Delta\bm{r}_{j}\\}$ should include the
cases where feature values are obtained. Typically, the
$\\{\Delta\bm{r}_{j}\\}$ include the cases where the atoms on the upper layer
and the lower layer are closest or farthest. We translate the
$\\{\Delta\bm{r}_{j}\\}$ to a position in the Moiré unit cell
$\\{\bm{r}_{j}\\}$ by using the equation
$\Delta\bm{r}_{j}=\bm{r}_{j}\times\theta\bm{e}_{z}$ (easily obtained from
Fig.8(b)). As a result, a sampling lattice (green crosses and blue arrows in
left figures in Figs.10 and 11) is obtained in the Moiré unit cell and a
sampling BZ is also defined in the reciprocal space (blue hexagons in the
right figures). The number of sampling points (the number of
$\Delta\bm{r}_{j}$, originally) decides a number of $\bm{G}_{l}$ that is taken
into account in Eq. (16) and Eq. (17). Moiré reciprocal lattice points
included in the sampling BZ (white dots) are adopted as $\bm{G}_{l}$.
Additionally, Moiré reciprocal lattice points on the boundary of the sampling
BZ (Cyan, magenta, and yellow dots) are also adopted as $\bm{G}_{l}$, but we
assume that
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$
have same value in the points that are connected by the reciprocal vectors of
the sampling lattice (blue arrows in the right figures). $\bm{G}_{l}$ outside
of the sampling BZ are neglected, i. e. components with frequencies higher
than the sampling points are neglected. By using this method, the number of
the sampling points and $\bm{G}_{l}$ are equal, and thus
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$
are uniquely determined by the discrete Fourier transform Eq. (16). Examples
of the way to choose $\bm{G}_{l}$ are shown in Figs.10 and 11 (for trigonal
and square cases, respectively). In each row, the left figure is the Moiré
lattice in the real space, and the right figure is the Moiré reciprocal
lattice in the reciprocal space. In the left figure, the Moiré unit cell (red)
and sampling points (green crosses) are shown. Blue arrows are basic lattice
vectors of the sampling lattice. In the right figure, the Moiré BZ (red
dashed) and the sampling BZ (blue dashed), which is defined by the sampling
lattice vectors, are shown. Moiré reciprocal lattice points that is adopted as
$\bm{G}_{l}$ are shown with dots. If it is included in the sampling BZ, the
dot is white. If it is on the boundary of the sampling BZ, the dot is colored
(cyan, magenta, or yellow) and the points that is connected by the reciprocal
vectors of the sampling lattice (blue arrows) have the same color.
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$
on those $\bm{G}_{l}$ have the same value. For example, in the case of the
three points sampling in Fig.10,

$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{1}}=t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{3}}=t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{5}}.$
(18)

Note that the number of the sampling points in the left figure equals to the
number of $\bm{G}_{l}$ in the right figure (dots with same color on the
boundary should be counted as one) in each column.

Note that depending on the symmetry of the system and the center of momentum
of the cutoff (See Sec.C.2), hopping parameters in some stacking orders can be
fixed to be 0 [65].

###  B.3 Specific formula for Bi2(Te1-xSex)3

In this subsection, we apply the method given in the previous subsection to
Bi2(Te1-xSex)3. In the case of Bi2(Te1-xSex)3, the in-plane atomic positions
are only three in the atomic unit cell ($C_{3}$ symmetry centers). Therefore,
the simplest choice is to calculate three untwisted stacking structures (the
three points sampling in Fig.10). Here, we define the three stacking orders as
AA, AB, and AC stacking. In this case, we take components of seven
$\bm{G}_{l}$,

$\bm{G}_{l}=\left\\{\bm{G}_{0},\bm{G}_{1},\bm{G}_{2},\bm{G}_{3},\bm{G}_{4},\bm{G}_{5},\bm{G}_{6}\right\\},$
(19)

and assume equivalences between them as

$\begin{split}t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{1}}&=t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{3}}=t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{5}},\\\
t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{2}}&=t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{4}}=t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{6}}.\end{split}$
(20)

Including a completely independent component

$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{0}},$
(21)

the degree of freedom is three, which is the same as the number of the
sampling points. Therefore, using the relation Eq. (20), we can determine the
seven
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}^{\prime},\bm{G}_{l}}$
from the three sampling points with the discrete Fourier transform.

The hopping parameters in the three sampling points,
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{AA}\right)$,
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{AB}\right)$,
and
$t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{AC}\right)$,
are calculated with the first-principle calculation of untwisted structures.
However, when the hopping parameters of the twisted structure are estimated
from untwisted structures, the Hamiltonian Eq. (17) will break the Hermitian
and $C_{3z}$ symmetry slightly. Therefore, we add a correction factor in Eq.
(16).

The restrictions by each symmetry are given as following. The detail of the
derivation of the restrictions is given in later subsection. The Hermitian
symmetry restriction is

$\left.t_{\bm{k}+\bm{G}_{l},\bm{k},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\right.^{*}=t_{\bm{k},\bm{k}+\bm{G}_{l},-\bm{G}_{l}}^{\beta^{\prime}\sigma^{\prime},\alpha\sigma}$
(22)

and the $C_{3z}$ symmetry restriction is

$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}=t^{\widetilde{\alpha\sigma},\widetilde{\beta^{\prime}\sigma^{\prime}}}_{C_{3}\bm{k},C_{3}\bm{k}-C_{3}\bm{G}_{l},C_{3}\bm{G}_{l}},$
(23)

where the orbital-spin index $\widetilde{\alpha\sigma}$ means the $C_{3z}$
rotated orbital-spin. The twisted structure also have time-reversal symmetry
and its restriction is

$(-1)^{1-\delta_{\sigma\sigma^{\prime}}}\left.t_{-\bm{k},-\bm{k}-\bm{G}_{l},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\right.^{*}=t_{\bm{k},\bm{k}+\bm{G}_{l},-\bm{G}_{l}}^{\alpha\bar{\sigma},\beta^{\prime}\bar{\sigma}^{\prime}},$
(24)

where $\bar{\sigma}$ means the opposite spin, but this equation is satisfied
when the three untwisted structures have the time-reversal symmetry.

To satisfy these equation, we add a correction factor
$e^{i\frac{\bm{G}_{m}}{2}\cdot\left(\bm{d}_{\alpha}(\bm{r}_{j})+\bm{d}_{\beta^{\prime}}(\bm{r}_{j})+\bm{a}_{pq}\right)}$
in the Eq. (16) as

$\begin{split}&\sum_{\bm{G}_{l}}e^{i\bm{G}_{l}\cdot\bm{r}_{j}}~{}t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}-\bm{G}_{m},\bm{G}_{l}}\\\
&=\left[\sum_{pq}t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)e^{-i\bm{k}\cdot\bm{d}_{\alpha}\left(\bm{r}_{j}\right)}e^{i\left(\bm{k}-\bm{G}_{m}\right)\cdot\left(\bm{d}_{\beta^{\prime}}\left(\bm{r}_{j}\right)+\bm{a}_{pq}\right)}\right.\\\
&\left.\phantom{~{}=\sum_{pq}t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)}\times
e^{i\frac{\bm{G}_{m}}{2}\cdot\left(\bm{d}_{\alpha}(\bm{r}_{j})+\bm{d}_{\beta^{\prime}}(\bm{r}_{j})+\bm{a}_{pq}\right)}\right]\\\
&=\left[\sum_{pq}t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)e^{i\left(\bm{k}-\frac{\bm{G}_{m}}{2}\right)\cdot\left(-\bm{d}_{\alpha}(\bm{r}_{j})+\bm{d}_{\beta^{\prime}}\left(\bm{r}_{j}\right)+\bm{a}_{pq}\right)}\right]\\\
&=~{}h_{\bm{r}_{j}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{m}/2).\end{split}$
(25)

Because the correction factor tends to 1 as the twist angle get small, this
correction is reasonable. Note that the right-hand-side
$h_{\bm{r}_{j}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{m}/2)$
is basically the matrix element of the Hamiltonian with momentum
$\bm{k}-\bm{G}_{m}/2$ given by the parameters of untwisted structure around
$\bm{r}_{j}$. The only difference is the momentum twisting due to the
different definition of $\bm{a}^{\prime}_{1,2}$ in Eq. (LABEL:eq.Model.a'def).
The explicit definition is given as, for
$\alpha\sigma,\beta^{\prime}\sigma^{\prime}\in$ Upper layer,

$\begin{split}&h_{\bm{r}_{j}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{m}/2)\\\
&=\sum_{pq}t_{\mathrm{DFT},pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{r}_{j})~{}e^{iC_{-\theta/2}\left(\bm{k}-\bm{G}_{m}/2\right)\cdot\left(-\bm{d}_{\alpha}(\bm{r}_{j})+\bm{d}_{\beta^{\prime}}(\bm{r}_{j})+\bm{a}_{pq}\right)},\end{split}$
(26)

for $\alpha\sigma,\beta^{\prime}\sigma^{\prime}\in$ Lower layer,

$\begin{split}&h_{\bm{r}_{j}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{m}/2)\\\
&=\sum_{pq}t_{\mathrm{DFT},pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{r}_{j})~{}e^{iC_{\theta/2}\left(\bm{k}-\bm{G}_{m}/2\right)\cdot\left(-\bm{d}_{\alpha}(\bm{r}_{j})+\bm{d}_{\beta^{\prime}}(\bm{r}_{j})+\bm{a}_{pq}\right)},\end{split}$
(27)

for $\alpha\sigma\in$ Upper Layer and $\beta^{\prime}\sigma^{\prime}\in$ Lower
layer,

$\begin{split}&h_{\bm{r}_{j}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{m}/2)\\\
&=\sum_{pq}t_{\mathrm{DFT},pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{r}_{j})~{}e^{i\left(\bm{k}-\bm{G}_{m}/2\right)\cdot\left(-\bm{d}_{\alpha}(\bm{r}_{j})+\bm{d}_{\beta^{\prime}}(\bm{r}_{j})+\bm{a}_{pq}\right)},\\\
&h_{\bm{r}_{j}}^{\beta^{\prime}\sigma^{\prime},\alpha\sigma}(\bm{k}-\bm{G}_{m}/2)={h_{\bm{r}_{j}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}}(\bm{k}-\bm{G}_{m}/2)^{\dagger},\end{split}$
(28)

where
$t_{\mathrm{DFT},pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{r}_{j})$
is the hopping parameter of the untwisted structure around $\bm{r}_{j}$
calculated by the first-principles calculation.

The discrete Fourier transform is performed as

$\begin{split}t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k},\bm{G}_{0}}=&\frac{1}{3}\left[h_{\bm{r}_{AA}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})+h_{\bm{r}_{AB}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})+h_{\bm{r}_{AC}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})\right],\\\
t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}=&\frac{1}{9}\left[h_{\bm{r}_{AA}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\right.\\\
&~{}~{}~{}~{}+e^{-i\frac{2\pi}{3}}h_{\bm{r}_{AB}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\\\
&~{}~{}~{}~{}\left.+e^{i\frac{2\pi}{3}}h_{\bm{r}_{AC}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\right]~{}~{}(l=1,3,5),\\\
t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}=&\frac{1}{9}\left[h_{\bm{r}_{AA}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\right.\\\
&~{}~{}~{}~{}+e^{i\frac{2\pi}{3}}h_{\bm{r}_{AB}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\\\
&~{}~{}~{}~{}\left.+e^{-i\frac{2\pi}{3}}h_{\bm{r}_{AC}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k}-\bm{G}_{l}/2)\right]~{}~{}(l=2,4,6).\end{split}$
(29)

With these definitions, Eqs. (20), (22), (23) and (24) are satisfied. The
$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k},\bm{G}_{0}}$
term can be understood as an averaged term, and the others as Moiré modulation
terms.

####  B.3.1 Gauge selection

In Eqs.(26)(27)(28), the phase factors are calculated with the explicit
difference of the orbital positions. There is another standard way to
calculate the phase factor, where a phase of the basic translations is given
only to inter-cell hopping terms. In an untwisted system, these two are
connected by a gauge transform, and thus have no physical difference. However,
in a Moiré system, they are physically different in a strict sense due to the
interpolation of the Hamiltonian. Although the difference is negligible in the
small angle limit with a small cutoff momentum, the former gives a better
interpolation for the large $\bm{k}$ region. Therefore, the former definition
is recommended, as we used.

####  B.3.2 Derivation of symmetry restrictions

We start from the Hamiltonian Eq.(10).

$\begin{split}H=&\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{pq}t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)~{}c^{\dagger}_{\alpha\sigma}\left(\bm{r}_{j}+\bm{d}_{\alpha}\left(\bm{r}_{j}\right)\right)c^{\phantom{\dagger}}_{\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}+\bm{d}_{\beta^{\prime}}\left(\bm{r}_{j}\right)+\bm{a}_{pq}\right)\end{split}$
(30)

  * •

Hermitian symmetry

The Hermite conjugate of the Hamiltonian is given as

$\begin{split}H^{\dagger}=&\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{pq}\left.t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\right.^{*}\left(\bm{r}_{i}\right)~{}c^{\dagger}_{\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}+\bm{d}_{\beta^{\prime}}\left(\bm{r}_{j}\right)+\bm{a}_{pq}\right)c^{\phantom{\dagger}}_{\alpha\sigma}\left(\bm{r}_{j}+\bm{d}_{\alpha}\left(\bm{r}_{j}\right)\right)\\\
=&\int
d\bm{k}d\bm{k}^{\prime}\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left[\sum_{pq}t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{r}_{j})e^{-i\bm{k}^{\prime}\cdot\bm{d}_{\alpha}(\bm{r}_{j})}e^{i\bm{k}\cdot\left(\bm{d}_{\beta^{\prime}}(\bm{r}_{j})+\bm{a}_{pq}\right)}\right]^{*}e^{i(-\bm{k}+\bm{k}^{\prime})\cdot\bm{r}_{j}}~{}c^{\dagger}_{\beta^{\prime}\sigma^{\prime},\bm{k}}c^{\phantom{\dagger}}_{\alpha\sigma,\bm{k}^{\prime}}\\\
=&\int
d\bm{k}d\bm{k}^{\prime}\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{\bm{G}_{l}}e^{-i\bm{G}_{l}\cdot\bm{r}_{j}}\left.t_{\bm{k}^{\prime},\bm{k},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\right.^{*}e^{i(-\bm{k}+\bm{k}^{\prime})\cdot\bm{r}_{j}}~{}c^{\dagger}_{\beta^{\prime}\sigma^{\prime},\bm{k}}c^{\phantom{\dagger}}_{\alpha\sigma,\bm{k}^{\prime}}\\\
=&\int
d\bm{k}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{\bm{G}_{l}}\left.t_{\bm{k}+\bm{G}_{l},\bm{k},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\right.^{*}c^{\dagger}_{\beta^{\prime}\sigma^{\prime},\bm{k}}c^{\phantom{\dagger}}_{\alpha\sigma,\bm{k}+\bm{G}_{l}}.\end{split}$
(31)

The conditions for this to be Hermite (Eq. (22)) is

$\left.t_{\bm{k}+\bm{G}_{l},\bm{k},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\right.^{*}=t_{\bm{k},\bm{k}+\bm{G}_{l},-\bm{G}_{l}}^{\beta^{\prime}\sigma^{\prime},\alpha\sigma}.$
(32)

When the Hamiltonian of the untwisted structures are Hermite,

$\left.h_{\bm{r}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})\right.^{*}=h_{\bm{r}}^{\beta^{\prime}\sigma^{\prime},\alpha\sigma}(\bm{k})$
(33)

is satisfied and the Eq. (22) is satisfied by using Eq. (29).

  * •

$C_{3z}$ symmetry

The $C_{3}$ rotated Hamiltonian is

$\begin{split}C_{3}HC_{3}^{\dagger}=&\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{pq}t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left(\bm{r}_{j}\right)~{}c^{\dagger}_{\widetilde{\alpha\sigma}}\left(C_{3}\bm{r}_{j}+C_{3}\bm{d}_{\alpha}(\bm{r}_{j})\right)c^{\phantom{\dagger}}_{\widetilde{\beta^{\prime}\sigma^{\prime}}}\left(C_{3}\bm{r}_{j}+C_{3}\bm{a}_{pq}+C_{3}\bm{d}_{\beta^{\prime}}(\bm{r}_{j})\right)\\\
&=\int
d\bm{k}d\bm{k^{\prime}}\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left[\sum_{pq}t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{r}_{j})e^{-iC_{3}^{-1}\bm{k}\cdot\bm{d}_{\alpha}(\bm{r}_{j})}e^{iC_{3}^{-1}\bm{k}^{\prime}\cdot\left(\bm{d}_{\beta^{\prime}}(\bm{r}_{j})+\bm{a}_{pq}\right)}\right]e^{iC_{3}^{-1}(-\bm{k}+\bm{k}^{\prime})\cdot\bm{r}_{j}}c^{\dagger}_{\widetilde{\alpha\sigma},\bm{k}}c_{\widetilde{\beta^{\prime}\sigma^{\prime}},\bm{k}^{\prime}}\\\
&=\int
d\bm{k}d\bm{k}^{\prime}\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left[\sum_{\bm{G}_{l}}e^{i\bm{G}_{l}\cdot\bm{r}_{j}}t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{C_{3}^{-1}\bm{k},C_{3}^{-1}\bm{k}^{\prime},\bm{G}_{l}}\right]e^{iC_{3}^{-1}(-\bm{k}+\bm{k}^{\prime})\cdot\bm{r}_{j}}c^{\dagger}_{\widetilde{\alpha\sigma},\bm{k}}c_{\widetilde{\beta^{\prime}\sigma^{\prime}},\bm{k}^{\prime}}\\\
&=\int
d\bm{k}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{\bm{G}_{l}}t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{C_{3}^{-1}\bm{k},C_{3}^{-1}\bm{k}-\bm{G}_{l},\bm{G}_{l}}~{}c^{\dagger}_{\widetilde{\alpha\sigma},\bm{k}}c_{\widetilde{\beta^{\prime}\sigma^{\prime}},\bm{k}-C_{3}\bm{G}_{l}}\\\
&=\int
d\bm{k}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{\bm{G}_{l}}t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}~{}c^{\dagger}_{\widetilde{\alpha\sigma},C_{3}\bm{k}}c_{\widetilde{\beta^{\prime}\sigma^{\prime}},C_{3}\bm{k}-C_{3}\bm{G}_{l}},\end{split}$
(34)

where the orbital-spin index $\widetilde{\alpha\sigma}$ means the $C_{3z}$
rotated orbital-spin.

The conditions for this to be $C_{3}$ symmetric (Eq. (23)) is

$t^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}=t^{\widetilde{\alpha\sigma},\widetilde{\beta^{\prime}\sigma^{\prime}}}_{C_{3}\bm{k},C_{3}\bm{k}-C_{3}\bm{G}_{l},C_{3}\bm{G}_{l}}.$
(35)

When the Hamiltonian of the untwisted structures are $C_{3}$ symmetric,

$h_{\bm{r}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})=h_{\bm{r}}^{\widetilde{\alpha\sigma},\widetilde{\beta^{\prime}\sigma^{\prime}}}(C_{3}\bm{k})$
(36)

is satisfied and the Eq. (23) is satisfied by using Eq. (29).

  * •

Time-reversal symmetry

The time-reversed Hamiltonian is

$\begin{split}{\cal T}H{\cal
T}^{\dagger}=&\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{pq}\left.t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\right.^{*}\left(\bm{r}_{i}\right)~{}c^{\dagger}_{\alpha\bar{\sigma}}\left(\bm{r}_{j}+\bm{d}_{\alpha}\left(\bm{r}_{j}\right)\right)c^{\phantom{\dagger}}_{\beta^{\prime}\bar{\sigma}^{\prime}}\left(\bm{r}_{j}+\bm{d}_{\beta^{\prime}}\left(\bm{r}_{j}\right)+\bm{a}_{pq}\right)(-1)^{1-\delta_{\sigma\sigma^{\prime}}}\\\
=&\int
d\bm{k}d\bm{k}^{\prime}\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\left[\sum_{pq}t_{pq}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{r}_{j})e^{i\bm{k}\cdot\bm{d}_{\alpha}(\bm{r}_{j})}e^{i\bm{k}^{\prime}\cdot\left(\bm{d}_{\beta^{\prime}}(\bm{r}_{j})+\bm{a}_{pq}\right)}\right]^{*}e^{i(-\bm{k}+\bm{k}^{\prime})\cdot\bm{r}_{j}}~{}c^{\dagger}_{\alpha\bar{\sigma},\bm{k}}c^{\phantom{\dagger}}_{\beta^{\prime}\bar{\sigma}^{\prime},\bm{k}^{\prime}}(-1)^{1-\delta_{\sigma\sigma^{\prime}}}\\\
=&\int
d\bm{k}d\bm{k}^{\prime}\sum_{\bm{r}_{j}}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{\bm{G}_{l}}e^{-i\bm{G}_{l}\cdot\bm{r}_{j}}\left.t_{-\bm{k},-\bm{k}^{\prime},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\right.^{*}e^{i(-\bm{k}+\bm{k}^{\prime})\cdot\bm{r}_{j}}~{}c^{\dagger}_{\alpha\bar{\sigma},\bm{k}}c^{\phantom{\dagger}}_{\beta^{\prime}\bar{\sigma}^{\prime},\bm{k}^{\prime}}(-1)^{1-\delta_{\sigma\sigma^{\prime}}}\\\
=&\int
d\bm{k}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{\bm{G}_{l}}\left.t_{-\bm{k},-\bm{k}-\bm{G}_{l},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\right.^{*}~{}c^{\dagger}_{\alpha\bar{\sigma},\bm{k}}c^{\phantom{\dagger}}_{\beta^{\prime}\bar{\sigma}^{\prime},\bm{k}+\bm{G}_{l}}(-1)^{1-\delta_{\sigma\sigma^{\prime}}},\end{split}$
(37)

where $\bar{\sigma}$ means the opposite spin.

The conditions for this to be time-reversal symmetric (Eq. (24)) is

$(-1)^{1-\delta_{\sigma\sigma^{\prime}}}\left.t_{-\bm{k},-\bm{k}-\bm{G}_{l},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\right.^{*}=t_{\bm{k},\bm{k}+\bm{G}_{l},-\bm{G}_{l}}^{\alpha\bar{\sigma},\beta^{\prime}\bar{\sigma}^{\prime}}.$
(38)

When the Hamiltonian of the untwisted structures are time-reversal symmetric,

$\left.h_{\bm{r}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})\right.^{*}=(-1)^{1-\delta_{\sigma\sigma^{\prime}}}~{}h_{\bm{r}}^{\alpha\bar{\sigma},\beta^{\prime}\bar{\sigma}^{\prime}}(-\bm{k})$
(39)

is satisfied and the Eq. (24) is satisfied by using Eq. (29).

##  Appendix C Numerical calculation methods

In this section, we describe methods used in our numerical calculations.

###  C.1 Matrix representation

We have obtained the Hamiltonian of the Moiré system,

$\begin{split}H&=\int
d\bm{k}\sum_{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\sum_{\bm{G}_{l}}t_{\bm{k},\bm{k}-\bm{G}_{l},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}~{}c_{\alpha\sigma,\bm{k}}^{\dagger}c_{\beta^{\prime}\sigma^{\prime},\bm{k}-\bm{G}_{l}}^{\phantom{\dagger}},\end{split}$
(40)

and here we give the matrix representation of this Hamiltonian so that it can
be implemented in numerical calculations.

Because all hopping processes have momentum difference $\bm{G}_{l}$, the basis
is given for $\bm{k}\in$ Moiré BZ as

$\left\\{c_{\alpha\sigma,\bm{k}+\bm{v}_{mn}}\left|\bm{v}_{mn}=m\bm{G}_{1}+n\bm{G}_{2},~{}m,n\in\mathbb{Z}\right.\right\\}.$
(41)

To describe low-energy physics with a finite dimension matrix, we take a
momentum cutoff $k_{c}$ as an approximation. In the Bi2Te3 system, the low-
energy physics are described around the $\Gamma$ point. Therefore, we take a
finite base set as

$\left\\{c_{\alpha\sigma,\bm{k}+\bm{v}_{mn}}\left|\bm{v}_{mn}=m\bm{G}_{1}+n\bm{G}_{2},~{}m,n\in\mathbb{Z},~{}|\bm{v}_{mn}|<k_{c}\right.\right\\}.$
(42)

With this basis, the Matrix representation of the Hamiltonian is

$\begin{split}H(\bm{k})=\left(\begin{array}[]{cccccc}t_{\bm{k},\bm{k},\bm{0}}&&&&\mathrm{if}~{}\bm{v}_{mn}-\bm{v}_{m^{\prime}n^{\prime}}=\bm{G}_{l}&\\\
&\ddots&&&\downarrow&\\\
&&t_{\bm{k}+\bm{v}_{mn},\bm{k}+\bm{v}_{mn},\bm{0}}&\cdots&t_{\bm{k}+\bm{v}_{mn},\bm{k}+\bm{v}_{m^{\prime}n^{\prime}},\bm{G}_{l}}&\\\
&&&\ddots&\vdots&\\\
&&&&t_{\bm{k}+\bm{v}_{m^{\prime}n^{\prime}},\bm{k}+\bm{v}_{m^{\prime}n^{\prime}},\bm{0}}&\\\
&&&&&\ddots\\\ \end{array}\right),\end{split}$ (43)

where $\bm{k}\in$ Moiré BZ.

Because we consider only seven $\bm{G}_{l}$, an off-diagonal block is non-zero
only when $\bm{v}_{mn}-\bm{v}_{m^{\prime}n^{\prime}}$ matches one of the seven
$\bm{G}_{l}$.

###  C.2 Validity of the cutoff $k_{c}$

The validity of the cutoff $k_{c}$ is evaluated as follows. Once we calculate
a wave function $\psi_{\bm{k}}$, and make a weight mapping onto Moiré BZs
within $k_{c}$ as

$\begin{split}\psi_{\bm{k}}=&\left(\begin{array}[]{c}\phi_{1}\\\ \phi_{2}\\\
\vdots\end{array}\right)\\\
&\mapsto\left(\begin{array}[]{c}\sum_{l=1}^{60}|\phi_{l}|^{2}\\\
\sum_{l=61}^{120}|\phi_{l}|^{2}\\\
\vdots\end{array}\right)~{}\begin{array}[]{l}\leftarrow\mathrm{weight~{}on~{}Moir\acute{e}~{}BZ\\#1}\\\
\leftarrow\mathrm{weight~{}on~{}Moir\acute{e}~{}BZ\\#2}\\\
\vdots\end{array}.\end{split}$ (44)

By plotting the weight mapping, we check the weight on the BZs around $k_{c}$.
If the BZs around $k_{c}$ have small weights, the cutoff $k_{c}$ is large
enough and the effect of $k_{c}$ is negligible. If the BZs around $k_{c}$ have
significant weight compared to other BZs, larger $k_{c}$ should be taken.
Generally, wave functions whose energies are far from Fermi level requires
larger $k_{c}$. Therefore, evaluating with wave functions that are distant
from the Fermi level ensures the validity of $k_{c}$ for wave functions that
are closer to the Fermi level. Figure 12 is an example of the weight mapping
for the 17th valence band of Bi2(Te0.92Se0.08)3 with a twist angle
$\theta=1.0^{\circ}$, cutoff $k_{c}=\pi/4$. Each hexagon represents a Moiré
BZ, and its color represents the weight (normalized with the maximum value).
It can be seen that the Moiré BZs around $k_{c}$ have almost no weights and
thus the cutoff is large enough at least for the 1st $\sim$ 17th valence bands
in this case. The validity of $k_{c}$ has been confirmed in all calculations
in this paper.

Figure 12: Weight mapping of the 17th valence band at the $\Gamma$ point of
twisted bilayer Bi2(Te0.92Se0.08)3 with a twist angle $\theta=1.0^{\circ}$ and
cutoff $k_{c}=\pi/4$. The hexagons are Moiré BZs and their colors represent
the weight on them (normalized with the maximum value). In this case, it can
be seen that the cutoff is large enough.

###  C.3 Periodicity in the momentum space

From the Bloch theorem, a wave function $\psi_{\bm{k}+\bm{G}_{l}}$ is
physically identical with $\psi_{\bm{k}}$, and thus
$|\braket{\psi_{\bm{k}+\bm{G}_{l}}}{\psi_{\bm{k}}}|=1$. This fact is used in
the numerical evaluation of the Wilson loop spectra.

In the effective model of a Moiré superlattice system with basis (42), the
momentum shift $\bm{k}\to\bm{k}+\bm{G}_{l}$ effectively result in a shift in
$\bm{v}_{mn}$ indices. Therefore, components of a numerically obtained
eigenvector $\psi_{\bm{k}+\bm{G}_{l}}$ are shifted from those of
$\psi_{\bm{k}}$, and $|\braket{\psi_{\bm{k}+\bm{G}_{l}}}{\psi_{\bm{k}}}|\neq
1$ due to the $\bm{k}$ dependence of the basis. Note that, as long as the
cutoff $k_{c}$ is large enough, i.e. no significant component is put outside
of the cutoff by the shift, $\psi_{\bm{k}+\bm{G}_{l}}$ and $\psi_{\bm{k}}$ are
physically identical states. To evaluate the inner product of wave functions
$\psi_{\bm{k}+\bm{G}_{l}}$ and $\psi_{\bm{k}}$ correctly with the numerically
obtained eigenvectors, we introduce a base-shift matrix $T$, which is defined
as

$T=\\{T_{(mn)(m^{\prime}n^{\prime})}\\}=\left\\{\begin{array}[]{ll}I_{\eta}&\bm{v}_{mn}+\bm{G}_{l}=\bm{v}_{m^{\prime}n^{\prime}}\\\
O_{\eta}&\mathrm{others}\\\ \end{array}\right.,$ (45)

where $\eta$ is the orbital-spin degree of freedom ($\eta=60$ in twisted
bilayer Bi2(Te1-xSex)3), and $I_{\eta}$ and $O_{\eta}$ are the $\eta$
dimension identity matrix and zero matrix, respectively. By using this $T$,
the inner product of wave functions $\psi_{\bm{k}+\bm{G}_{l}}$ and
$\psi_{\bm{k}}$ is evaluated as
$\braket{\psi_{\bm{k}+\bm{G}_{l}}}{T}{\psi_{\bm{k}}}$, which approximately
satisfies $|\braket{\psi_{\bm{k}+\bm{G}_{l}}}{T}{\psi_{\bm{k}}}|=1$ when the
cutoff is large enough.

###  C.4 Wannierization and surface state calculation

To calculate surface states, we construct Wannier functions from wave
functions and a band dispersion obtained from the effective model. We use the
option use_bloch_phases=T of the WANNIER90 [46] to set the Bloch functions as
the initial guess for the projections. In calculations with this option, we
need a list of eigenenergies $E_{n\bm{k}}$ in file prefix.eig and a list of
wave function overlaps $\braket{\psi_{\bm{k}+\delta\bm{k}}}{\psi_{\bm{k}}}$ in
file prefix.mmn for the target bands. These lists are numerically obtained
from the effective model calculation. Both of them should be calculated on a
finite $\bm{k}$ mesh. Note that in the calculation of overlaps, the base-shift
matrix given in the previous subsection should be used for overlaps that cross
the end of the $\bm{k}$ mesh. With this method, Wannier functions and hoppings
between them are obtained for Moiré bands. With the obatined parameters, edge
state spectra are calculated by the Green’s function method implemented in the
WannierTools package [50].

For the parameters obtained by the wannierization, we make a correction to
recover the time-reversal symmetry. It is because the energy scale of a Moiré
superlattice system is so small that numerical errors in the first-principles
calculations in untwisted systems get more noticeable. The correction is done
by replacing the hopping parameters as

$t_{\bm{a}_{nm}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}\to\frac{1}{2}\left(t_{\bm{a}_{nm}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}+(-1)^{1-\delta_{\sigma\sigma^{\prime}}}\left.t_{\bm{a}_{nm}}^{\alpha\bar{\sigma},\beta^{\prime}\bar{\sigma}^{\prime}}\right.^{*}\right),$
(46)

where $\bar{\sigma}$ is the opposite spin of $\sigma$. The
$\delta_{\sigma\sigma^{\prime}}$ is the Kronecker delta for spin $\sigma$ and
$\sigma^{\prime}$, which is essentially comes from ${\cal T}^{2}=-1$ in a
spin-$1/2$ system. We have confirmed that this correction works only as a
minor correction in the edge state spectra.

###  C.5 Symmetry in wave function plot

Twisted bilayer Bi2(Te1-xSex)3 has the time-reversal symmetry and the layer
group symmetry #67, which include $C_{3z}$ rotation symmetry. Due to the time-
reversal symmetry, all bands at the $\Gamma$ point appear as doubly degenerate
Kramers’ pairs. Because the numerical diagonalization does not necessarily
give simultaneous eigenvectors of the Hamiltonian and the $C_{3z}$ rotation,
we need to symmetrize the Kramers’ pair to obtain a $C_{3z}$ symmetric wave
function plot at the $\Gamma$ point.

For a given Kramers’ pair $\ket{\psi_{1}},\ket{\psi_{2}}$ and $C_{3z}$
operator, a matrix representation of $\Gamma_{C_{3z}}$ is given as

$\Gamma_{C_{3z}}=\left(\begin{array}[]{cc}\braket{\psi_{1}}{C_{3z}}{\psi_{1}}&\braket{\psi_{1}}{C_{3z}}{\psi_{2}}\\\
\braket{\psi_{2}}{C_{3z}}{\psi_{1}}&\braket{\psi_{2}}{C_{3z}}{\psi_{2}}\end{array}\right).$
(47)

With a Unitary transform $U$ that diagonalize $\Gamma_{C_{3z}}$ as
$U^{\dagger}\Gamma_{C_{3z}}U=\mathrm{diag}(\lambda_{1},\lambda_{2})$,
symmetrized wave functions $\ket{\tilde{\psi}_{1}},\ket{\tilde{\psi}_{2}}$ are
obtained as

$\left(\begin{array}[]{cc}\ket{\tilde{\psi}_{1}}&\ket{\tilde{\psi}_{2}}\end{array}\right)=\left(\begin{array}[]{cc}\ket{\psi_{1}}&\ket{\psi_{2}}\end{array}\right)U.$
(48)

##  Appendix D Effect of Se-doping on twisted bilayer Bi2(Te1-xSex)3

In this section, we discuss the effect of Se-doping in in twisted bilayer
Bi2(Te1-xSex)3. First, we see the Se-doping effect in parameters of untwisted
Bi2(Te1-xSex)3 systems. Next, we see the Moiré band dispersions for various
Se-doping values.

###  D.1 Effect of Se-doping on parameters in untwisted bilayer Bi2(Te1-xSex)3

The discrete Fourier transform Eq. (25) interpolates the untwisted
Hamiltonians to the intermediate area. Therefore, we can define “bandgap in
each real space position $\bm{r}$” in the Moiré unit cell, and it can be used
to estimate where the topological domain boundary is. To evaluate the bandgap
in the $\Gamma$ point in the small twist angle limit, we consider
$\bm{k}=\bm{0}$ in Eq. (25) as

$\sum_{\bm{G}_{l}}e^{i\bm{G}_{l}\cdot\bm{r}}~{}t_{\bm{0},-\bm{G}_{l},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}=h_{\bm{r}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{0}),$
(49)

where
$t_{\bm{0},-\bm{G}_{l},\bm{G}_{l}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}$
are given in Eq. (29). Note that $\bm{G}_{l}$ are negligible compared to the
atomic BZ in the small twist angle limit. By diagonalizing the right-hand-
side, we obtain a bandgap $\Delta E(\bm{r})$. The position where $\Delta
E(\bm{r})=0$ is satisfied is the topological domain boundary. In Fig.13, the
energies of the valence top band and the conduction bottom band at the
$\Gamma$ point calculated with Eq.(49) are plotted along the AA-AB-AC-AA
stacking line for Bi2Te3, Bi2(Te1-xSex)3 with various $x$, Bi2Se3, and Sb2Te3.
Here, we calculate the electronic band structure of Sb2Te3 under the same
condition as Bi2Te3 and Bi2Se3. Note that Bi2Se3 and Sb2Te3 cases are plotted
in a different energy range. $\Delta E(\bm{r})=0$ is the gap between the two
band in $\bm{r}$ (horizontal axis).

Figure 13: Energies of the valence top band and the conduction bottom band at
the $\Gamma$ point plotted along the AA-AB-AC-AA stacking line for Bi2Te3,
Bi2(Te1-xSex)3 with various $x$, Bi2Se3, and Sb2Te3. Note that the Bi2Se3 and
Sb2Te3 cases are plotted in a different energy range.

As shown in the main article, the valence and conduction bands in Bi2Te3 are
inverted only around the AB-stacking case. As Se is doped ($x$ gets larger),
the (inverted) bandgap in the AB-stacking area gets smaller, and the bandgap
in the AA- and AC-stacking are get larger. Around Bi2(Te0.60Se0.40)3
($x=0.40$), a band-touching occurs at the AB-stacking point, and the band
structure in the AB-stacking area become trivial in the $x>0.40$ region.
Because the plots of Bi2Se3 and Sb2Te3 are almost identical, it can be seen
that replacing Bi with Sb has a similar effect as replacing Te with Se.

###  D.2 Effect of Se-doping on Moiré band dispersion

Figure 14 shows Moiré band dispersions of twisted bilayer Bi2(Te1-xSex)3 for
various $x$. The twist angle is fixed at $\theta=1.00^{\circ}$. The values of
$x$ are the same as those in Fig.13. After the gap-closing at the AB-stacking
occurs ($x>0.40$), the edge-state-originated bands disappear and the band gap
around the Fermi level get larger. For twisted bilayer Bi2Se3, a wave function
of the valence top band at the $\Gamma_{M}$ point is plotted. The wave
function is found to be a bulk-originated state around the AB-stacking area.

Figure 14: Moiré band dispersions of twisted bilayer Bi2(Te1-xSex)3 for
various $x$. The twist angle is fixed at $\theta=1.00^{\circ}$. Note that the
Bi2Se3 case is plotted in a different energy range. A wave function of the
valence top band at the $\Gamma_{M}$ point (upper layer, lower Bi atom,
$p_{z}$ orbital, spin-up component) is also shown.

##  Appendix E Supplemental Information of twisted bilayer Bi2(Te1-xSex)3

In this section, we give supplemental information of twisted bilayer
Bi2(Te1-xSex)3.

###  E.1 Angular momentum of edge-state-originated bands

Figure 15: Example of the wave function and angular momentum in the VB2 of
twisted bilayer Bi2(Te0.85Se0.15)3 with $\theta=0.50^{\circ}$. The wave
function plots shows the components of spin-up and spin-down on (upper layer,
lower Bi atom, $p_{z}$ orbital). The total angular momentum of this wave
function is calculated as 3/2.

When a wave function is plotted for an edge-state-originated band, it appears
to have an integer angular momentum. However, the electronic state is coupled
with spin degree-of-freedom and thus the total angular momentum is a half-
integer. Actually, when both of spin-up and spin-down components are plotted
for a wave function, the phase winding number ($\approx$ orbital angular
momentum) is different by 1 (Fig.15). By coupling the orbital and spin angular
momentum, a total angular momentum (or a rotation eigenvalue of the wave
function) is uniquely defined.

###  E.2 Twist angle dependence of the bandwidth

Figure 16: Twist angle dependence of the bandwidth of VB1 (orange), VB2 (red),
VB3 (green), and VB4 (purple) in twisted bilayer Bi2(Te0.85Se0.15)3. The band
width is estimated with the band dispersion along the MM-$\Gamma_{M}$-KM-MM
line.

Figure 16 shows the twist angle dependence of the bandwidth of VB1-4 in
twisted bilayer Bi2(Te0.85Se0.15)3. The bandwidth is estimated as
$\mathrm{Max}(E(\bm{k}))-\mathrm{Min}(E(\bm{k})$)
($\bm{k}\in\\{$MM-$\Gamma_{M}$-KM-MM line$\\}$). As a general trend, the
bandwidth tends to increase as the twist angle increases. This tendency is
explained by the band folding in the Moiré BZ, as explained in the main
article. More detailed behavior is determined by a combination of multiple
factors such as Rashba splitting and hybridization with other bands.

###  E.3 Edge dependence of Moiré edge state spectra

Figure 17: Edge dependence of edge state spectra. (a) The edge BZ used in the
edge state calculations. (b),(c), Edge state spectra of the right edge and
left edge in twisted bilayer Bi2(Te0.92Se0.08)3.

Here, we compare the edge state spectra for the left and right edge. Note that
the left and right edges are defined by using the Moiré unit cell as a unit to
build the half-infinite plane, i.e. the the left edge (right edge) of the
Moiré unit cell appears on the left edge (right edge) of the half-infinite
plane. Figure 17 shows the edge BZ (a) and edge state spectra of the right
edge (b) and left edge (c) in twisted bilayer Bi2(Te0.92Se0.08)3. The edge
state spectra are topologically identical for the both edge, but the energies
of the Dirac cones of the edge states are different. It is because there is no
symmetry to guarantee the equivalence of their energies in twisted bilayer
Bi2(Te0.92Se0.08)3. Although the Moiré superlattice of twisted bilayer
Bi2(Te1-xSex)3 has in-plane $C_{2}$ rotation symmetries along the
$\Gamma_{M}$-MM lines, they do not exchange the two edges. Therefore, each
edge state does not necessarily have the same spectrum. This edge dependence
is not Moiré-specific but generally found in a noncentrosymmetric 2D
topological insulator. See Sec.F for general and more detailed discussion.
Note that in more realistic conditions, the energy difference between the two
edge states may be suppressed by the charge neutral condition.

###  E.4 Topological phase transition in Moiré bands

Figure 18: Bulk and edge band spectra of twisted bilayer Bi2(Te0.92Se0.08)3
with $\theta=1.35^{\circ}$. The gap-closing of the bulk band spectra occurs at
general momenta (white dashed circles).

Here, we show how a gap-closing and a topological phase transition occur in
Moiré bands. As the twist angle get smaller, more edge-state-originated nearly
flat bands appear around the Fermi level due to the band folding with the
small Moiré BZ. For example in the case of twisted bilayer Bi2(Te0.92Se0.08)3
with $\theta\sim 1.55^{\circ}$, topological phase transitions occur on VB8 and
VB9 and eventually they become topologically trivial bands as the twist angle
is decreased. Generally for the topological phase transition in twisted
bilayer Bi2(Te1-xSex)3, the gap-closing occurs on the $\Gamma_{M}$-KM lines
due to the in-plane $C_{2}$ rotation symmetry [66]. Note that the
$\Gamma_{M}$-KM lines are not high-symmetry lines, although the position of
the gap-closing is determined by the symmetry restriction. When edges are made
so that there is no symmetry between the two edges, the gap-closing is also
seen in the 1D edge BZ as a gap-closing of bulk-band spectra at general
momenta. Figure 18 shows a bulk band spectra and an edge band spectra of
twisted bilayer Bi2(Te0.92Se0.08)3 with the twist angle $\theta=1.35^{\circ}$,
which is close to a phase transition point. Gap-closing points are indicated
with white dashed circles. Note that the edge states appear from the points
where the gap-closing will occur and make a Dirac-cone that is close to the
bulk-band crossing of a Kramers pair at the $\bar{\Gamma}_{M}$ point. The
topological phase transition that occurs as the twist angle is decreased is
not simple as a pair annihilation between VB8 and VB9. As the twist angle is
decreased, the gap between edge-state-originated bands gets smaller and thus
they shift relatively upward, while the bulk-state-originated band stays.
Therefore, the nontrivial structure moves to the next and next band as the
twist angle is decreased. Actually, another gap-closing occurs between VB9 and
VB10 and there is midgap edge states in the lower energy region VB9 in Fig.18.

##  Appendix F Noncentrosymmetric 2D topological phase transition

A Moiré system does not have the inversion symmetry and we have found that a
topological phase transition occurs in twisted Bi2(Te1-xSex)3. In this
section, we explain how a topological phase transition occurs in a general
noncentrosymmetric 2D system. We consider a topological phase transition
between a topological insulator (TI) phase and a normal insulator (NI) phase
in time-reversal ${\cal T}$ symmetric systems. Generally, when a topological
phase transition occurs, a gap-closing must occurs in the bulk band dispersion
[51, 52]. First, we briefly review the case of a centrosymmetric TI. In a
phase transition occurs from a NI to a TI in a centrosymmetric system, the
bulk gap-closing occurs at one of (or an odd number of) the TRIM. When a
surface band spectrum is drawn, the gap-closing is seen as a gap-closing of
the bulk continuous spectra at one of the TRIM of the surface BZ, because TRIM
of 3D BZ are projected on TRIM of the surface BZ. In the TI phase, a
corresponding surface state emerges around the TRIM in which the gap-closing
occurred. In noncentrosymmetric systems, the bulk gap-closing occurs at a
generic $\bm{k}$ point. In 3D case, the gapless phase has finite width in the
parameter space, which is Weyl semimetal (WSM) phase. In 2D case, the gap-
closing occurs only on the critical point in the parameter space. Because the
bulk gap-closing occurs at a generic $\bm{k}$ point, a gap-closing in the
surface band spectrum should also be seen at a generic $\bm{k}$ in the surface
BZ. Ref.[66] clarified that, in noncentrosymmetric 2D systems, $\bm{k}$ points
where the gap-closing can occur are restricted by crystalline symmetries.
However, it has not been clearly described how a corresponding surface state
emerges from the gap-closing point in the surface band spectrum. Therefore,
here we demonstrate the emergence of the surface state in topological phase
transition in noncentrosymmetric 2D systems, and explain the relation between
the surface state and the crystalline symmetries.

As explained above, in 3D system, a topological phase transition in a
noncentrosymmetric system is described as a NI-WSM-TI transition. Therefore,
we start from a 3D model that describe a NI-WSM-TI transition, and by focusing
a particular 2D momentum plane in the model we discuss the bulk and surface
band spectra in a noncentrosymmetric 2D system. We use the inversion broken
TI/NI multilayer model given in Refs.[67] and [68], and it is a four by four
model written as

$\begin{split}H=&v\tau_{z}\left(k_{y}\sigma_{x}-k_{x}\sigma_{y}\right)+V\tau_{z}\\\
&+\Delta_{N}(\bm{k})\tau_{x}+\Delta_{T}(\bm{k})\left(e^{ik_{z}}\tau_{+}+e^{-ik_{z}}\tau_{-}\right).\end{split}$
(50)

The four components consist of surface states on the top surface of the TI
layer with spin up and down ($t_{\uparrow\bm{k}},t_{\downarrow\bm{k}}$), and
those on the bottom surface ($b_{\uparrow\bm{k}},b_{\downarrow\bm{k}}$). The
$\sigma_{x,y,z}$ and $\tau_{x,y,z}$ are Pauli matrices for spin and top/bottom
surfaces, respectively, where $\tau_{\pm}=\left(\tau_{x}\pm
i\tau_{y}\right)/\sqrt{2}$. We assume a two-fold rotation symmetry about the
$z$-axis ($C_{2z}$), and thus the $\Delta_{N}(\bm{k})$ and
$\Delta_{T}(\bm{k})$ are written as

$\begin{split}\Delta_{T}\left(\bm{k}\right)&=\Delta_{T}^{0}+\delta_{T}^{x}k_{x}^{2}+\delta_{T}^{y}k_{y}^{2},\\\
\Delta_{N}\left(\bm{k}\right)&=\Delta_{N}^{0}+\delta_{N}^{x}k_{x}^{2}+\delta_{N}^{y}k_{y}^{2}.\\\
\end{split}$ (51)

To draw the surface band spectrum, we transform the model into an equivalent
lattice periodic model. The lattice periodic model is

$\begin{split}H=&v\tau_{z}\left(\sin k_{y}\sigma_{x}-\sin
k_{x}\sigma_{y}\right)+V\tau_{z}\\\
&+\Delta_{T}(\bm{k})\tau_{x}+\Delta_{N}(\bm{k})\left(e^{ik_{z}}\tau_{+}+e^{-ik_{z}}\tau_{-}\right),\end{split}$
(52)

where $\Delta_{N}(\bm{k})$ and $\Delta_{T}(\bm{k})$ are written as

$\begin{split}\Delta_{T}\left(\bm{k}\right)&=\Delta_{T}^{0}+2\delta_{T}^{x}(1-\cos
k_{x})+2\delta_{T}^{y}(1-\cos k_{y}),\\\
\Delta_{N}\left(\bm{k}\right)&=\Delta_{N}^{0}+2\delta_{N}^{x}(1-\cos
k_{x})+2\delta_{N}^{y}(1-\cos k_{y}).\\\ \end{split}$ (53)

By taking $\Delta_{T}^{0}$ as a tuning parameter, this model shows the NI-WSM-
TI transition. We fix the other parameters as

$\begin{split}&V=0.2~{},~{}v=0.4~{},~{}\delta_{N}^{x}=0.60~{},~{}\delta_{T}^{x}=0.20,\\\
&\delta_{N}^{y}=\delta_{N}^{x}\left(1+\frac{0.1}{\delta_{T}^{x}-\delta_{N}^{x}}\frac{1/2}{1-\sqrt{1-V^{2}/v^{2}}}\right),\\\
&\delta_{T}^{y}=\delta_{T}^{x}\left(1+\frac{0.1}{\delta_{T}^{x}-\delta_{N}^{x}}\frac{1/2}{1-\sqrt{1-V^{2}/v^{2}}}\right),\\\
&\Delta_{N}^{0}=1.0+2\left(\delta_{T}^{y}-\delta_{N}^{y}\right)\left(1-\sqrt{1-V^{2}/v^{2}}\right).\end{split}$
(54)

Here, $\delta^{y}_{N}$, $\delta^{y}_{T}$, and $\Delta^{0}_{N}$ are determined
to fix the transition points at $\Delta^{0}_{T}=1.0$ and $1.1$. When
$\Delta_{T}^{0}<1.0$, the model is a normal insulator (Fig.19). When
$\Delta_{T}^{0}=1.0$, Weyl points are created as pairs on the
$(k_{z}=\pi,k_{x}=0)$ line (Fig.19). Due to the time-reversal symmetry, a pair
of Weyl points are created each of positive and negative $k_{y}$ region and
thus four Weyl points are created in total. When $1.0<\Delta_{T}^{0}<1.1$, the
Weyl points move on the $k_{z}=\pi$ plane, and when $\Delta_{T}^{0}=1.1$, they
annihilate on the $(k_{z}=\pi,k_{y}=0)$ line. When $\Delta_{T}^{0}>1.1$, the
model is a strong topological insulator.

Figure 19: Phase diagram of 3D WSM model (52). A WSM phase appears between
$\Delta_{T}^{0}=1.0$ and $\Delta_{T}^{0}=1.1$ as a gapless phase between NI
and TI. At $\Delta_{T}^{0}=1.0$, Weyl points are created on the
$(k_{x}=0,k_{z}=\pi)$ line. As $\Delta_{T}^{0}$ increases, the Weyl points
move in the $k_{z}=\pi$ plane and annihilate on the $(k_{y}=0,k_{z}=\pi)$ line
at $\Delta_{T}^{0}=1.1$.

By focusing on a particular 2D momentum space in the 3D BZ of this model, we
can discuss a topological phase transition and surface states in a
noncentrosymmetric 2D system. To consider a time-reversal symmetric system,
the focused 2D plane must be a time-reversal invariant plane
($\forall\bm{k}\in\left\\{\mathrm{2D~{}plane}\right\\},-\bm{k}\in\left\\{\mathrm{2D~{}plane}\right\\}$).
Here we focus on two planes, $k_{x}=0$ and $k_{y}=0$. To draw the surface band
spectrum with Green’s function method, we transform the lattice periodic model
into a real space representation as

$\begin{split}H=\sum_{\bm{r}}&\left[i\frac{v}{2}\left(-t_{\uparrow\bm{r}}^{\dagger}t_{\downarrow\bm{r}+\bm{y}}^{\phantom{\dagger}}+t_{\uparrow\bm{r}}^{\dagger}t_{\downarrow\bm{r}-\bm{y}}^{\phantom{\dagger}}+b_{\uparrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}+\bm{y}}^{\phantom{\dagger}}-b_{\uparrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}-\bm{y}}^{\phantom{\dagger}}\right)\right.\\\
&~{}+\frac{v}{2}\left(t_{\uparrow\bm{r}}^{\dagger}t_{\downarrow\bm{r}+\bm{x}}^{\phantom{\dagger}}-t_{\uparrow\bm{r}}^{\dagger}t_{\downarrow\bm{r}-\bm{x}}^{\phantom{\dagger}}-b_{\uparrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}+\bm{x}}^{\phantom{\dagger}}+b_{\uparrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}-\bm{x}}^{\phantom{\dagger}}\right)\\\
&~{}+V\left(t_{\uparrow\bm{r}}^{\dagger}t_{\uparrow\bm{r}}^{\phantom{\dagger}}+t_{\downarrow\bm{r}}^{\dagger}t_{\downarrow\bm{r}}^{\phantom{\dagger}}-b_{\uparrow\bm{r}}^{\dagger}b_{\uparrow\bm{r}}^{\phantom{\dagger}}-b_{\downarrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}}^{\phantom{\dagger}}\right)\\\
&~{}+\left(\Delta_{N}^{0}+2\delta_{N}^{x}+2\delta_{N}^{y}\right)\left(t_{\uparrow\bm{r}}^{\dagger}b_{\uparrow\bm{r}}^{\phantom{\dagger}}+t_{\downarrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}}^{\phantom{\dagger}}\right)\\\
&~{}-\delta_{N}^{x}\left(t_{\uparrow\bm{r}}^{\dagger}b_{\uparrow\bm{r}+\bm{x}}^{\phantom{\dagger}}+t_{\downarrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}+\bm{x}}^{\phantom{\dagger}}+t_{\uparrow\bm{r}}^{\dagger}b_{\uparrow\bm{r}-\bm{x}}^{\phantom{\dagger}}+t_{\downarrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}-\bm{x}}^{\phantom{\dagger}}\right)\\\
&~{}-\delta_{N}^{y}\left(t_{\uparrow\bm{r}}^{\dagger}b_{\uparrow\bm{r}+\bm{y}}^{\phantom{\dagger}}+t_{\downarrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}+\bm{y}}^{\phantom{\dagger}}+t_{\uparrow\bm{r}}^{\dagger}b_{\uparrow\bm{r}-\bm{y}}^{\phantom{\dagger}}+t_{\downarrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}-\bm{y}}^{\phantom{\dagger}}\right)\\\
&~{}+\left(\Delta_{T}^{0}+2\delta_{T}^{x}+2\delta_{T}^{y}\right)\left(t_{\uparrow\bm{r}}^{\dagger}b_{\uparrow\bm{r}+\bm{z}}^{\phantom{\dagger}}+t_{\downarrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}+\bm{z}}^{\phantom{\dagger}}\right)\\\
&~{}-\delta_{T}^{x}\left(t_{\uparrow\bm{r}}^{\dagger}b_{\uparrow\bm{r}+\bm{x}+\bm{z}}^{\phantom{\dagger}}+t_{\downarrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}+\bm{x}+\bm{z}}^{\phantom{\dagger}}\right.\\\
&~{}~{}~{}~{}~{}~{}~{}\left.+t_{\uparrow\bm{r}}^{\dagger}b_{\uparrow\bm{r}-\bm{x}+\bm{z}}^{\phantom{\dagger}}+t_{\downarrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}-\bm{x}+\bm{z}}^{\phantom{\dagger}}\right)\\\
&~{}-\delta_{T}^{y}\left(t_{\uparrow\bm{r}}^{\dagger}b_{\uparrow\bm{r}+\bm{y}+\bm{z}}^{\phantom{\dagger}}+t_{\downarrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}+\bm{y}+\bm{z}}^{\phantom{\dagger}}\right.\\\
&~{}~{}~{}~{}~{}~{}~{}\left.\left.+t_{\uparrow\bm{r}}^{\dagger}b_{\uparrow\bm{r}-\bm{y}+\bm{z}}^{\phantom{\dagger}}+t_{\downarrow\bm{r}}^{\dagger}b_{\downarrow\bm{r}-\bm{y}+\bm{z}}^{\phantom{\dagger}}\right)\right]\\\
&~{}+\mathrm{h.c.},\end{split}$ (55)

where $\bm{r}$ is the unit cell index and a cubic cell with a lattice constant
1 is assumed.

Figure 20: Phase diagrams of noncentrosymmetric 2D systems and surface band
spectra obtained by focusing on the $k_{y}=0$ plane (top) and $k_{x}=0$ plane
(bottom) of the 3D BZ of the WSM (52). In each phase diagram, the top and
bottom density plots are the surface band spectra of the left and right edges,
respectively.

We obtain surface band spectrum with the real space lattice periodic
Hamiltonian as shown in Fig.20. First, we focus on the $k_{y}=0$ plane. The
$y$ direction is removed by fixing $k_{y}$, and the surface band spectrum is
calculated in a ribbon with a periodic $x$ direction and a finite $z$
direction. The bulk gap closing occurs at $\Delta_{T}^{0}=1.1$, which is the
point where the Weyl points in 3D BZ touch the $k_{y}=0$ plane. This is a
topological phase transition point between 2D NI and TI, and when
$\Delta_{T}^{0}>1.1$, we can see topological edge states that connect the
valence and conduction band spectra. Next we focus on the $k_{x}=0$ plane. In
this case, we consider a ribbon with a periodic $y$ direction and a finite $z$
direction. Also in this case, a topological phase transition occurs when the
Weyl points in 3D BZ touch the $k_{x}=0$ plane, $\Delta_{T}^{0}>1.0$. These
results correspond to looking a particular $\bm{k}$ slice of the surface state
of the 3D WSM. It should be noted that the energies of the surface states are
different in the left and right edges in both cases. The emergence of the
surface state with edge dependent energies is consistent with two natural
requirement in properties of surface states: (i) edge state should emerges
from a gap-closing point of the bulk band spectrum, (ii) a Dirac cone of a
time-reversal protected surface state should be located on a TRIM. At the same
time, this edge dependence indicates asymmetry of the left and right edges. If
the gap-closing of the bulk band spectrum occurs on a generic point in the
surface BZ, the bulk Hamiltonian is not allowed to have a symmetry that
exchange the left and right edges. Consequently, the position in 2D BZ where a
gap-closing point can appear is restricted by crystalline symmetries of the
bulk Hamiltonian, which is consistent with Ref.[66].

In conclusion of this section, we have shown the edge dependence of the
emergence of the surface states and its relation with symmetries of the bulk
Hamiltonian, and demonstrated it with an easy model.

##  Appendix G Twisted BHZ model

In this section, we propose a twisted Bernevig-Hughes-Zhang (BHZ) model as a
simple model to describe the essential behavior of twisted bilayer
Bi2(Te1-xSex)3.

Figure 21: Example of the twisted BHZ model ($\theta=0.80^{\circ},A=1,B=0.1$).
The stacking dependent $M$ is set to $(M_{AA},M_{AB},M_{AC})=(0.2,-0.3,0.4)$.
(a) Moiré band dispersion. (b)-(f) real space plot of the wave functions
(Upper layer, spin-up) at the $\Gamma$ point for VB1, VB2, VB3, VB6, and VB7.
(g) Magnified view of the VB6 and VB7 and Wilson loop spectra of them. (h)
Moiré edge state spectrum.

The BHZ model is a well-known effective model of a 2D topological insulator,
i.e. a quantum spin Hall insulator. The BHZ model is given with two orbitals
and two spin, and the Hamiltonian is typically written as

$\begin{split}&H_{\mathrm{BHZ}}(\bm{k})\\\
&=\left(M-Bk^{2}\right)\tau_{z}\sigma_{0}+Ak_{x}\tau_{x}\sigma_{z}+Ak_{y}\tau_{y}\sigma_{0}\\\
&=\left(\begin{array}[]{cccc}(M-Bk^{2})&Ak^{-}&&\\\ Ak^{-}&-(M-Bk^{2})&&\\\
&&(M-Bk^{2})&-Ak^{+}\\\ &&-Ak^{-}&-(M-Bk^{2})\end{array}\right),\end{split}$
(56)

where $k^{\pm}=k_{x}\pm ik_{y}$, $k^{2}=k_{x}^{2}+k_{y}^{2}$, and the basis is
ordered as $\ket{1,\uparrow}$, $\ket{2,\downarrow}$, $\ket{1,\downarrow}$,
$\ket{2,\uparrow}$. This model has the time-reversal symmetry ${\cal T}$ and a
rotation symmetry around z-axis $C_{\theta}$ (perpendicular to the 2D system),
and the operators are given as

${\cal T}=i\sigma_{y}{\cal
K}~{},~{}C_{\theta}=\left(\begin{array}[]{cccc}e^{i\frac{\theta}{2}}&&&\\\
&e^{-i\frac{\theta}{2}}&&\\\ &&e^{-i\frac{\theta}{2}}&\\\
&&&e^{i\frac{\theta}{2}}\\\ \end{array}\right)$ (57)

This model describes a topological insulator when $MB>0$ and a normal
insulator when $BM<0$.

By introducing a Moiré scale oscillation in $M$, we obtain a twisted BHZ
model. To make the model realistic, we transform the basis. We assume that one
of the two orbitals is located on the upper layer, and the other on the lower
layer, and spin is located on each orbital. This assumption restricts the
time-reversal and the rotation operators to be block-diagonal with intralayer
blocks. We also assume that the Moiré oscillating term $M$ should be an
interlayer element. To satisfy these assumptions, we transform the Hamiltonian
with a unitary matrix

$\begin{split}U&=\frac{1}{2\sqrt{2}}\left[(\tau_{0}+i\tau_{y})(\sigma_{0}-i\sigma_{x})+(\tau_{0}-i\tau_{y})(\sigma_{y}+\sigma_{z})\right]\\\
&=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cccc}1&0&-i&0\\\ 0&1&0&-i\\\
0&-i&0&1\\\ i&0&-1&0\\\ \end{array}\right),\end{split}$ (58)

and obtain a new representation of Hamiltonian

$\begin{split}&\tilde{H}_{\mathrm{BHZ}}(\bm{k})\\\
&=U^{\dagger}H_{\mathrm{BHZ}}(\bm{k})U\\\
&=\left(M-Bk^{2}\right)\tau_{z}\sigma_{y}+Ak_{x}\tau_{x}\sigma_{0}+Ak_{y}\tau_{y}\sigma_{0}\\\
&=\left(\begin{array}[]{cccc}&Ak^{-}&-i(M-Bk^{2})&\\\ Ak^{+}&&&i(M-Bk^{2})\\\
i(M-Bk^{2})&&&Ak^{-}\\\ &-i(M-Bk^{2})&Ak^{+}&\\\
\end{array}\right),\end{split}$ (59)

where the new basis is written as $\ket{1,\uparrow}-i\ket{2,\uparrow}$,
$\ket{2,\downarrow}+i\ket{1,\downarrow}$,
$i\ket{1,\uparrow}-\ket{2,\uparrow}$,
$i\ket{2,\downarrow}+\ket{1,\downarrow}$. In this representation, the time-
reversal and the rotation operators are written as

$\tilde{{\cal T}}=-\tau_{y}\sigma_{0}{\cal
K}~{},~{}\tilde{C}_{\theta}=\left(\begin{array}[]{cccc}e^{i\frac{\theta}{2}}&&&\\\
&e^{-i\frac{\theta}{2}}&&\\\ &&e^{i\frac{\theta}{2}}&\\\
&&&e^{-i\frac{\theta}{2}}\\\ \end{array}\right).$ (60)

By assigning $\tilde{H}_{\mathrm{BHZ}}(\bm{k})$ to
$h_{\bm{r}}^{\alpha\sigma,\beta^{\prime}\sigma^{\prime}}(\bm{k})$ in Eq. (29),
we obtain a Hamiltonian of the twisted BHZ model. By tuning $M$ for three
sampling points, we can design a Moiré superlattice system with topological
and normal insulator domains.

Note that the basis ordering makes a physical difference in a strict sense due
to the momentum twisting in Eqs. (26)(27)(28), although it makes no physical
difference in an untwisted case. However, the difference is negligible in the
small twist angle limit when we consider an effective model around the
$\Gamma$ point in the atomic BZ. We calculate with
$\tilde{H}_{\mathrm{BHZ}}(\bm{k})$, but it is confirmed that almost same
results can be obtained by calculating with $H_{\mathrm{BHZ}}(\bm{k})$.

Figure 21 shows an example of the twisted BHZ model with
$\theta=0.80^{\circ}$, $A=1$, $B=0.1$, and stacking-dependent $M$ that is set
to $(M_{AA},M_{AB},M_{AC})=(0.2,-0.3,0.4)$. Because $MB<0$ only in the AB-
stacking area, the topological insulator domain should appear around the AB-
stacking area, as is in twisted bilayer Bi2(Te1-xSex)3. It can be seen that
the expected ring-shape edge-state-originated states VB1-3 and their angular
momentum ordering are obtained as shown in Fig.21(b)-(d). There are also Moiré
topological bands VB6 and VB7 (Fig.21(g)), and a corresponding Moiré-scale
helical edge state is obtained as shown in Fig.21(h). Because the wave
functions of the VB6 (Fig.21(e)) and VB7 (Fig.21(f)) are a edge-state-
originated and a bulk-state-originated, respectively, the obtained Moiré-scale
helical edge state is a “edge state from edge state”. These results indicate
that the twisted BHZ model is a good simple model that reproduces the
properties of twisted bilayer Bi2(Te1-xSex)3. Furthermore, it is also shown
that the “edge state from edge state” is not specific for twisted bilayer
Bi2(Te1-xSex)3, but is a general phenomenon in similar Moiré systems.

## References

  * 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._ , Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018b).
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Origin of mott insulating behavior and superconductivity in twisted bilayer graphene, Phys. Rev. X 8, 031089 (2018).
  * 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, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (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, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Jiang  _et al._ [2019] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule, J. Mao, and E. Y. Andrei, Charge order and broken rotational symmetry in magic-angle twisted bilayer graphene, Nature 573, 91 (2019).
  * Polshyn  _et al._ [2019] H. Polshyn, M. Yankowitz, S. Chen, Y. Zhang, K. Watanabe, T. Taniguchi, C. R. Dean, and A. F. Young, Large linear-in-temperature resistivity in twisted bilayer graphene, Nature Physics 15, 1011 (2019).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, _et al._ , Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * 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._ , Maximized electron interactions at the magic angle in twisted bilayer graphene, 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._ , Electronic correlations in twisted bilayer graphene near the magic angle, Nature Physics 15, 1174 (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, Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene, Nature 572, 101 (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, Strange metal in magic-angle graphene with near planckian dissipation, Phys. Rev. Lett. 124, 076801 (2020).
  * Serlin  _et al._ [2020] M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Intrinsic quantized anomalous hall effect in a moiré heterostructure, Science 367, 900 (2020).
  * Stepanov  _et al._ [2020] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Untying the insulating and superconducting orders in magic-angle graphene, Nature 583, 375 (2020).
  * Saito  _et al._ [2020] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F. Young, Independent superconductors and correlated insulators in twisted bilayer graphene, Nature Physics 16, 926 (2020).
  * Hunt  _et al._ [2013] B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, _et al._ , Massive dirac fermions and hofstadter butterfly in a van der waals heterostructure, Science 340, 1427 (2013).
  * 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._ , Hofstadter’s butterfly and the fractal quantum hall effect in moiré superlattices, Nature 497, 598 (2013).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene, Phys. Rev. X 8, 031087 (2018).
  * Moon and Koshino [2012] P. Moon and M. Koshino, Energy spectrum and quantum hall effect in twisted bilayer graphene, Phys. Rev. B 85, 195458 (2012).
  * 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, Unraveling the intrinsic and robust nature of van hove singularities in twisted bilayer graphene by scanning tunneling microscopy and theoretical analysis, Phys. Rev. Lett. 109, 196802 (2012).
  * 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._ , Simulation of hubbard model physics in wse 2/ws 2 moiré superlattices, Nature 579, 353 (2020).
  * Regan  _et al._ [2020] E. C. Regan, D. Wang, C. Jin, M. I. B. Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, _et al._ , Mott and generalized wigner crystal states in wse 2/ws 2 moiré superlattices, Nature 579, 359 (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._ , Correlated electronic phases in twisted bilayer transition metal dichalcogenides, Nature materials 19, 861 (2020).
  * Wu  _et al._ [2019] F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, Topological insulators in twisted transition metal dichalcogenide homobilayers, Phys. Rev. Lett. 122, 086402 (2019).
  * Chu  _et al._ [2020] Z. Chu, E. C. Regan, X. Ma, D. Wang, Z. Xu, M. I. B. Utama, K. Yumigeta, M. Blei, K. Watanabe, T. Taniguchi, S. Tongay, F. Wang, and K. Lai, Nanoscale conductivity imaging of correlated electronic states in ${\mathrm{wse}}_{2}/{\mathrm{ws}}_{2}$ moiré superlattices, Phys. Rev. Lett. 125, 186803 (2020).
  * Xu  _et al._ [2020] Y. Xu, S. Liu, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, V. Elser, K. F. Mak, and J. Shan, Correlated insulating states at fractional fillings of moiré superlattices, Nature 587, 214 (2020).
  * Huang  _et al._ [2021] X. Huang, T. Wang, S. Miao, C. Wang, Z. Li, Z. Lian, T. Taniguchi, K. Watanabe, S. Okamoto, D. Xiao, _et al._ , Correlated insulating states at fractional fillings of the ws2/wse2 moiré lattice, Nature Physics 17, 715 (2021).
  * Padhi  _et al._ [2021] B. Padhi, R. Chitra, and P. W. Phillips, Generalized wigner crystallization in moiré materials, Phys. Rev. B 103, 125146 (2021).
  * Zhang  _et al._ [2021] Y. Zhang, T. Liu, and L. Fu, Electronic structures, charge transfer, and charge order in twisted transition metal dichalcogenide bilayers, Phys. Rev. B 103, 155142 (2021).
  * Liu  _et al._ [2021] E. Liu, T. Taniguchi, K. Watanabe, N. M. Gabor, Y.-T. Cui, and C. H. Lui, Excitonic and valley-polarization signatures of fractional correlated electronic phases in a ${\mathrm{wse}}_{2}/{\mathrm{ws}}_{2}$ moiré superlattice, Phys. Rev. Lett. 127, 037402 (2021).
  * Wu  _et al._ [2018a] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Hubbard model physics in transition metal dichalcogenide moiré bands, Phys. Rev. Lett. 121, 026402 (2018a).
  * Jin  _et al._ [2021] C. Jin, Z. Tao, T. Li, Y. Xu, Y. Tang, J. Zhu, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, _et al._ , Stripe phases in wse 2/ws 2 moiré superlattices, Nature Materials , 1 (2021).
  * Zhao  _et al._ [2020] X.-J. Zhao, Y. Yang, D.-B. Zhang, and S.-H. Wei, Formation of bloch flat bands in polar twisted bilayers without magic angles, Phys. Rev. Lett. 124, 086401 (2020).
  * Lian  _et al._ [2020] B. Lian, Z. Liu, Y. Zhang, and J. Wang, Flat chern band from twisted bilayer mnbi 2 te 4, Physical review letters 124, 126402 (2020).
  * Zhang  _et al._ [2009] H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Topological insulators in bi 2 se 3, bi 2 te 3 and sb 2 te 3 with a single dirac cone on the surface, Nature physics 5, 438 (2009).
  * Xia  _et al._ [2009] Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, _et al._ , Observation of a large-gap topological-insulator class with a single dirac cone on the surface, Nature physics 5, 398 (2009).
  * Chen  _et al._ [2009] Y. Chen, J. G. Analytis, J.-H. Chu, Z. Liu, S.-K. Mo, X.-L. Qi, H. Zhang, D. Lu, X. Dai, Z. Fang, _et al._ , Experimental realization of a three-dimensional topological insulator, bi2te3, science 325, 178 (2009).
  * Yazyev  _et al._ [2010] O. V. Yazyev, J. E. Moore, and S. G. Louie, Spin polarization and transport of surface states in the topological insulators ${\mathrm{bi}}_{2}{\mathrm{se}}_{3}$ and ${\mathrm{bi}}_{2}{\mathrm{te}}_{3}$ from first principles, Phys. Rev. Lett. 105, 266806 (2010).
  * Liu  _et al._ [2010] C.-X. Liu, H. Zhang, B. Yan, X.-L. Qi, T. Frauenheim, X. Dai, Z. Fang, and S.-C. Zhang, Oscillatory crossover from two-dimensional to three-dimensional topological insulators, Phys. Rev. B 81, 041307 (2010).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Jung  _et al._ [2014] J. Jung, A. Raoux, Z. Qiao, and A. H. MacDonald, Ab initio theory of moiré superlattice bands in layered two-dimensional materials, Physical Review B 89, 205414 (2014).
  * Momma and Izumi [2008] K. Momma and F. Izumi, Vesta: a three-dimensional visualization system for electronic and structural analysis, Journal of Applied Crystallography 41, 653 (2008).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Peng  _et al._ [2016] H. Peng, Z.-H. Yang, J. P. Perdew, and J. Sun, Versatile van der waals density functional based on a meta-generalized gradient approximation, Phys. Rev. X 6, 041005 (2016).
  * Becke [1993] A. D. Becke, Density‐functional thermochemistry. iii. the role of exact exchange, The Journal of Chemical Physics 98, 5648 (1993), https://doi.org/10.1063/1.464913 .
  * Pizzi  _et al._ [2020] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, Wannier90 as a community code: new features and applications, Journal of Physics: Condensed Matter 32, 165902 (2020).
  * Liu  _et al._ [2012] Y. Liu, G. Bian, T. Miller, M. Bissen, and T.-C. Chiang, Topological limit of ultrathin quasi-free-standing bi2te3 films grown on si(111), Phys. Rev. B 85, 195442 (2012).
  * Zhang  _et al._ [2010] Y. Zhang, K. He, C.-Z. Chang, C.-L. Song, L.-L. Wang, X. Chen, J.-F. Jia, Z. Fang, X. Dai, W.-Y. Shan, S.-Q. Shen, Q. Niu, X.-L. Qi, S.-C. Zhang, X.-C. Ma, and Q.-K. Xue, Crossover of the three-dimensional topological insulator bi2se3 to the two-dimensional limit, Nature Physics 6, 584 (2010).
  * Förster  _et al._ [2016] T. Förster, P. Krüger, and M. Rohlfing, $gw$ calculations for ${\mathrm{bi}}_{2}{\mathrm{te}}_{3}$ and ${\mathrm{sb}}_{2}{\mathrm{te}}_{3}$ thin films: Electronic and topological properties, Phys. Rev. B 93, 205442 (2016).
  * Wu  _et al._ [2018b] Q. Wu, S. Zhang, H.-F. Song, M. Troyer, and A. A. Soluyanov, Wanniertools: An open-source software package for novel topological materials, Computer Physics Communications 224, 405 (2018b).
  * Murakami  _et al._ [2007] S. Murakami, S. Iso, Y. Avishai, M. Onoda, and N. Nagaosa, Tuning phase transition between quantum spin hall and ordinary insulating phases, Phys. Rev. B 76, 205304 (2007).
  * Murakami [2007] S. Murakami, Phase transition between the quantum spin hall and insulator phases in 3d: emergence of a topological gapless phase, New Journal of Physics 9, 356 (2007).
  * Chalker and Coddington [1988] J. Chalker and P. Coddington, Percolation, quantum tunnelling and the integer hall effect, Journal of Physics C: Solid State Physics 21, 2665 (1988).
  * Obuse  _et al._ [2007] H. Obuse, A. Furusaki, S. Ryu, and C. Mudry, Two-dimensional spin-filtered chiral network model for the ${\mathbb{z}}_{2}$ quantum spin-hall effect, Phys. Rev. B 76, 075301 (2007).
  * Obuse  _et al._ [2008] H. Obuse, A. Furusaki, S. Ryu, and C. Mudry, Boundary criticality at the anderson transition between a metal and a quantum spin hall insulator in two dimensions, Phys. Rev. B 78, 115301 (2008).
  * Obuse  _et al._ [2014] H. Obuse, S. Ryu, A. Furusaki, and C. Mudry, Spin-directed network model for the surface states of weak three-dimensional ${\mathbb{z}}_{2}$ topological insulators, Phys. Rev. B 89, 155315 (2014).
  * Ryu  _et al._ [2010] S. Ryu, C. Mudry, H. Obuse, and A. Furusaki, The network model for the quantum spin hall effect: two-dimensional dirac fermions, topological quantum numbers and corner multifractality, New Journal of Physics 12, 065005 (2010).
  * Sanchez-Yamagishi  _et al._ [2017] J. D. Sanchez-Yamagishi, J. Y. Luo, A. F. Young, B. M. Hunt, K. Watanabe, T. Taniguchi, R. C. Ashoori, and P. Jarillo-Herrero, Helical edge states and fractional quantum hall effect in a graphene electron–hole bilayer, Nature nanotechnology 12, 118 (2017).
  * Tong  _et al._ [2017] Q. Tong, H. Yu, Q. Zhu, Y. Wang, X. Xu, and W. Yao, Topological mosaics in moiré superlattices of van der waals heterobilayers, Nature Physics 13, 356 (2017).
  * Rickhaus  _et al._ [2018] P. Rickhaus, J. Wallbank, S. Slizovskiy, R. Pisoni, H. Overweg, Y. Lee, M. Eich, M.-H. Liu, K. Watanabe, T. Taniguchi, _et al._ , Transport through a network of topological channels in twisted bilayer graphene, Nano letters 18, 6725 (2018).
  * San-Jose and Prada [2013] P. San-Jose and E. Prada, Helical networks in twisted bilayer graphene under interlayer bias, Phys. Rev. B 88, 121408 (2013).
  * Efimkin and MacDonald [2018] D. K. Efimkin and A. H. MacDonald, Helical network model for twisted bilayer graphene, Phys. Rev. B 98, 035404 (2018).
  * Huang  _et al._ [2018] S. Huang, K. Kim, D. K. Efimkin, T. Lovorn, T. Taniguchi, K. Watanabe, A. H. MacDonald, E. Tutuc, and B. J. LeRoy, Topologically protected helical states in minimally twisted bilayer graphene, Phys. Rev. Lett. 121, 037702 (2018).
  * Ramires and Lado [2018] A. Ramires and J. L. Lado, Electrically tunable gauge fields in tiny-angle twisted bilayer graphene, Phys. Rev. Lett. 121, 146801 (2018).
  * Kariyado and Vishwanath [2019] T. Kariyado and A. Vishwanath, Flat band in twisted bilayer bravais lattices, Phys. Rev. Research 1, 033076 (2019).
  * Yu and Liu [2020] J. Yu and C.-X. Liu, Piezoelectricity and topological quantum phase transitions in two-dimensional spin-orbit coupled crystals with time-reversal symmetry, Nature communications 11, 1 (2020).
  * Burkov and Balents [2011] A. A. Burkov and L. Balents, Weyl semimetal in a topological insulator multilayer, Phys. Rev. Lett. 107, 127205 (2011).
  * Halász and Balents [2012] G. B. Halász and L. Balents, Time-reversal invariant realization of the weyl semimetal phase, Phys. Rev. B 85, 035103 (2012).