Geometric origins of topological insulation in twisted layered semiconductors

Hao Tang,1 Stephen Carr,2 and Efthimios Kaxiras3, 4

1
Department of Physics, Peking University, Beijing, 100871, P.R.China
2
Brown Theoretical Physics Center and Department of Physics,
Brown University, Providence, Rhode Island 02912-1843, USA
3
Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA
4
John A. Paulson School of Engineering and Applied Sciences,
Harvard University, Cambridge, Massachusetts 02138, USA
(Dated: January 14, 2021)
Twisted bilayers of two-dimensional (2D) materials are proving a fertile ground for investigating
arXiv:2101.04867v1 [cond-mat.mtrl-sci] 13 Jan 2021

strongly correlated electron phases. This is because the moiré pattern introduced by the relative
twist between layers introduces long-wavelength effective potentials which lead to electron localiza-
tion. Here, we develop a generalized continuum model for the electronic structure of moiré patterns,
based on first-principles calculations and tailored to capture the physics of twisted bilayer 2D semi-
conductors. We apply this model to a database of eighteen 2D crystals covering a range of atomic
relaxation and electronic structure features. Many of these materials host topologically insulating
(TI) moiré bands in a certain range of twist angles, which originate from the competition between
triangular and hexagonal moiré patterns, tuned by the twist angle. The topological phases occur in
the same range as the maximally flat moiré bands.

Quantum materials, engineered by creative manipula- It is difficult to model twistronic systems at small twist
tion of the features of conventional crystals, offer new angles (θ ' 1◦ ) using first-principles calculations, be-
possibilities for breakthroughs in understanding electron cause the number of atoms in the MSL scales as θ−2 .
correlations and superconductivity. What is intrigu- To overcome this limitation, continuum models with a
ing is that these phenomena emerge at a length scale low-energy effective Hamiltonian based on first-principles
much larger than the underlying crystal lattice constant (density-functional theory, DFT) calculations were devel-
(by a factor of 10 to 1000), through features due to oped for electronic structures of TBG; this approach can
strain patterns, controlled by geometric constraints. A accurately describe flat bands and magic angles [9, 15–
new platform for such studies are systems comprising 17]. Although continuum models have also been applied
few-layer two-dimensional (2D) crystals, like twisted bi- to the twisted bilayer semiconductors [14, 18], they have
layers of graphene or transition metal dichalcogenides yet to include the effect of atomic relaxations which play
(TMDCs) [1–5]. The slight lattice mismatch between two an important role at small angles [11, 12, 19, 20].
layers of a 2D material at a relative twist angle results in
In this letter, we present results from a DFT-based
a moiré superlattice (MSL) and a long-wavelength peri-
generalized continuum method designed specifically for
odic modulation of the effective electronic potential [6, 7].
twisted bilayer semiconductors. The computed electronic
In a narrow range of the twist angle, the moiré poten-
structures are consistent with full-DFT results [2, 11, 12]
tials act as confining wells for the electrons of the con-
but only require a relatively inexpensive set of bilayer cal-
stituent monolayers, causing isolated flat bands and lo-
culations. These calculations involve systems with only
calized wave functions near the Fermi surface [8, 9].
a handful of atoms per unit cell, in contrast to the many
thousands of atoms necessary in the calculation of a full
In the moiré flat bands, the kinetic energy is heavily
MSL. We derive a database of relaxation and the corre-
suppressed and electronic interactions play a dominant
sponding coefficients of tight-binding electronic structure
role, with the intensity of the interactions controlled by
hamiltonians for eighteen materials with various lattice
the twist angle; this effect has been dubbed “twistron-
symmetries and band edge momenta. These coefficients
ics” [10]. Compared to twisted bilayer graphene (TBG),
do not capture the material’s entire band structure, but
the twisted bilayer semiconductors can host flat bands in
rather focus on the details of the parabolic band edges.
a large range of twist angles [11, 12] instead of at pre-
Each layer contributes one band to the full twisted bi-
cisely a magic angle [1]. This makes it possible to over-
layer model, and these bands are coupled through a set
come some experimental challenges in twisted bilayers of
of stacking-dependent electronic interactions.
semiconductors; thus, the twist angle becomes an addi-
tional degree of freedom for fine-tuning other physical ef- In Fig. 1 we provide an overview of the different phases
fects in the strongly correlated regime [13]. Intriguingly, and their geometric origins, as the twist angle is changed
topological insulator (TI) moiré bands were predicted in in a moiré bilayer (Fig. 1a). The full interaction be-
a twisted bilayer of MoTe2 [14]; this work introduced tween the bands can be decomposed into two comple-
a possible candidate for observing concurrent correlated mentary parts. The first describes the tendency for elec-
and TI phases in the same material. trons in one layer to tunnel to the other, as shown in
 2

(a) L2 (b) Interlayer Tunneling ΔT potentials naturally lead to quantum harmonic oscillator
AA
(QHO) states at small twist angles [21–23], whose pres-
ence underlies the appearance of the moiré flat bands
and whose competition explains the large number of TI
(c) Electrostatic Repulsion Δt/Δb phases possible in twisted semiconductors.
We briefly introduce the methodology here. To cap-
BA ture the stacking-dependent electronic and atomic de-
AB tails, we perform DFT calculations on aligned bilayers
L1
over a grid sampling all possible interlayer displacements.
Using the effective mass approximation, we treat the dis-
(d) Dominant interaction (e) QHO-TR QHO-HC persion around the monolayer band extrema as a kinetic
1 TI
Band Energy
unrelaxed 3 2
TR energy term in a continuum Hamiltonian. Including the
ΔT

SC/MI relaxed 2 TI
2 4
three stacking-dependent potentials described in the in-
TI 3 troduction, we obtain the effective Hamiltonian:
Δt/Δb

Metallic Bands
HC Δt/Δb !
ΔT − ~ (∇−ik
2 2
0)
0° θ 5° Dominant Potential 2m∗ + ∆t (r) ∆∗T (r)
H=
− ~ (∇−ik
2 2
0)
∆T (r) 2m∗ + ∆b (r)
FIG. 1. (a) A moiré pattern of a twisted TMDC with areas (1)
of aligned (AA) stacking and eclipsed (AB or BA) stacking where m∗ is the effective mass, ∆t,b are the electrostatic
highlighted. (b) At AA stacking, the local electronic Hamil-
potentials for the top and bottom layer, and ∆T is the
tonian is described by interlayer tunneling between the layers
(∆T ), forming a triangular (TR) lattice. (c) At AB or BA interlayer tunneling strength. At small twist angles we
stacking, the electronic Hamiltonian is described by interlayer also include the all-important relaxation effects by mini-
electrostatic potentials (∆t, ∆b), forming a honeycomb (HC) mizing the total mechanical energy of the moiré patterns
lattice. (d) Dependence of the dominant moiré interaction on [20]. The local electronic structures can then be derived
twist angle, with atomic relaxation of the superlattice (red from the DFT through the expansion:
line) or without it (blue line): evolution of the lattice charac-
ter from TR to HC with decreasing θ induces TI states and ~2 (k − k0 )2
flat bands at small angles, leading to Mott insulator (MI) and E (±) (r, k) = E (±) (θ × r + 2u) + (2)
2m∗
superconducting (SC) behavior. (e) The TI phase is caused by
band hybridization as the electrons transition between quan- To determine the ∆t/b , ∆T (r) from the local electronic
tum harmonic oscillator (QHO) states on a TR to a HC lat- structure, the Bloch wave functions at band extrema are
tice.
extracted from DFT calculations to assess the layer po-
larization of each band. The moiré bands can then be
Fig. 1b – labeled ∆T for “Tunneling”. The tunneling co- calculated by diagonalizing the matrix form of the Hamil-
efficients are strongest at AA (aligned) stacking regions tonian for a truncated basis set in k-space (see SM for
and form a triangular (TR) lattice across the MSL. The details).
second contribution captures the stacking dependence of The aligned (AA) stacking configuration has higher en-
the monolayer bands’ on-site energies, which depends on ergy than the partially eclipsed (AB/BA) configuration
the electrostatic potential from the opposite layer, shown for TMDCs and hBN homo-bilayers (θ ' 0◦ ). Conse-
in Fig. 1c) – labeled ∆t and ∆b for “top/bottom” layer. quently, relaxation tends to reduce the in-plane area of
The electrostatic potentials have maxima at the AB and AA stacking region and to increase that of the AB/BA
BA stacking regions, forming a honeycomb (HC) lattice. stacking region, thus minimizing the total energy [16, 20].
In the majority of materials studied here, the in- Upon relaxation, the large values of ∆t,b , at AB and
plane atomic relaxation of the twisted bilayers enhances BA stacking, expand to cover a larger area, while the
the effects of the electrostatic potential fluctuations and peak regions of ∆T , at AA stacking, shrink, with relaxed
reduces that of the interlayer tunneling, especially at structures showing stronger electrostatic potential effects
small twist angles. Electron localization transitions, from and weaker tunneling effects, which causes a clear angle-
the tunneling-dominated TR lattice to the electrostatic- dependent transition of the moiré electronic structure.
dominated HC lattice, occur at low twist angle in various In Fig. 2a,b we present the angle-dependent bands and
materials, as illustrated in Fig. 1d. In the low-angle re- real-space localization of twisted bilayer MoTe2 as a rep-
gion, ultra-flat moiré bands are predicted, likely to host resentative case. The top-most valence bands for both
superconducting (SC) and Mott insulating (MI) phases. the relaxed and unrelaxed moiré systems have a band-
The topologically insulating (TI) moiré bands appear for width less than 10 meV when the twist angle is below
intermediate values of the twist angle, during the transi- 2.5◦ . Using the Coulomb repulsion energy (U ) in TBG
tion between the two types of lattice geometries for the and twisted hBN as a guide [4, 11], this small kinetic
interlayer electronic interactions (Fig. 1e). The electronic energy implies that strongly correlated states could exist
 3

(a) MoTe2 (c) Unrelaxed: TR QHO erate), and double excitation states (3-fold degenerate).
Unrelaxed (0, 0) (2, 0) LDOS
0 n = 8
n = 1
After relaxation, the competition between the two
-20 n = 2 types of lattices is greatly altered, as shown in Fig. 2d.
ave[En] (meV)

10 meV
(1, 0) (1, 1) 4 The two top bands in the low-angle case (θ = 0.7◦ ) lo-
-40
1 calize in the AB/BA HC potential wells to from 2-fold
-60 10 0
(0, 1) (0, 2)
degenerate ground states of the HC QHO model, while
-80 2 the third band is one of the 4-fold degenerate first ex-
cited states. In contrast, the localization of the upper-
(b) Relaxed (d) Relaxed: TR vs HC most band gradually changes to the AA stacking in the
0 0.70 20 LDOS
3 high-angle case (θ = 2◦ ) with the second band localized
ave[En] (meV)

-20 1 in the AB/BA region. This indicates that the energy

2
level of the HC QHO model and the triangular QHO
-40 0.70
2 2 1 model cross as the angle changes, leading to a reordering
-60 4 of the bands’ localization. Therefore, the transition be-
0
-80 3
tween the HC and TR electronic states can be controlled
y by varying the twist angle in the presence of atomic re-
θ (o) x laxations. The intermediate twist angles correspond to
competing real-space distribution and band reordering,
FIG. 2. Angle-dependent band structures and real-space lo- making these materials excellent candidates for hosting
calization. (a, b) The band average energy (ave[En ]), band- non-trivial topological properties.
width (color saturation, deeper color for flatter bands), and To study the topological properties associated with the
localization (red for AB/BA stacking and blue for AA stack- TR-HC transition, we calculate the Chern numbers of the
ing) of the three top moiré valence bands as a function of the
moiré bands for a generic twisted bilayer of a hexagonal
twist angle θ for twisted bilayers of MoTe2 . A fitted QHO
model is given by the dashed lines, and the first twist angles semiconductor, and use these results to generate a TI
with bandwidths less than 10 meV are denoted with green phase diagram. The condition for the existence of topo-
stars. The 1◦ band structures are shown in the insets, with logical bands is represented by the electrostatic potential
the first two QHO levels highlighted in color. (c) LDOS of and tunneling coefficients in Fig. 3a for a generic sys-
the top six moiré bands of unrelaxed bilayer MoTe2 twisted at tem, which are approximated by their first order Fourier
0.5◦ . The quantum numbers in the QHO model are denoted coefficients labeled V and w, respectively. In this simpli-
as n = (nx , ny ). (d) LDOS of the moiré bands for a 0.7◦ and a
fied model, the topological properties can be completely
2◦ twisted bilayer of MoTe2 . The average LDOS in the MSL
are normalized to one in both (c) and (d). described by just two parameters, α = V m∗ a2 /θ2 and
β = wm∗ a2 /θ2 , with the Hamiltonian

for any twist angle below a certain critical value in some θ2  

TMDC materials [2, 24]. The relaxation effects drive the H ≈ K̂ + V̂ + Ŵ = ∗ 2
K̂0 + αVˆ0 + β Ŵ0 (3)
m a
small-angle (< 2◦ ) valence moiré bands from TR-type
to HC-type (Fig. 2). In the unrelaxed case, the band where a is the lattice parameter, K̂, V̂ , Ŵ are the ki-
structure at 1◦ shows a single uppermost flat band and netic, potential, and tunneling energy, respectively, and
a pair of flat bands under it, consistent with a AA stack- K̂0 , Vˆ0 , Ŵ0 are the material-independent quantities. The
ing QHO model [21, 22]. The energy levels of the 2D topologically trivial TR and HC phases appear in the w-
TR QHO model are given by En = −~ω(n1 + n2 + 1), dominated and V -dominated regions of the phase space
ω = ωθ θ, and the first three levels −~ω,−2~ω, and −3~ω of Fig. 3a, respectively. Interesting band structures occur
are accurately represented by the computed band aver- in the intermediate region, which includes two types of
ages for θ < 1◦ [21]. In contrast, a pair of top moiré topological non-trivial valence moiré bands. The Chern
bands and the four lower bands at 0.7◦ correspond to the numbers of the three top-most moiré bands are (1, -1,
first (n = 0) and second (n = 1) QHO states, respectively 0) and (1, 1, -2) for X-type and Y-type topological in-
of the HC lattice [22]. sulators, respectively. Transitions between topological
The real space localization of the tunneling-dominated and trivial insulators for the uppermost band will occur
(TR) QHO state is evident in the local density of states on the boundary between the Y-type topological phase
(LDOS) of the moiré bands for the unrelaxed condi- and TR phase where the top two bands overlap and then
tion, as shown in Fig. 2c. The electronic states local- reopen, hosting protected edge modes at the transition.
ize around the AA stacking center and show s, p, and Similar topological transitions between the second and
d -orbital distribution for the first three energy levels, re- the third bands appears in the boundary between the X-
spectively. The six moiré bands shown correspond to the type and Y-type topological phases. The transition on
QHO ground state, single excitation states (2-fold degen- the boundary of the X-type topological phase and the
 4

FIG. 3. Topological phase diagram for moiré bands. (a) Phase of the moiré bands as a function of the intensity of potential
and tunneling fluctuations (1st Fourier components). The light (dark) green areas represent X (Y) type TI phase, while the
orange and light blue areas represent topologically trivial HC and TR QHO states, respectively. The twist-angle dependence
of the relaxed band structures is shown by the solid lines for selected materials. The band structures and Chern numbers for
unrelaxed and relaxed moiré bands in 1◦ twisted bilayer WSe2 are compared. (b) Berry curvature of the uppermost band for
X-type (1◦ ) and Y-type (1.5◦ ) twisted bilayer WSe2 in the MBZ.

HC phase is a TI-to-semimetal transition. likely to induce a spin or valley swapping.

In Fig. 3a we show the curves in the phase diagram for
4 TMDC’s, indicating transitions between phases with TABLE I. θHC/X , θX/Y , and θY/TR are the critical twist an-
increasing twist angle (without relaxation these curves gles (in degrees) for transitions of the top bands between the
would be straight lines like the one shown for WSe2 ). indicated phases. Egap is the maximum band gap between
For comparison, we include the curve for hBN, which the top and lower bands (in meV), and ωtop is the minimum
remains in one phase (TR) as it shows very weak relax- band width of the top bands in the topologically non-trivial
ation. An explicit example of the topological transitions regime (in meV); θgap and θtop are the angles where these
extrema occur.
is also shown Fig. 3a for WSe2 , including the unrelaxed
case. Material θHC/X θX/Y θY/TR Egap (θgap ) ωtop (θtop )
MoSe2 0.9 1.4 2.4 2.98 (1.67◦ ) 0.18 (1.21◦ )
The band structures and Chern numbers of the bands MoTe2 0.7 1.1 1.5 1.43 (1.20◦ ) 0.30 (1.10◦ )
confirm that the topological transition is primarily driven
WSe2 0.8 1.2 2.1 3.24 (1.40◦ ) 0.19 (1.15◦ )
by atomic relaxation. In the unrelaxed 1◦ twisted WSe2 ,
WTe2 0.7 1.4 2.4 3.70 (1.68◦ ) 0.18 (1.31◦ )
although the uppermost valence band is topologically
trivial, the second valence band has a nonzero Chern
number. The topological features of these lower valence In Fig. 3b we show the Berry curvature of the upper-
bands are not reflected by the phase diagram, but we list most valence band for two values of the twist angle. In
them in our database for each material. In realistic ex- the HC phase in the small-angle region, the first transi-
perimental conditions, atomic-scale defects, bending, and tion appears at θHC/X = 1◦ where the two top bands sep-
local strain introduced during fabrication cause twisted arate at the MBZ’s K points, indicated by the Berry cur-
bilayers to include different twist angles in different areas vature’s concentration there. With increasing twist an-
within a single sample [25]. Therefore, different domains gle, the Berry curvature of the top band gradually trans-
of twist angle occur and introduce both topological and fers to the Γ point accompanied by an X/Y transition at
trivial bands. Topological edge states can appear along θX/Y = 1.5◦ . During this process, the topologically non-
domain boundaries in the twist angle, according to the trivial top-most band takes on its minimal bandwith ,
critical values in Tab. I. These edge states would not be ωtop ' 0.19 meV, and its maximal topological band gap
defined by a sharp change from one crystal to another, Egap ' 3.24 meV. At larger twist angle, the Berry cur-
or from material to vacuum, but rather by slow varia- vature concentrates at the Γ point during the transition
tion in the twist angle. For this reason, these “internal” to the TR phase, where the two top-most bands merge.
protected edge states would be excellent candidates for The numerical results for the critical angles, maximal
observing spin or valley-polarized states, as the only dis- gap, and minimal band width for the TMDCs with topo-
order comes from twist-angle variations, which are un- logical valence bands are presented in Table 1. We have
 5

also verified that the bands and real space localization of C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim,
MoS2 [12], WSe2 [2], and hBN [11] from our model are et al., Nat. Mater. , 1 (2020).
consistent with previous full DFT calculations. Twisted [3] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watan-
bilayers of WTe2 show the highest gap and one of the low- abe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean,
Science 363, 1059 (2019).
est bandwidths, making this material the best candidate [4] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken,
for experimental study. J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe,
Including spin degrees of freedom, the spin-dependent T. Taniguchi, E. Kaxiras, et al., Nature 556, 80 (2018).
moiré Hamiltonian decomposes into two copies [14]. At [5] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao,
the K valleys of the aligned bilayer, the two copies will R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von
Oppen, A. Stern, et al., Nature 582, 203 (2020).
be split by the spin-orbit coupling term, ∆SOC , which [6] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H.
for most TMDCs is on the scale of 100 meV, and the Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
opposite spin-ordering will occur at the K 0 valley. The [7] C.-H. Park, L. Yang, Y.-W. Son, M. L. Cohen, and S. G.
Γ and M points tend to have very weak spin-splitting, Louie, Nat. Phys. 4, 213 (2008).
leading to two copies of the spin-independent Hamilto- [8] R. Bistritzer and A. H. MacDonald, PNAS 108, 12233
nian, as is the case in continuum treatments of twisted (2011).
bilayer graphene [8]). However, as the two nonequiva- [9] E. S. Morell, J. Correa, P. Vargas, M. Pacheco, and
Z. Barticevic, Phys. Rev. B 82, 121407 (2010).
lent K valleys are related by time reversal symmetry, the [10] S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin,
spin-up bands at K valley and spin-down bands at K 0 val- and E. Kaxiras, Phys. Rev. B 95, 075420 (2017).
ley have opposite Berry curvature and Chern numbers. [11] L. Xian, D. M. Kennes, N. Tancogne-Dejean,
Topologically non-trivial uppermost bands with opposite M. Altarelli, and A. Rubio, Nano Lett. 19, 4934
spin and Chern number give rise to the helical edge states (2019).
protected by the TR symmetries at the boundaries, and [12] M. H. Naik and M. Jain, Phys. Rev. Lett. 121, 266401
could lead to observable quantum spin Hall (QSH) ef- (2018).
[13] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian,
fects [26]. M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi,
The origin of the TI phases has been ascribed [14] to J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy, Na-
a skyrmionic texture in the moiré ∆ potentials, which ture 572, 95 (2019).
is an alternate interpretation of the TR/HC competition [14] F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. MacDon-
we presented here but captures the same key features. ald, Phys. Rev. Lett. 122, 086402 (2019).
[15] S. Carr, S. Fang, Z. Zhu, and E. Kaxiras, Phys. Rev. Res.
As the flatness of the bands causes stronger correlation 1, 013001 (2019).
effects, the topological phases in twisted bilayer semi- [16] J. Jung, A. Raoux, Z. Qiao, and A. H. MacDonald, Phys.
conductors could exhibit the combination of TI and su- Rev. B 89, 205414 (2014).
perconductivity, which has been the subject of an intense [17] S. Fang and E. Kaxiras, Phys. Rev. B 93, 235153 (2016).
decades-long search for fractional statistics and Majorana [18] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Phys.
fermions [26]. Rev. Lett. 121, 026402 (2018).
[19] N. N. T. Nam and M. Koshino, Phys. Rev. B 96, 075311
The generality of the relaxation-induced TI phases (2017).
is partially explained by the connection between the [20] S. Carr, D. Massatt, S. B. Torrisi, P. Cazeaux, M. Luskin,
stacking-dependence of the electronic structure and that and E. Kaxiras, Phys. Rev. B 98, 224102 (2018).
of the ground state energy: the lowest stacking configu- [21] S. Carr, D. Massatt, M. Luskin, and E. Kaxiras, Phys.
ration of these materials seem to correspond to maxima Rev. Res. 2, 033162 (2020).
of the electrostatic interaction and minima of the tun- [22] M. Angeli and A. H. MacDonald, arXiv e-prints (2020),
arXiv:2008.01735 [cond-mat.str-el].
neling strength. We note that this trend might not be
[23] M. H. Naik, S. Kundu, I. Maity, and M. Jain, Phys. Rev.
true across all 2D materials, and finding exceptions to B 102, 075413 (2020).
the rule could be a valuable endeavor. [24] Z. Zhang, Y. Wang, K. Watanabe, T. Taniguchi,
We thank Ziyan Zhu and Daniel Larson for useful dis- K. Ueno, E. Tutuc, and B. J. LeRoy, Nature Physics 16,
cussions. This work was supported by the National Sci- 1093 (2020).
[25] A. Uri, S. Grover, Y. Cao, J. A. Crosse, K. Bagani,
ence Foundation under grant No. OIA-1921199. The
D. Rodan-Legrain, Y. Myasoedov, K. Watanabe,
calculations in this work were performed in part on the T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero,
FAS Research Computing cluster of Harvard University. and E. Zeldov, Nature 581, 47 (2020).
[26] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057
(2011).

[1] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi,

E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43
(2018).
[2] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes,
 Geometric Origins of topological insulation in twisted semiconductors
(supplementary materials)

Hao Tang,1 Stephen Carr,2 and Efthimios Kaxiras3, 4

1
Department of Physics, Peking University, Beijing, 100871, PRC
2
Brown Theoretical Physics Center and Department of Physics,
Brown University, Providence, Rhode Island 02912-1843, USA
3
Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA
4
John A. Paulson School of Engineering and Applied Sciences,
Harvard University, Cambridge, Massachusetts 02138, USA
(Dated: January 14, 2021)
arXiv:2101.04867v1 [cond-mat.mtrl-sci] 13 Jan 2021

METHODOLOGY (a) (b)

∆t b
The DFT calculations of commensurate bilayers use ∆b
K-edge
the projector-augmented wave (PAW) basis sets with cut- Γ-edge
Kt 2|∆T|
off energy of 400 eV implemented by the Vienna ab ini-
Kb
Μsc Γsc
tio simulation package (VASP) [1, 2]. To overcome the
band gap problem of the generalized gradient approxima-
tion (GGA), we perform a meta-GGA calculation with Mb Mt
Strongly constrained and appropriately normed (SCAN) M-edge
semilocal density functional [3, 4], with the vdW-DF
functional rVV10 to more accurately capture the van (c) GSFE (meV) (d) Moire Potential (meV)
der Waals interaction [5]. A vacuum layer of 15 Å in AA
the c (z) direction is set and the dipole corrections are 20
∆t
switched on to cancel any residual inter-slab interactions. 10
The k-point mesh is sampled by the Monkhorst package BA 0
with separation of 0.2 Å [1]. All the electronic iterations -10
are converged to 10−5 eV and the atomic relaxations are AA AA ∆b
-20
converged to 10−4 eV. AB 24
The mechanical properties (generalized stacking fault 18
energy (GSFE)) [6] and electronic structure properties |∆T| 12
are calculated for displacement d = u1 a1 = u2 a2 on a y
AA 6
9 × 9 grid in the unit cell for each bilayer material, where
a1,2 is the unit cell lattice vectors. For each displacement, x 0

the in-plane position of the atoms is fixed while the layer

FIG. 1. Electronic model for the twisted bilayers. (a) Bril-
distance is relaxed. The typical variation ranges of z co- louin zones for a twisted bilayer of hexagonal lattices. The
ordinates of TMDCs for different stacking are about 0.4 black hexagons show the MBZ expanded around selected high
Å. The strain modulus, GSFE, band edge momenta, en- symmetry points. (b) Formation of moiré flat bands (black)
ergy level, effective mass, the corresponding Kohn-Sham from top (red) and bottom (blue) single-layer bands near the
wave functions, and their projection weight to the top band (c) GSFE landscape (color) and atomic displacements
and bottom layers, are extracted for each atomic dis- (arrows) for a 2◦ twisted bilayer of WSe2 . (d) ∆t,b and |∆T |
placement independently, and then applied in the con- values from DFT without (left) and with (right) atomic relax-
ation; stars represent AA (red), AB (green), and BA (yellow)
tinuum model. stacking.
The continuum model for the moiré superlattice (MSL)
of twisted bilayers is described as follow. The transla-
tion symmetry of the moiré superlattice defines the moiré
Brillouin zone (MBZ), as illustrated in Fig. 1a. The moiré
bands closest to the band gap can be derived from the
band extrema in the original Brillouin zone, which are constructed by these potentials as Eq. (1) in the main
typically located at high symmetry k-points. In the local body.
electronic structures, the bands of individual layers are
shifted by the electrostatic potential ∆t,b and split by the At small twist angles we also include the all-important
magnitude of the interlayer tunneling 2|∆T |, as shown in relaxation effects by minimizing the total mechanical en-
Fig. 1b. The electronic Hamiltonian in the MSL is then ergy of the moiré patterns [6]. This energy is a summa-
 2

tion of the (GSFE) and intra-layer strain energy: (c+ ∗ +

t ) cb
∆T = (E + − E−) (7)
XZ 1 Z |cb | + |c+
+ 2
t |2

E= j Cj j d2 r + VGSFE (d(r))d2 r (1a)

2
j=t,b The c coefficients can be derived from DFT calcula-
d(r) = d0 (r) + ut (r) − ub (r) (1b) tions by two methods. In the first method, one projects
the Bloch wave functions onto atomic orbitals in the
where ut/b (r) is the in-plane displacement of the top/bottom layer, so that the value of |c± t/b | is deter-
top/bottom layer, and , C, VGSFE , d0 are the strain
mined by the projection weight of ψ±,k to each of the
tensor, strain modulus matrix, GSFE density, and dis-
layers. In this case, the phase factor in the tunneling in-
placement under rigid twist, respectively. The effect of
tegral must be derived from the symmetry, like F = eiφ
atomic relaxations can be easily incorporated into the
with φ = arg (1 + eiG1 ·d + eiG2 ·d ) for a K band extrema,
electronic Hamiltonian. In the unrelaxed case, the DFT-
or F = 1 for a Γ point band extrema in a hexagonal lat-
extracted values of the stacking-dependent potentials are
tice [7, 8]. The second method is more straightforward,
unwrapped onto the moiré pattern with a straightfor-
but computationally more intensive. The c coefficients
ward linear mapping. To include relaxations we add one
are calculated by taking the inner product of the bilayer
step after the linear mapping which updates the spatial
Bloch wave functions with a single (top/bottom) layer
dependence of the potentials to accurately reflect the
projected Bloch function (derived from a separate ML
relaxed configurations of the bilayer. As the stacking
DFT calculation and shifted to the corresponding posi-
displacement is a function of the real space coordinate,
tion of the bilayer) at the same band edge momenta:
d(r) = θ × r + 2u (u = ut = −ub , considering the
layer-exchange symmetry for the in-plane relaxations), c± t ± ± b ±
t = hφk (r)|ψk (r)i , cb = hφk (r)|ψk (r)i (8)
we derive the local band energy at arbitrary position r
expanded in the effective mass approximation: Although the derived coefficients c still have a phase fac-
tor uncertainty due to the gauge-freedom of the Bloch
~2 (k − k0 )2 phase, that factor is the same for c+ +
E (±) (r, k) = E (±) (θ × r + 2u) + (2) b and ct and will be
2m∗ canceled when calculating ∆T .
where the band extrema energy E ± and effective mass After deriving the Hamiltonian from the DFT, the next
m∗ are determined by the DFT. step is to solve the moiré electronic structures. The Bloch
To determine the ∆t/b , ∆T (r) from the local elec- wave functions across the moiré superlattice are
tronic structure, the Bloch wave functions at band ex-   X
unk,t (r)
trema are extracted from DFT calculations to assess ψnk (r) = eik·r , unk,t/b (r) = m
γnk,t/b eiGm ·r
unk,b (r)
the layer polarization of each band. At the band ex- m
trema of bilayer structures, any non-degenerate mono- (9)
layer band edge splits into the bonding (+) and anti- where m runs over all reciprocal lattice vectors (Gm ) of
bonding (−) combinations, with energies E (±) and wave- the MBZ. We then derive the matrix form of the Hamil-
(±) tonian
functions ψk (r), containing contributions from the top
 
(φtk (r)) and bottom (φbk (r)) layers, with coefficients c±
t/b λk + ∆t (Gu − Gv ) ∆∗T (Gu − Gv )
H(k) =
obtained from DFT calculations: ∆T (Gu − Gv ) λk + ∆b (Gu − Gv )
(10)
with λk = −~2 (k − k0 + Gu )2 δuv /2m∗ . To derive the
ψk+ (r) = c+ b + t moiré electronic structure, we solve the eigenvalue prob-
b φk (r) + ct φk (r) (3)
lems of the effective Hamiltonian matrix in Eq. 10. The
ψk− (r) = c− b
b φk (r) + c− t
t φk (r) (4) MBZ reciprocal vectors Gm represent the truncated ba-
sis for the continuum expansion about k0 . We find that
The coefficients c±t,b satisfy the eigenvalue equation of lo- a |Gm | < 7|Gsc | is a sufficient number of Fourier compo-
cal electronic structures: nents to yield convergent eigenenergies. The ∆ potentials
    ±
∆b ∆ T c± ± cb
have been moved to Fourier space by a 2D fast Fourier
b = E . (5)
∆∗T ∆t c± c± transform (FFT), but the grid-sampling of the moiré po-
t t
tentials are first up-scaled to a 27 × 27 real space grid
Consequently, ∆t , ∆b and ∆T can be solved for at each (divisible by 3 to capture the AB/BA stacking points
displacement as follows: and to improve the accuracy) by using a harmonic ex-
pansion of the original 9 × 9 grid. The relaxations are
|c+ 2 + + 2 −
t | E + |cb | E |c+ 2 + + 2 −
b | E + |ct | E also computed in real space on a 27 × 27 grid, and the
∆t = , ∆ b =
|c+ 2 + 2
t | + |cb | |c+ 2 + 2
t | + |cb | updated potentials (given by ∆(d0 (r) + u(r))) are prop-
(6) erly symmetrized before applying the FFT. The Chern
 3

FIG. 2. Comparison of relaxation and electronic coefficients. (a) First order Fourier components of the GSFE energy c1 vs
anisotropic strain coefficient G. The relative intensity of intralayer and interlayer energy coefficients from stiff to soft is remarked
by the color from dark to bright. (b) Inverse effective mass vs the range of potential energy fluctuation in the MBZ. The relative
importance of energy terms is remarked by white and red color for kinetic energy dominant and potential energy dominant
condition, respectively.

TABLE I. Strain and GSFE coefficients and critical twist angles for atomic relaxation. The coefficients K and G are strain
modulus for isotropic and anisotropic distortion with the unit of eV per unit cell, respectively. c1 is 2 times of the 1st order
Fourier component of GSFE as defined in ref.[6]. The critical twist angle θc for evident relaxation effect is evaluated as equ. 11,
with regarding the critical twist when the separate domain walls begin to appear.
Materials G (eV/u.c.) K (eV/u.c.) c1 (meV/u.c.) θc (0 )
MoS2 27.36 45.45 14.25 1.7
MoSe2 28.52 43.17 13.54 1.7
MoTe2 22.16 35.63 20.18 2.2
WS2 35.27 52.33 14.93 1.6
WSe2 31.67 46.22 13.11 1.6
WTe2 23.61 33.82 17.67 2.1
hBN 40.66 59.65 4.82 0.8
SiH 19.33 31.47 1.13 0.6
Ars. 12.62 36.04 36.71 3.0
BP 25.14 43.23 25.67 2.3
Bi2 10.15 16.81 43.71 4.8
InSe 9.1 25.47 22.01 2.8
Gr. 47.78 67.87 4.9 0.8
GaS 24.11 38.9 13.53 1.7
SnS 20.62 34.1 11.04 1.7
MnBr2 13.79 20.41 10.68 2.1
ZnBr2 5.3 20.93 5.74 1.6
GaSe 22.98 36.67 13.2 1.8

number of each bands can then be computed as the in- laxation and electronic structures. The coefficients for
tegral of the Berry curvature through the 1st Brillouin mechanical properties related to atomic relaxation are
zone [9]. shownR in Fig I. Strain modulus K and G are fitted as
E = 21 K(uxx + uyy )2 + 12 G[(uxx − uyy )2 + 4u2xy ], repre-
senting the isotropic and anisotropic strain modulus, re-
RELAXATION AND ELECTRONIC spectively. The GSFE
COEFFICIENTS OF MATERIALS P energy i(ubis expanded as the Fourier
1 +vb2 )·r
series: GSFE(r)= u,v Fuv e . For materials
st
with hexagonal lattices, the 1 order components of Fij
Applying our algorithm to 18 different twisted bilayer are equal and we adopted two times of them as c1 to rep-
materials, we derive the coefficients for both atomic re-
 4

TABLE II. Coefficients of the band extrema properties. Most band extrema locate at high-symmetry k-points of K, Γ, and M.
If there are two competing band extrema, we denote K1 /K2 (K1 ) where K1 is the band extrema for untwisted structures, and
the following data is for K1 edge. K1 − K2 − K3 means the band extrema locate on the high symmetry line between the referred
high symmetry points. m∗ , V , and w denote the effective mass, the 1st -order Fourier coefficients of electrostatic potential and
tunneling matrix elements. θc is the critical twist angle for the electronic localization when the width of the uppermost band
is 10 meV.
Materials Band Extrema m∗ (me ) V (meV) w (meV) θc (0 )
MoS2 conduction K 0.54 5.64 4.8 1.95
valence Γ 2.6 12.01 95.8 4.11
MoSe2 conduction K 0.53 6.48 7.12 2.07
valence Γ 5.74 20.29 107.96 >6
MoTe2 conduction K 0.57 8.01 4.57 2.30
valence K 0.68 9.34 10.3 2.81
WS2 conduction K 0.34 2.02 1.52 1.48
valence Γ 3.21 2.67 64.22 4.85
WSe2 conduction K 0.35 9.18 8.75 1.70
valence K/Γ (K) 0.46 8.55 11.1 2.14
WTe2 conduction K 0.32 6.85 7.28 1.74
valence K 0.45 10.27 11.4 2.34
hBN conduction K/M (K) 0.69 99.21 134.85 3.48
valence K 0.52 45.42 58.53 2.55
SiH conduction Γ 0.23 1.35 3.5 1.34
valence Γ (degenerate) – – – –
InSe conduction Γ 0.32 15.3 138.37 1.70
valence Γ (degenerate) – – – –
BP conduction Γ 1.18 26.77 174.95 3.44
valence Γ 3.09 19.93 109.75 1.57
Ars. conduction Γ 2.01 37.55 249.52 >6
valence Γ 1.95 29.62 156.15 >6
GaS conduction Γ 0.29 11.78 67.93 1.37
valence K-Γ-M – – – –
GaSe conduction Γ 0.21 15.59 91.95 1.25
valence K-Γ-M – – – –

resent the GSFE intensity.[6] For materials with other als can be generally divided into three class: the materi-
lattices, we use 29 (max GSFE -min GSFE) for compari- als with high strain modulus and low GSFE coefficients
son, as this is consistent with the previous definition for (graphene, h-BN, SiH) are stiffer during relaxation, and
hexagonal lattices. The relative intensity of the GSFE their characteristic angles θc are about 1◦ ; the materials
and strain energy determines the characteristic twist an- with the intermediate proportion between the two coeffi-
gle when the in-plane relaxation effects play important cients (TMDCs, BP, etc.) have θc in the range of 1 − 3◦ ;
role. Large GSFE leads to strong in-plane relaxation, or while the materials with high GSFE coefficients and low
namely, large upper-limit of the twist angle where the strain modulus (InSe, Asenene, Bi) have low stiffness and
relaxation plays dominant role on the atomic structure. displays evident relaxation around or above 3◦ . This es-
According the analysis in previous work [10], the width timation is also well confirmed by the computed results
of the domain wall wd and critical twist angle θc are es- of relaxation.
timated as follow: The electronic properties of twisted bilayers are di-
r r vided into Moire molecule and Moire crystals according
a K 4 c1 to the degree of real space localization[6], which is de-
wd = , θc = √ (11)
4 c1 3 K termined by the relative importance of the kinetic and
potential energy in the MBZ. The comparison of the ef-
where a is the unit cell parameter (of hexagonal lattices).
fective mass and band edge energy fluctuation in config-
We set the criterion of θc as the condition when the do-
uration space are visualized in Fig 2b. As the kinetic
main walls overlap at AB/BA stacking.
energy scales with θ as T ∝ θ2 (m∗ )−1 , electronic √ local-
a
The energy coefficients of these materials are also plot- ization appears in the low θ region where θ < π~ m∗ V .
ted in Fig 2a to visualize their comparison. The materi- The critical twist angles for each materials when the up-
 5

FIG. 3. Moire electronic properties of TB MoS2 . Band structures of (a) 5.08◦ and (b) 5.08◦ MoS2 . The VBM and CBM derive
from Γ and K-edges, respectively. The red lines represent the results with full relaxation, and the black dash lines represent
the results without in-plane relaxation. (c) Realspace distribution of the MoS2 conduction and valance bands represented by
the local density of states (LDOS) with and without the in-plane relaxation. The LDOS are normalized to have the average
value of 1 in the MSL and plotted with the same color scale from 0 to 2. (d) Valence bands of twisted bilayer WSe2 . The red
and blue dots are the projection weight to the top and bottom layers. (e) Angle-dependence of moiré bands for twisted bilayer
h-BN. The format of this panel is the same as Fig. 3b in the main text. The small panel is the band structure at θ = 10 .

permost band width is about 10 meV are listed in the CONSISTENCY WITH PREVIOUS WORK
table. The h-BN, V group elementary substances and
Γ-edge TMDs typically have a high ratio of the poten- To validate our methods, we also do calculations that
tial to kinetic energy coefficients, indicating the evident could compare with the previous full DFT calculations,
localization appears for θ at about 3 − 4◦ ; while the com- as shown in Fig 3. The bands of 3.50 and 5.080 twisted
mon TMDs with K-edge and some other materials (like bilayer MoS2 in Fig 3a,b are consistent with the full DFT
InSe, SiH, ect.) display significant localization only for results in ref [12]. The real-space localization of the up-
small twist angle of about 1 − 2◦ . Otherwise, the prod- permost band states at AB/BA stacking in Fig 3c is
uct of the two coefficients determines the frequency ωθ in also the same with the full DFT results of the relaxed
the harmonic oscillator model[11], which determines the system. Valence bands of 5.080 twisted bilayer WSe2
slope of the average energy of the top bands to the twist shown in Fig 3d are accurately consistent with ref [13]
angle. regarding both the band energy and projection weight
to the top and bottom layer. This panel is also a good
example to show how the moiré bands originate from the
The complete data is listed in Tab II including the ML bands as illustrated in Fig 2b in the main text. The
band extrema and relative intensity of electrostatic po- angle-dependent bands of twisted bilayer h-BN are shown
tential, tunneling, and kinetic energy (represented by ef- in Fig 3e. We derive the consistent results of the pairing
fective mass). Specifically, the band extrema at Γ point behavior (single uppermost band and pairing 2nd and 3rd
typically host stronger tunneling and energy fluctuations, bands) and angle dependence with the moiré bands in ref
as well as a larger scale of bandgap variation during the [14].
inter-layer displacement and twisting.

BANDS REORDERING DURING THE

Our present calculations have not covered the degener- TOPOLOGICAL PHASE TRANSITION
ate band edge or band extrema out of the high symmetry
k-points. However, there is no substantial obstacle for Here we show the band reordering process during the
our methods to be generalized to different k-points. topological transition as illustrated in Fig 1d in the main
 6

(a)
(b)

X-type
wm*/θ2 (meV)

HC HC HC

HC TI TR
φ

Energy (meV)
Vm*/θ2 (meV)

Y/TR transition

X-type X/Y transition Y-type Y-type

X Y
TR TR TR TR TR

Angle φ

FIG. 4. Band reordering during the transition from TL to HL QHO conditions. (a) Calculated moire band structures for the
points on the 1/4 circle in the phase diagram with a radius of 4 meV and twist angle φ as shown in the left-top panel. The
angle φ gradually changes from 00 to 900 and the band structures change from HL, A-type, B-type, to TL condition through
the closing and reopening of the band gaps. (b) Evolution of the moire band maximum and minimum as a function of φ. The
1st , 2nd , 3rd , and 4th bands are denoted by black, red, blue, and purple lines, respectively. For each color, the upper line is the
band maximum and the lower line is the band minimum.

text. To present the change from V dominant to w dom-

inant moiré potential, we select a 1/4 circle route in the
phase diagram as shown in the first panel of Fig 4a. On
[1] G. Kresse and J. Furthmüller, Efficient iterative schemes
this route, the kinetic energy term keeps constant, and for ab initio total-energy calculations using a plane-wave
only the relative intensity of V and w changes, as rep- basis set, Phys. Rev. B 54, 11169 (1996).
resented by the variation of φ. The HL QHO bands [2] G. Kresse and D. Joubert, From ultrasoft pseudopoten-
appear at the small φ region. When tuning up φ, the tials to the projector augmented-wave method, Phys.
two top bands separate at K-point at about φ = 300 to Rev. B 59, 1758 (1999).
form the A-type topological phase. And then, the third [3] J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly con-
strained and appropriately normed semilocal density
band lifts and contact with the second band at Γ-point at
functional, Phys. Rev. Lett. 115, 036402 (2015).
φ = 510 , exchanging the Chern number by 2 and experi- [4] J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsin-
encing phase transition to B-type topological phase. The szky, H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu,
second band then contacts with the first band at Γ-point et al., Accurate first-principles structures and energies of
at φ = 560 and exchange the Chern number by 1, making diversely bonded systems from an efficient density func-
the Chern number of the uppermost band zero. The sys- tional, Nature chemistry 8, 831 (2016).
tem finally enters the TL phase where the topologically [5] H. Peng, Z.-H. Yang, J. P. Perdew, and J. Sun, Versa-
tile van der waals density functional based on a meta-
trivial uppermost bands separate from the lower bands.
generalized gradient approximation, Phys. Rev. X 6,
041005 (2016).
[6] S. Carr, D. Massatt, S. B. Torrisi, P. Cazeaux, M. Luskin,
This process is cohesively shown in Fig 4b, where the
and E. Kaxiras, Relaxation and domain formation in in-
maximum and minimum energy of each band are plot- commensurate two-dimensional heterostructures, Physi-
ted against the parameter φ. Ideal band degeneracy se- cal Review B 98, 224102 (2018).
quences of (1,2,3,...) and (2,4,6,...) display on the right [7] F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. MacDon-
and left side, respectively. Their reordering leads to the ald, Topological insulators in twisted transition metal
connection and separation in the intermediate region, dichalcogenide homobilayers, Physical review letters 122,
where the TI phases appear. 086402 (2019).
 7

[8] R. Bistritzer and A. H. MacDonald, Moiré bands in tons in moiré patterns of twisted bilayer transition metal
twisted double-layer graphene, Proceedings of the Na- dichalcogenides, Physical review letters 121, 266401
tional Academy of Sciences 108, 12233 (2011). (2018).
[9] X.-L. Qi and S.-C. Zhang, Topological insulators and su- [13] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes,
perconductors, Rev. Mod. Phys. 83, 1057 (2011). C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim,
[10] N. N. T. Nam and M. Koshino, Lattice relaxation and en- et al., Correlated electronic phases in twisted bilayer
ergy band modulation in twisted bilayer graphene, Phys. transition metal dichalcogenides, Nature Materials , 1
Rev. B 96, 075311 (2017). (2020).
[11] S. Carr, D. Massatt, M. Luskin, and E. Kaxiras, Du- [14] L. Xian, D. M. Kennes, N. Tancogne-Dejean,
ality between atomic configurations and bloch states M. Altarelli, and A. Rubio, Multiflat bands and strong
in twistronic materials, Phys. Rev. Research 2, 033162 correlations in twisted bilayer boron nitride: Doping-
(2020). induced correlated insulator and superconductor, Nano
[12] M. H. Naik and M. Jain, Ultraflatbands and shear soli- Letters 19, 4934 (2019).